4

这个问题是基于之前的一个类似问题。

我有以下等式和调整后的(一些随机数据):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 更准确,而不必仅对频率使用最小二乘法?

4

7 回答 7

5

我不认为 FFT 对(准)周期性信号的高分辨率频率测量有好处 - 见下文。

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

很多努力。

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

而是尝试测量波长。这可以像测量样本中过零之间的距离一样简单(可能需要多个周期以获得更高的精度 - 哎呀,如果你有那么多,测量 1000 个周期),然后将采样率除以得到频率。更简单、更快、更精确。

示例:48000 Hz 采样率、4.77 Hz 信号仅通过使用最粗略的方法测量一个周期的长度即可产生约 0.0005 Hz 的分辨率。(如果使用n个周期,频率分辨率也会乘以n。)

于 2010-01-21T21:16:40.557 回答
5

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

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)

周期图

Note that this function uses a logarithmic scale: plot(fr,10*log10(Pxx)) instead of plot(fr,Pxx)

于 2010-01-22T00:03:46.000 回答
1

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

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

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

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

于 2010-01-21T21:32:40.163 回答
1

What you are looking for is a frequency estimation method, and there are many. An FFT is one component of several estimation methods. Just using the peak magnitude bin, as in your example, gives you the worst resolution (but the greatest noise immunity to any other exactly periodic sinusoids). In low noise situations, you can interpolate. Parabolic interpolation of the log magnitude is one common estimator, but Sync interpolation of the FFT results may be better for a rectangular window. Zero-padding and doing a longer FFT is basically equivalent to interpolation.

对于零噪声的精确正弦曲线,忘记 FFT,只需求解 3 个未知数的方程,这可能只涉及 3 或 4 个非混叠样本点,此处此处的算法。

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

于 2012-04-17T22:45:55.900 回答
0

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

于 2010-01-21T21:19:12.207 回答
0

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

于 2010-01-21T21:28:59.793 回答
-1

试试开窗功能

于 2010-01-21T21:16:21.877 回答