2012-04-02 124 views
3

我試圖展示簡單正弦波的頻譜,因爲我們知道固定頻率的單個正弦波必須在其譜中的峯值上我寫這段代碼但我不能得到這個峯值什麼是錯在我的代碼:簡單正弦波在matlab中的傅立葉變換

clc 
nsteps=200;%number of signal elements in time domain 
i=sqrt(-1); 
NFREQS=100;%number of elements in frequency domain 
ddx=1e-9; 
dt=ddx/(6e8);%separation between each time domain elements 
lambdai=150e-9; 
lambdaf=500e-9; 
freqi=3e8/lambdai; 
freqf=3e8/lambdaf; 
freq=zeros(1,NFREQS); 
for j=1:NFREQS 
freq(j)=freqi-j*(freqi-freqf)/NFREQS;%desired frequency domain 
end 
arg=2*pi*freq*dt; 
et=zeros(nsteps,1); 
    for j=1:nsteps 
     et(j)=sin(2*pi*3e15*j*dt);%sin wave in time domain 
    end 
e=zeros(NFREQS,1); 
for n=1:NFREQS 
for j=1:nsteps 
e(n)=e(n)+et(j)*exp(-i*arg(n)*n);%sin wave in frequency domain 
end 
end 
lambda=linspace(lambdai,lambdaf,NFREQS); 
plot(lambda,abs(e)) 

The result is here

+0

什麼是你得到的結果呢? – 2012-04-02 14:49:21

+0

@OliCharlesworth我添加了編輯結果。 – peaceman 2012-04-02 14:57:24

回答

0

這裏的關鍵是使用窗口函數。將得到的代碼是這樣的:

clc 
clear 
nsteps=40000; 
i=sqrt(-1); 
Sc=0.5; 
ddx=1e-9; 
dt=ddx*Sc/(3e8); 
lambdai=100e-9; 
lambdaf=700e-9; 
lambda=lambdai:1e-9:lambdaf; 
[m,NFREQS]=size(lambda); 

    freq=3e8./lambda; 
et=zeros(nsteps,1); 
    for j=1:nsteps 
     et(j)=sin(2*pi*3e8/(300e-9)*dt*j);%sin wave in time domain 
    end 
    w=blackman(nsteps);%window function for periodic functions 
    for j=1:nsteps 
     et(j)=et(j)*w(j); 
    end 
    e=zeros(1,NFREQS); 
    for n=1:NFREQS 
     for j=1:nsteps 
      e(n)=e(n)+et(j)*exp(-i*2*pi*freq(n)*dt*j); 
     end 
    end 
    plot(lambda,abs(e)) 

,其結果是:

enter image description here

2

您可以考慮使用內置的傅立葉變換MATLAB提供編寫自己的,而不是。見http://www.mathworks.se/help/techdoc/math/brentm1-1.html

更新

大約有你應該知道的fft功能的一些特殊的東西。首先,它生成的數組需要歸一化因子1/N,其中N是數組的大小(這是因爲MATLAB函數的實現)。如果你不這樣做,頻域中的所有幅度將比它們「實際」大N倍。

第二件事是你需要以某種方式找到輸出數組中每個元素對應的頻率。下面代碼中的第三行將採樣時間數組轉換爲fourier數組對應的頻率。

該函數在時域中採用一個信號並將其繪製在頻域中。 t是時間數組,y是信號幅度。

function plotSpectrum(t, y) 
freqUnit = 1/(abs(t(2) - t(1)) * length(t)); 
f = -length(t) * freqUnit : freqUnit : (length(t) - 1) * freqUnit; 
oneSidedFFT = fft(y); 
fourier = horzcat(oneSidedFFT, oneSidedFFT); 
plot(f, abs(fourier)); 
xlabel('Frequency'); 
ylabel('Magnitude'); 
+0

當我在我的代碼中寫入fft(et)時,我遇到了MATLAB的fft函數的問題什麼是頻率範圍,我該如何改變這個範圍?你能否改變我的代碼來繪製傅里葉變換的「e」 .6e15到1e15 Hz)區域 – peaceman 2012-04-02 15:23:32