FFT of Gaussian filter in 1D

31 views (last 30 days)
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:

Accepted Answer

Matt J
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)
  5 Comments
Donya Khaledyan
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.

Sign in to comment.

More Answers (1)

Paul
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)
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.

Community Treasure Hunt

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

Start Hunting!