# plotting dirac delta function using inverse Fourier transform

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\$$.

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;


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

• 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. Commented May 2, 2020 at 6: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 Commented May 2, 2020 at 9:30
• @PG1995 No, just element-by-element multiplication of w, then multiplying by 0.0001. Sorry for the untidy code notation! Commented May 2, 2020 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! Commented May 3, 2020 at 3:09
• 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. Commented May 3, 2020 at 3:17