正弦波频率拟合答案

【问题标题】：Sine wave frequency fitting正弦波频率拟合
【发布时间】：2023-04-06 15:52:01
【问题描述】：

我有以下等式和调整后的（一些随机数据）：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 可以是found here。基本上它使 FFT 图居中并切断负信号。

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

【问题讨论】：

我认为问题在于您的预期频率（7hz）太接近您的频率下限（1hz，我假设），因此在频率维度（你将有不到三个 bits 留在那里）。我认为您需要将整个事物重新校准到一个允许很多低频域的域，例如大约 0.1 到 0.01hz。
一个猜测：你有没有简单地在 bin 索引上绘制频率？因为你的频率和时间绝对不匹配
零填充你的数据说256，以提高频率分辨率。

标签： math matlab fft octave data-fitting

【解决方案1】：

我认为FFT对于（准）周期性信号的微分辨率频率测量有益 - 见下文。

每个离散的FFT在非整数箱频率上扩展（即任何频率，该频率不适用于特定FFT的频率步骤之一）;这些“中间”频率将在最近的整数箱周围涂抹/展开。这种扩展的形状（“扩展功能”）取决于用于FFT的窗口函数。这种传播功能 - 简化和概括事情 - 要么非常狭窄，要么非常窄，但非常非常粗糙（非常高的山峰/非常低的山谷），或更宽但不那么粗糙。从理论上讲，您可以做一个非常精细的正弦波频率扫描并为每个频率计算FFT，然后通过将所有FFT的输出与导致该输出的频率一起保存所有FFT，您可以“校准”功能的形状和行为，然后通过比较要测量的信号的FFT输出到先前保存的结果并找到“最接近”的输出，找到更精确的频率。

很多努力。

但如果只需要测量单个信号的频率，则不执行此操作。

而是尝试测量波长。这可以像测量零十字之间的距离一样简单（可能是多个周期，以获得更多精度 - 如果您在样品中有许多），则测量1000个周期），并将采样率除以频率。更简单，更快，更精确。

示例：48000Hz采样率，4.77 Hz信号导致〜0.0005 Hz分辨率，仅通过测量具有粗糙方法的一个强>循环的长度。（如果使用 N EM>循环，则频率分辨率也乘以 N EM>。）

【讨论】：

我想到了这一点。唯一的问题是，频率可能在稍后在路上改变。但是，它看起来可能是我最好的选择，因为毕竟我将数据归一成到0居中。 span>
我没有在那里看到问题：如果频率会改变波长的变化，则测量新波长，执行划分 - 爆炸，结果：新频率。或者也许我不明白你想要做什么。 span>
零点是非常敏感的现实生活数据：即使是简单的传感器噪声也可以为您提供额外的零点 - 平滑是最简单的解决方案如果 i>你知道寻找什么频率。许多系统产生额外的频率分量（失真），易于电源嗡嗡声，或在显着级别捕获环境噪声。 span>
所有这一切都是真的，而且我也想要在op告诉更多关于问题后进入过滤，因为上面问题中所述的问题没有噪音。 span>
我只是猜测这是对某些或另一个方式的“真实”数据的定制的测试，其中一个或另一个方式 span>

【解决方案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's theorem 是否成立：

( 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)

【讨论】：

我就是这么干的，有什么不同？
就像我提到的，您可能错误地解释了信号频率。我的意思是从我从问题中了解到的，即使你随机噪声（以合理的方式）你仍然应该能够准确地恢复正弦的频率（尝试通过在开头添加：x = x + 0.5*randn(size(t));）

【解决方案3】：

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

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

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

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

【讨论】：

这里并没有真正涉及到 nyquist，因为我们可以在样本中看到几个完整的循环，并且预期的答案实际上低于生成的答案。
图上红色和绿色的波是实际数据，蓝色和黄色是FFT结果的最佳拟合。
@Phil：它定义了采样频率的下限，对于给出的结果来说“OK”，但仍然必须考虑。 -- @steven：我不确定发生了什么，正如其他人所指出的那样，无论如何你应该得到截然不同的值。

【解决方案4】：

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

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

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

【讨论】：

【解决方案5】：

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

【讨论】：

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

【解决方案6】：

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

【讨论】：

我认为您错误地读取了 FFT 最大的频率。

【解决方案7】：

试试windowing function？

【讨论】：

窗口函数如何解决这个问题？还有它与 FFT 和准确性有什么关系？
@steven：目前信号的边缘被简单地截断了，这会引入很多新的频率。应用窗口函数将平滑开始和结束，并根据函数删除/减少这些添加的频率。
@Jay：您想详细说明一下，以方便那些不知道窗户有什么作用的人吗？
@Phil H：在这个答案中链接的维基百科文章中讨论得比较好。我建议先咨询一下。
图上数据的结束位置就是数据的结束位置。这是要检查的数据范围。