3

I was trying to plot shifted dirac delta function in Matlab. $$\begin{align}\mathscr{F}\left(\delta(t-t_0)\right)&=\mathcal{F}(\omega)=e^{-j\omega t_0} \\ e^{-j\omega t_0}&=\cos\omega t_0-j\sin\omega t_0\end{align}$$

Using Trigonometric form of Fourier transform:

$$\begin{align}\mathrm{A}(\omega)-j\mathrm{B}(\omega)&=\cos\omega t_0-j\sin\omega t_0 \\ \delta(t-t_0)&=\frac 1{\pi}\int_0^{\infty}{\cos\omega t_0 \cos\omega t+\sin\omega t_0 \sin\omega t \space \mathrm{d}\omega}\end{align}$$

Let's assume \$t_0=\pi\$.

For more information on trigonometric form of FT: https://imagizer.imageshack.com/img922/8960/Z9xkMw.jpg

Given below is Matlab code and the plot generated using it. Where am I going wrong? Thank you for the help.

clear all; close all; clc;
t=linspace(-5,5,800);

for it=1:800

f=@(w)(1/pi).*(cos(w.*pi).*cos(w.*t(it))+ sin(w.*pi).*sin(w.*t(it)));
F(it)=integral(f,0,3000);

end

figure('Name','inverse Fourier transform');
plot(t,F,'red');
hold on;

enter image description here

a concerned citizen
  • 21,445
  • 1
  • 23
  • 40
PG1995
  • 245
  • 1
  • 10

1 Answers1

4

I don't think you did anything wrong. The problem is probably that a delta function is quite pathological, and the Matlab integral function can't handle things to sufficient numerical precision. Instead of a delta function, do a very narrow Gaussian pulse, which you can also Fourier transform analytically. Below is a shifted Gaussian pulse of width \$\Delta t=0.01\$

clear all; close all; clc;
t=linspace(-5,5,800);

for it=1:800

f=@(w)(1/pi).*(cos(w.*pi).*exp(-(w.*w*.0001/2)).*cos(w.*t(it))+ sin(w.*pi).*exp(-(w.*w*.0001/2)).*sin(w.*t(it)));
F(it)=integral(f,0,3000);

end

figure('Name','inverse Fourier transform');
plot(t,F,'red');
hold on;

Gives

enter image description here

rpm2718
  • 2,812
  • 1
  • 6
  • 14
  • Thank you for confirming that my work is correct. As you say it might only be limitations of Matlab's integral function. I'd try to code it in Mathematica though I have never coded in Mathematica before. I'd get back to you soon. – PG1995 May 02 '20 at 06:51
  • I wanted to clarify something. I used Wolframalpha to evluate Gaussian curve; I used second expression, g(x), from the top. I used this result for my code. For some reason I'm not getting expected result. Also, in you code, "w.w.0001". Are you take transpose of second "w" by putting "*" after it? Thanks – PG1995 May 02 '20 at 09:30
  • 1
    @PG1995 No, just element-by-element multiplication of w, then multiplying by 0.0001. Sorry for the untidy code notation! – rpm2718 May 02 '20 at 10:02
  • I re-run the code to see if I can make it working by changing values. There is no doubt that it has to do with the limitations of Matlab; more specifically 'integral' function. I was able to make it work and you can see a reasonable plot dirac delta here. In theory the dirac delta should have an infinite amplitude but its area/energy should be unity. I don't think I'd try to code it into Mathematica because, as you indicated, my work was correct. t=linspace(0,5,1400); for it=[1:1400]; F(it)=integral(f,0,6000); Thanks! – PG1995 May 03 '20 at 03:09
  • By the way, it was a trial and error method. If you change the numbers, the 'reasonable plot' of dirac delta might not appear to be that reasonable. – PG1995 May 03 '20 at 03:14
  • 1
    I accidentally discovered that the Matlab integral function was using only 150 points to calculate the integral function F.....way too few for the highly oscillatory nature of the function it was integrating. – rpm2718 May 03 '20 at 03:17