matlab - Matlab 计算频谱的方法

Question

我在 Matlab 中计算时间序列的频谱时遇到一个问题。我已阅读有关“fft”函数的文档。然而，我看到了两种实现方式，并且都给了我不同的结果。我很高兴能得到一些关于这种差异的答案:

第一种方法:

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);

第二种方法:

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)));

我认为Matlab中的“fft”函数中没有必要使用“nextpow”函数。最后，哪个好？

谢谢

Answer 1

简短回答:您需要windowing用于频谱分析。

现在给出长答案...在第二种方法中，您使用的是优化的 FFT 算法，当输入向量的长度是 2 的幂时，该算法非常有用。假设您的原始信号有 401 个来自无限长信号的样本(如下例所示)； nextpow2() 将为您提供 NFFT=512 个样本。当您将较短的 401 样本信号输入到 fft() 函数时，它会隐式补零以匹配请求的 512 长度 (NFFT)。但是(棘手的部分来了):对信号进行零填充相当于将无限长的信号乘以 rectangular function ，在频域中转换为具有 sinc function 的卷积的操作。这就是半对数图底部本底噪声增加的原因。

避免这种噪声增加的方法是使用平滑器 window function 手动创建要输入到 fft() 的 512 个样本信号。而不是默认的矩形。加窗意味着将信号乘以锥形、对称的信号。关于选择良好的窗函数有大量文献，但具有低旁瓣(低噪声增加)的典型准确窗函数是 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' )