正弦波频率拟合

Question

这个问题基于一个类似的问题：previous similar 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)

正向快速傅里叶变换（Positive FFT）可以在这里找到。它基本上将FFT图形居中并切断负信号。

我的问题是，如何在不必仅针对频率使用最小二乘法的情况下使FFT更加准确？

- monksy

---

我认为问题在于您期望的频率（7赫兹）与下限频率（我假设是1赫兹）太接近，因此在频率维度上失去了大部分的表征精度（在那里你只剩下不到三个“比特”）。我认为您需要重新校准整个系统，以适应更低的频率范围，例如0.1至0.01赫兹左右。 - RBarryYoung

猜测一下：您是否只是在 bin 索引上简单绘制了频率？因为您的频率和时间绝对不匹配。 - peterchen

将您的数据填充为256，以增加频率分辨率。 - Paul

Answer 1

正如其他人所提到的，您误解了信号的频率。让我举个例子来澄清一些事情：

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's theorem是否成立：

( sum(x.^2) - sum(X) )/nfft < 1e-6

---

选项二

或者，我们可以使用信号处理工具箱函数：

%# 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)。

- Amro

---

这基本上就是我所做的，有什么不同吗？ - monksy

就像我之前提到的，你可能错误地解释了信号频率。我的意思是，从问题中我所理解的来看，即使你以合理的方式添加随机噪声，你仍然应该能够精确地恢复正弦波的频率（在开头尝试添加：x = x + 0.5*randn(size(t));）。 - Amro

Answer 2

我认为FFT不适合用于(quasi)periodic信号的精细分辨率频率测量 - 请参见下文。

每个离散FFT都会在非整数频率上扩展(即在任何不完全对应于特定FFT频率步骤之一的频率上); 这些“中间”频率将在最近的整数bin周围被涂抹/扩展。这种扩展的形状（“扩展函数”）取决于用于FFT的窗口函数。简化和概括事情，这个扩展函数要么非常窄但非常不平滑（非常高的峰值/非常低的谷值），要么更宽但不那么不平滑。理论上，您可以对正弦波进行非常精细的频率扫描，并为每个正弦波计算FFT，然后通过保存所有FFT的输出以及导致该输出的频率来“校准”函数的形状和行为，然后通过比较要测量的信号的FFT输出与之前保存的结果，并找到“最接近”的结果以找到更精确的频率。

需要大量努力。

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

相反，尝试测量波长。这可以简单地通过测量样本中的零交叉点之间的距离（也许对于多个周期以获得更高的精度 - 如果您有那么多，甚至可以测量1000个周期），并将采样率除以该距离来到达频率。更简单、更快速且更精确。

例如：48000 Hz采样率，4.77 Hz信号通过使用最简单的方法测量一个周期的长度就能得到大约0.0005 Hz的分辨率。（如果您取n个周期，则频率分辨率也会乘以n。）

- Bandi-T

---

我考虑过那样做。唯一的问题是频率以后可能会改变。但是，看起来这可能是我最好的选择，因为毕竟我已经将数据归一化为0中心化。 - monksy

我在这里没有看到问题：如果频率改变，波长也会随之改变，你测量新的波长，做除法 - 就可以得到新的频率。或者我可能不明白你想做什么。 - Bandi-T

零交叉对现实数据非常敏感：即使是简单的传感器噪声也可能会给你额外的零交叉 - 如果您知道要查找的频率，平滑是最简单的解决方案。许多系统会产生额外的频率成分（失真），容易受到电源杂音的影响，或者捕获显着水平的环境噪声。 - peterchen

这一切都是真的，当 OP 更多地讲述问题后，我也想进入过滤器中，因为在上面的问题陈述中，问题没有噪音。 - Bandi-T

我猜这是一个测试，他最终以某种方式得到“真实”的数据。 - peterchen

图表中的红色信号几乎没有噪音（大部分已被归一化），虽然不足以影响结果。 - monksy

Answer 3

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

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

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

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

- peterchen

---

这里实际上与奈奎斯特没有关系，因为我们可以在样本中看到几个完整的周期，并且预期的答案实际上比生成的答案要低。 - Phil H

图表上的红色和绿色波形是实际数据，蓝色和黄色是来自FFT结果的最佳拟合。 - monksy

@Phil：它定义了样本频率的下限，对于给出的结果来说“还好”，但仍需考虑。-- @steven：我不确定发生了什么，正如其他人指出的那样，你应该得到完全不同的值。 - peterchen

Answer 4

你所需要的是一种频率估计方法，而这种方法有很多。FFT是几种估计方法中的一个组成部分。就像你的例子中只使用峰值幅度二进制一样，这会给你最差的分辨率（但对于任何精确周期正弦波具有最大的噪声免疫力）。在低噪声情况下，你可以进行插值。对数幅度的抛物线插值是一种常见的估计器，但对于矩形窗口，FFT结果的同步插值可能更好。零填充并进行更长的FFT基本上相当于插值。

对于零噪声中的精确正弦波，忘记FFT，只需解决三个未知数的方程，这可能涉及到至少3或4个非混叠采样点，做这件事的算法在这里和在这里。

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

- hotpaw2

Answer 5

如果你是从一个函数中生成数据，而不是使用样本数据，那么你可以生成大量的点并运行一个大的FFT，这样高精度下频率分辨率会非常小。但这并不能解决基本问题。

- Joe Koberg

---

样本数据[如图所示]是有限的点，但“最佳拟合数据”是由一个函数生成的。 - monksy

Answer 6

首先，对你的问题进行更正：(30/2PI)不是频率。你的信号频率是1/30乘以你使用的采样率。 其次，请问您的sampledata向量长度是多少？当FFT返回一个值向量时，第i个值将对应于f_i = i/N，其中N是向量的长度，i∈[0,N-1]。 您希望i/N恰好等于1/30，对于某些整数i。换句话说，N应该等于30*i，即N应该是30的倍数。现在，您使用的向量长度是否是30的倍数？如果不是，请尝试使它成为30的倍数，这应该可以解决问题。

- morpheus

---

我认为你可能错误地读取了FFT最大值的频率。 - morpheus

Answer 7

尝试一下窗口函数？

- Jay Kominek

---

窗函数如何解决这个问题？它与FFT和精度有什么关系？ - monksy

@Jay：你想详细说明一下窗口对于那些不知道其作用的人有什么好处吗？ - Phil H

@Phil H：在这个答案中链接的维基百科文章中已经相对详细地讨论了这个问题。我建议先参考那篇文章。 - Bandi-T

由于信号的峰值可能不会完全在一个bin中，我认为振铃可能会在相邻的bin中累加，将峰值向左拉了一点。（尽管与峰值相比微不足道。）此外，如果Steven正在使用Matlab，则尝试应该只需要大约30秒钟。 http://www.mathworks.com/access/helpdesk/help/toolbox/signal/f9-131178c.html#f9-140546 - Jay Kominek

-1：窗口技术有助于消除相邻干扰和频谱“扩散”，但在零和非常低噪声情况下，会使峰值变得更宽，从而降低频率估计的准确性。 - hotpaw2