2015-12-19 4 views
3

Я пытаюсь получить gaussian curve с помощью mft. Проблема в том, что в одном случае мой attemp для уменьшения шума путем деления F=dt.*F./exp(-1i.*nu.*T/2) не работает (img 1), а во втором случае, если я пытаюсь принять абсолютное значение результата fft, у меня нет подходящего масштаба в графе (img 2).Matlab FFT для функции gaussian

N=512; 
T=10; 
dt=T/(N-1); 
t=linspace(-5,5,N); 
f=exp(-t.^2); 
F=fft(f); 

F1=F(1:N/2+1); 
F2=F(N/2+1:N); 
F=[F2,F1]; 

dnu=(N-1)/(N*T); 
nuNyq=1/(2*dt); 
nu=-nuNyq+dnu*(0:N); 
F=dt.*F; 
%F=dt.*F./exp(-1i.*nu.*T/2); 


y=linspace(-5,5,N); 
F2=pi.^(1/2).*exp(-y.^2/4); 

hold on 
plot(y,F2); 
%plot(nu,real(F),'r'); 
plot(nu,abs(F),'r'); 
legend('analiticFT','FFT') 
xlim([-5 5]) 
hold off 

IMG 1 img. 1

img2 img. 2

+1

Просто к сведению, это гораздо проще, чтобы помочь вам, когда вы даете переменные значимые имена. (Примечание 2: для достижения того, что вы делаете с F1 и F2, вы можете вызвать 'fftshift (F)') – dangom

+0

Что вы имеете в виду с _is не работает? _ –

ответ

3

кажется масштабирование в формуле для аналитического преобразования Фурье не вполне корректно. Согласно this Fourier Transform table on Wikipedia, преобразование непрерывного сигнала во временной области

y\left(t\right) = e^{-a t^2}

является

Y\left(f\right) = \sqrt{\frac{\pi}{a}} \cdot e^{-\frac{\left(\pi f\right)^2}{a}}

где в вашем случае a=1. Соответственно, вы должны сравнить БПФ сигнала временной области

t=linspace(-5,5,N); 
f=exp(-t.^2); 

с аналитическим преобразованием Фурье

F2 = sqrt(pi)*exp(-(pi*y).^2); 

Таким образом, откладывая сравнение с:

hold off; 
plot(y,F2); 
hold on; 
plot(nu,abs(F),'r'); 
legend('analiticFT','FFT') 
xlim([-5 5]) 

урожайности:

enter image description here

Теперь, когда мы установили правильную основу для сравнения, мы можем посмотреть, почему вы получаете колебания в img 1. Проще говоря, опорный гауссовский импульс f=exp(-t.^2);, который вы создали, имеет пик при t=0. Соответствующий «нулевой» дискретный момент времени, естественно, является первым индексом в массиве (индекс 1). Однако в вашем массиве этот пик отображается с индексом N/2. Под номером Shift theorem это вызывает дополнительный exp(-pi*j*k) член в частотной области, отвечающий за колебания, которые вы видите. Чтобы исправить это, вы должны перейти обратно ваш гауссов импульс с ifftshift:

F=fftshift(fft(ifftshift(f))); 
Смежные вопросы