Borrar filtros
Borrar filtros

FFT of Gaussian filter in 1D

32 visualizaciones (últimos 30 días)
2NOR_Kh
2NOR_Kh el 10 de Abr. de 2022
Comentada: 2NOR_Kh el 11 de Abr. de 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:

Respuesta aceptada

Matt J
Matt J el 10 de Abr. de 2022
Editada: Matt J el 11 de Abr. de 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 comentarios
Matt J
Matt J el 11 de Abr. de 2022
I'm not sure if that's true. You haven't provided enough to allow us to run your code.
2NOR_Kh
2NOR_Kh el 11 de Abr. de 2022
My time axis is a parameter of a large dataset, however, your code fixed my problem and I really appreciate your help.

Iniciar sesión para comentar.

Más respuestas (1)

Paul
Paul el 10 de Abr. de 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.

Productos


Versión

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by