# FFT of Gaussian filter in 1D

31 views (last 30 days)
Donya Khaledyan on 10 Apr 2022
Commented: Donya Khaledyan on 11 Apr 2022
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:

Matt J on 10 Apr 2022
Edited: Matt J on 11 Apr 2022
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)
Donya Khaledyan on 11 Apr 2022
My time axis is a parameter of a large dataset, however, your code fixed my problem and I really appreciate your help.

Paul on 10 Apr 2022
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)
figure
hold on
fplot(angle(G(w,sigmaval)),[dftfreq(1) dftfreq(end)],'-x')
stem(dftfreq,angle(dt*Gdft),'r')
xlim([-1e7 1e7]/2)
The angle of the DFT is zero or pi depending on the numerical noise in the imaginary part of the Gdft.

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!