正弦波频率拟合

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 的输出以及导致该输出的频率来“校准”函数的形状和行为，然后通过将待测信号的 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)

Answer 3

假设 N 是以秒为单位的时间，则频率为 1/30Hz (y=A * sin( 2* PI * f * t))

频率分辨率 = 采样率 / FFT 点数

采样率由下式确定根据奈奎斯特标准，采样率（样本/秒）必须至少是要分析的最大频率的两倍，例如 48kHz 用于分析高达 24kHz 的频率。 （对于“现实生活”数据，最好有一点缓冲区）。

因此，您可能需要增加 FFT 的大小。

Answer 4

您正在寻找的是一种频率估计方法，并且有很多方法。 FFT 是多种估计方法的组成部分。仅使用峰值幅度箱（如您的示例中所示）会给您带来最差的分辨率（但对任何其他完全周期性的正弦曲线具有最大的抗噪能力）。在低噪声情况下，您可以进行插值。对数幅度的抛物线插值是一种常见的估计器，但 FFT 结果的同步插值对于矩形窗口可能更好。补零并进行更长的 FFT 基本上相当于插值。

对于零噪声中的精确正弦曲线，忘记 FFT，只需求解 3 个未知数的方程，这可能涉及少至 3 或 4 个非混叠样本点，执行此操作的算法 此处 和 在这里。

我在 DSP 网页 上列出了一些其他频率估计方法。

Answer 5

如果您从函数生成，而不是使用样本，则可以生成很多点并运行大 fft，因此频率仓非常小，可以实现高精度。但这并不能解决根本问题。

Answer 6

首先，更正您的问题：（30/2PI）不是频率。信号的频率是 1/30 * 无论您使用什么采样率。
其次，你能告诉我样本数据向量的长度是多少吗？当 FFT 返回值向量时，第 i 个值将对应于 f_i = i/N，其中 N 是向量长度，i \in [0,N-1]
对于某个整数 i，您希望 i/N 恰好等于 1/30。换句话说，N 应该等于 30*i，即 N 应该是 30 的倍数。 现在，您使用的向量的长度是 30 的倍数吗？ 如果不是，请尝试制作它，并且那应该可以解决问题。

Answer 7

尝试一下窗口函数？