Как уже упоминалось, вы неправильно интерпретируете частоту сигнала. Позвольте мне привести пример, чтобы прояснить несколько вещей:
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')
Теперь вы можете видеть, что пик в БПФ соответствует исходной частоте сигнала при 6 Гц
[v idx] = max(X);
fr(idx)
ans =
6
Мы можем даже проверить, что Теорема Парсеваля справедлива:
( 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)
Обратите внимание, что эта функция использует логарифмическую шкалу: plot(fr,10*log10(Pxx)) вместо plot(fr,Pxx)