正如其他人所提到的,您误解了信号的频率。让我举一个例子来澄清一些事情:
Fs = 200; %# sampling rate
t = 0:1/Fs:1-1/Fs; %# time vector of 1 second
f = 6; %# frequency of signal
x = 0.44*sin(2*pi*f*t); %# sine wave
N = length(x); %# length of signal
nfft = N; %# n-point DFT, by default nfft=length(x)
%# (Note: it is faster if nfft is a power of 2)
X = abs(fft(x,nfft)).^2 / nfft; %# square of the magnitude of FFT
cutOff = ceil((nfft+1)/2); %# nyquist frequency
X = X(1:cutOff); %# FFT is symmetric, take first half
X(2:end -1) = 2 * X(2:end -1); %# compensate for the energy of the other half
fr = (0:cutOff-1)*Fs/nfft; %# frequency vector
subplot(211), plot(t, x)
title('Signal (Time Domain)')
xlabel('Time (sec)'), ylabel('Amplitude')
subplot(212), stem(fr, X)
title('Power Spectrum (Frequency Domain)')
xlabel('Frequency (Hz)'), ylabel('Power')
现在您可以看到 FFT 中的峰值对应于 6Hz 处信号的原始频率
[v idx] = max(X);
fr(idx)
ans =
6
我们甚至可以检查Parseval 定理是否成立:
( sum(x.^2) - sum(X) )/nfft < 1e-6
选项 2
或者,我们可以使用信号处理工具箱功能:
%# estimate the power spectral density (PSD) using the periodogram
h = spectrum.periodogram;
hopts = psdopts(h);
set(hopts, 'Fs',Fs, 'NFFT',nfft, 'SpectrumType','onesided')
hpsd = psd(h, x, hopts);
figure, plot(hpsd)
Pxx = hpsd.Data;
fr = hpsd.Frequencies;
[v idx]= max(Pxx)
fr(idx)
avgpower(hpsd)
Note that this function uses a logarithmic scale: plot(fr,10*log10(Pxx))
instead of plot(fr,Pxx)