FFT of Gaussian filter in 1D

I have a gaussian filter and its fourier transform in matlab is just located at zero. I am sure I did everything right but still dont know what is the problem, maybe it is the right fft of the signal but I cant explain why.

This is my code:

sigma = 3.76e-6;
u=0;
t1=(-max(t(:)):t(2):max(t(:)));
Exp_comp = -((t1-u).^2)/(2*sigma*sigma);%tn/(sigma*sqrt(2))
G = exp(Exp_comp)/((sqrt(2*pi))*sigma);
plot(t1,G);title('G');
FG= abs(fftshift(fft(G)));
freq=1/t(2)*(-6783:6783)/13567;
figure;plot(freq/1e6,FGHE2);title('FGHE2');

This is the plot of my gaussian and fourier of it:

I zoomed the fourier signal in to take this picture:

Answer 1

[Matt J]

sigma = 3.76e-6;
u=0;

%sampling parameters
dt=sigma/30;5         %time sampling interval
df=1/sigma/30;        %frequency sampling interval
N=ceil(1/df/dt);      %number of samples
df=1/N/dt;            %adjust df so that N=1/(df*dt)

%sampling axes
nAxis=(0:N-1)-ceil((N-1)/2);
tAxis=nAxis*dt;        %time axis sample locations
fAxis=nAxis*df;        %frequency axis sample locations

%signal samples
Exp_comp = -(tAxis-u).^2/(2*sigma^2);
G = exp(Exp_comp)/((sqrt(2*pi))*sigma);
FG=abs( fftshift( (fft(ifftshift(G*dt))) ) );

%plots 
figure
plot(tAxis,G);title('G'); xlim([-1,1]*4*sigma)

figure;
plot(fAxis,FG);title('FGHE2'); xlim([-1,1]/sigma)

[2NOR_Kh]

Matt: Thank you, what are those dt and df, and why we need them?

[Matt J]

dt is the time sampling interval and df is the frequency sampling interval. You should decide on those first. I just picked some values that seemed reasonable in relation to sigma.

[2NOR_Kh]

I put my time matrix in the code, everything exactly the same as the code I put here and you know what? that xlim([-1,1]/sigma) part made everything work. why and how?

[Matt J]

I'm not sure if that's true. You haven't provided enough to allow us to run your code.

[2NOR_Kh]

My time axis is a parameter of a large dataset, however, your code fixed my problem and I really appreciate your help.

Answer 2

[Paul]

Hello @Donya Khaledyan

The code in the question doesn't define t nor FGHE2, so the plots can't be recreated. Code below might be helpful

First define the symbolic solution

syms t w
syms sigma positive
g(t,sigma) = exp(-t^2/2/sigma/sigma)/sqrt(2*sym(pi))/sigma;
G(w,sigma) = simplify(fourier(g(t,sigma),t,w));

Define g as a function to evaluate numerically.

gfunc = matlabFunction(g(t,sigma),'Vars',{'t','sigma'});

Define the value of a sigma, the time vector for sampling, and the values of g(t)

sigmaval = 3.76e-6;
dt = 1e-7;
tvals = -3e-5:dt:3e-5;
gvals = gfunc(tvals,sigmaval);

Plot the symbolic g(t) and the samples

figure;
hold on
fplot(g(t,sigmaval),[tvals(1) tvals(end)]);
plot(tvals,gvals,'o')

Find the shift that puts t(1) at t = 0

nshift = round(-tvals(1)/dt)

nshift = 300

Compute the DFT the sequence of gvals. The multiplication by exp() is only needed to get the correct phase.

dftfreq = (0:numel(tvals)-1)/numel(tvals)*2*pi; 
Gdft = fft(gvals).*exp(1j*nshift*dftfreq);

Shift the DFT to the negative frequencies and get the corresponding frequency vector

Gdft = fftshift(Gdft);
n = numel(tvals);
dftfreq = ceil(-n/2:(n/2-1))/(n-1)*2*pi;
dftfreq = dftfreq/dt;  % convert to rad/sec

Plot the CTFT and the (scaled) DFT

figure;
hold on
fplot(abs(G(w,sigmaval)),[dftfreq(1) dftfreq(end)])
ylim([0 1])
stem(dftfreq,abs(dt*Gdft),'r')
xlim([-1e7 1e7]/2)
xlabel('Freq (rad/sec)'); ylabel('Magnitude')

figure
hold on
fplot(angle(G(w,sigmaval)),[dftfreq(1) dftfreq(end)],'-x')
stem(dftfreq,angle(dt*Gdft),'r')
xlim([-1e7 1e7]/2)
xlabel('Freq (rad/sec)');ylabel('Angle (rad)')

The angle of the DFT is zero or pi depending on the numerical noise in the imaginary part of the Gdft.