計算譜Matlab的方法

Question

在Matlab中計算時間序列的譜圖時，我有個問題。我已閱讀有關'fft'功能的文件。但是我已經看到了兩種實現方式，並且都給了我不同的結果。我很感激有這個差一些答案：

第1種方法：

nPoints=length(timeSeries);  
Time specifications: 
Fs = 1; % samples per second 
Fs = 50; 
freq = 0:nPoints-1; %Numerators of frequency series 
freq = freq.*Fs./nPoints; 
% Fourier Transform: 
X = fft(timeSeries)/nPoints; % normalize the data 
% find find nuquist frequency 
cutOff = ceil(nPoints./2); 
% take only the first half of the spectrum 
X = abs(X(1:cutOff)); 
% Frequency specifications: 
freq = freq(1:cutOff); 
%Plot spectrum 
semilogy(handles.plotLoadSeries,freq,X); 

第2種方法：

NFFT = 2^nextpow2(nPoints); % Next power of 2 from length of y 
Y = fft(timeSeries,NFFT)/nPoints; 
f = 1/2*linspace(0,1,NFFT/2+1); 
% % Plot single-sided amplitude spectrum. 
% plot(handles.plotLoadSeries, f,2*abs(Y(1:NFFT/2+1)))  
semilogy(handles.plotLoadSeries,f,2*abs(Y(1:NFFT/2+1))); 

我認爲這是沒有必要使用'nextpow'功能在' fft'函數。最後，哪個好？

感謝

Answer 1

簡短的回答：你需要windowing頻譜分析。

現在爲長時間的答案...在第二種方法中，當輸入向量的長度是2的冪時，使用優化的FFT算法很有用。假設您的原始信號具有401個樣本（如下面的示例所示）來自無限長的信號; nextpow2()會給你NFFT = 512個樣本。將較短的401採樣信號送入fft()函數時，它將被隱式填充爲零，以匹配請求的512（NFFT）長度。但是（這裏出現了棘手的部分）：對信號進行零填充相當於將無限長信號乘以rectangular function，該操作在頻域中轉化爲與sinc function的卷積。這將是在半對數圖底部增加背景噪音的原因。

避免噪聲增加的一種方法是手動創建要採用fft()的512採樣信號，使用平滑的window function而不是默認的矩形信號。開窗意味着將信號乘以錐形，對稱的信號。關於選擇一個好的開窗功能，有大量的文獻，但是低旁瓣（低噪聲增加）的典型精確的文獻是Hamming function，在MATLAB中被實現爲hamming()。

這裏是說明該問題的圖（在頻域和時域）：

...，並且生成該圖中，代碼：

clear 

% Create signal 
fs = 40;   % sampling freq. 
Ts = 1/fs;   % sampling period 
t = 0:Ts:10;  % time vector 
s = sin(2*pi*3*t); % original signal 
N = length(s); 

% FFT (length not power of 2) 
S = abs(fft(s)/N); 
freq = fs*(0:N-1)/N; 

% FFT (length power of 2) 
N2 = 2^nextpow2(N); 
S2 = abs(fft(s, N2)/N2); 
freq2 = fs*(0:N2-1)/N2; 
t2 = (0:N2-1)*Ts;  % longer time vector 
s2 = [s,zeros(1,N2-N)]; % signal that was implicitly created for this FFT 

% FFT (windowing before FFT) 
s3 = [s.*hamming(N).',zeros(1,N2-N)]; 
S3 = abs(fft(s3, N2)/N2); 

% Frequency-domain plot 
figure(1) 
subplot(211) 
cla 
semilogy(freq,S); 
hold on 
semilogy(freq2,S2,'r'); 
semilogy(freq2,S3,'g'); 
xlabel('Frequency [Hz]') 
ylabel('FFT') 
grid on 
legend('FFT[401]', 'FFT[512]', 'FFT[512] with windowing') 

% Time-domain plot 
subplot(212) 
cla 
plot(s) 
hold on 
plot(s3,'g') 
xlabel('Index') 
ylabel('Amplitude') 
grid on 
legend('Original samples', 'Windowed samples')