FFT of Gaussian filter in 1D
58 次查看(过去 30 天)
显示 更早的评论
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:
0 个评论
采纳的回答
Matt J
2022-4-10
编辑:Matt J
2022-4-11
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 个评论
Matt J
2022-4-11
I'm not sure if that's true. You haven't provided enough to allow us to run your code.
更多回答(1 个)
Paul
2022-4-10
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)
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.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Bartlett 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!