正弦波频率拟合

Question

此问题基于之前的类似问题.

我有以下等式和调整后的(一些随机数据):0.44*sin(N*2*PI/30)

我试图使用FFT从生成的数据中获取频率.然而,频率最终接近但不等于频率(这使得波比预期的大一点)

FFT的最大频率为7hz,但预期频率为(30/2PI)4.77hz.

我已经包含了FFT和绘图值的图表.

我使用的代码是:

[sampleFFTValues sFreq] = positiveFFT(sampledata, 1);
sampleFFTValues = abs(sampleFFTValues);
[v sFFTV]= max(sampleFFTValues)

可在此处找到正FFT .基本上它将FFT图中心并切断负信号.

我的问题是如何才能使FFT更准确,而不必仅针对频率采用最小二乘法？

Answer 1

我认为 FFT 不适合（准）周期信号的精细分辨率频率测量 - 见下文。

每个离散 FFT 都在非整数 bin 频率上进行扩展（即在与特定 FFT 的频率步长之一不完全对应的任何频率上）；这些“中间”频率将被涂抹/散布在最近的整数箱周围。这种扩展的形状（“扩展函数”）取决于用于 FFT 的窗口函数。这种扩展函数 - 简化和概括事物 - 要么非常窄但非常参差不齐（非常高的峰/非常低的谷），或者更宽但不那么参差不齐。从理论上讲，您可以对正弦波进行非常精细的频率扫描并为它们中的每一个计算 FFT，然后您可以通过将所有 FFT 的输出与导致该输出的频率一起保存来“校准”函数的形状和行为，

很多努力。

但是如果您只需要测量单个信号的频率，请不要这样做。

而是尝试测量波长。这可以像测量样本中的零交叉之间的距离一样简单（可能需要多个周期以获得更高的精度 - 哎呀，如果您有这么多周期，则测量 1000 个周期），然后将采样率除以得出频率。更简单、更快、更精确。

示例：48000 Hz 采样率，4.77 Hz 信号仅通过使用最粗略的方法测量一个周期的长度就可以产生 ~0.0005 Hz 的分辨率。（如果采用n个周期，频率分辨率也会乘以n。）

Answer 2

正如其他人所说,你错误地解释了信号的频率.让我举一个例子来澄清一些事情:

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)

请注意,此函数使用对数刻度:plot(fr,10*log10(Pxx))而不是plot(fr,Pxx)