3
\$\begingroup\$

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

\$\endgroup\$

1 Answer 1

4
\$\begingroup\$

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

\$\endgroup\$
6
  • \$\begingroup\$ 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. \$\endgroup\$
    – PG1995
    May 2, 2020 at 6:51
  • \$\begingroup\$ 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 \$\endgroup\$
    – PG1995
    May 2, 2020 at 9:30
  • 1
    \$\begingroup\$ @PG1995 No, just element-by-element multiplication of w, then multiplying by 0.0001. Sorry for the untidy code notation! \$\endgroup\$
    – rpm2718
    May 2, 2020 at 10:02
  • \$\begingroup\$ 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! \$\endgroup\$
    – PG1995
    May 3, 2020 at 3:09
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – rpm2718
    May 3, 2020 at 3:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.