-
22-09-2019 - |
题
这个问题是基于一个 以前的类似问题。
我有以下方程和调整后的(一些随机数据):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 更准确,而不必仅采用最小二乘法来计算频率?
解决方案
我不认为FFT是用于精细分辨率频率测量(准)周期性信号良好 - 见下文
每个离散FFT已经在非整数仓频率扩展(即在其上不完全对应的特定FFT的频率的步骤之一的任何频率);这些“中间”的频率将被涂抹/展开围绕最接近的整数箱。此扩展(“扩展功能”)的形状取决于用于FFT窗函数。此扩展函数 - 来简化和泛化的事情 - 要么非常窄,但非常-非常参差不齐(非常高的峰值/非常低的谷),或更宽但小于粗糙。从理论上讲,你可以通过保存所有的FFT输出与导致该输出频率一起做正弦波,并计算FFT的一个非常精细的频率扫描为他们每个人,然后你就可以“校准”功能的形状和行为,然后通过比较信号的FFT输出被测定为先前保存的结果并发现“最接近”一个找到一个更精确的频率。
的大量努力。
但不这样做,如果你只需要测量的单个信号的频率。
相反尝试测量波长。这可以简单到作为测量零个交叉之间的距离(可能为多个循环以获得更多的精度 - 赫克,测量1000次循环,如果你有许多)中的样品,并通过将样本率在频率到达。更简单,更快速的和更精确。
例:48000赫兹的采样速率,4.77赫兹信号导致〜0.0005 Hz分辨率只是通过测量的一个强>周期的长度与最原始的方法。 (如果拿名词的周期,频率分辨率由乘法名词的为好。)
其他提示
如所提到的由他人,你曲解信号的频率。让我举一个例子来清除几件事情:
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')
现在可以看到在6HZ
在FFT对应于信号的原始频率的峰值[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)
假设N是时间(秒),您的频率为1 / 30Hz的(y=A * sin( 2* PI * f * t)
)
频率分辨率=采样速率/ FFT点数
采样率是由奈奎斯特绕圈,采样率确定(采样/秒)必须至少两倍要被分析的最高频率,例如48kHz的分析高达24kHz时。 (有关“现实生活”的数据,这是很好的有一个位的缓冲液)。
所以,你可能需要增加你的FFT的大小。
如果您是从一个函数生成,对与样本的工作,可以生成点的LOT和运行BIG FFT所以频率区间是用于高精度非常小。但它不会解决的基本问题。
首先,校正你的问题:(30 / 2PI)不是频率。你的信号的频率为1/30 *任何采样率你已经使用。 其次,你可以告诉我什么是sampleData在向量的长度?当FFT返回值的向量,第i个值将对应于f_i = I / N,其中N是向量的长度和I \在[0,N-1] 你想,我/ N正好等于1/30对于一些整数i。换句话说,N应该等于30 * 1,即,N应为30的现在,是你所用载体,为30的倍数?的长度如果不尝试使,和多应解决的问题。
尝试窗函数?