首頁猿問為什么使用 FFT 在信號(hào)中舍入頻率值？

為什么使用 FFT 在信號(hào)中舍入頻率值？

Python

浮云間 2021-11-09 20:42:43

所以，我試圖弄清楚如何在實(shí)踐中使用 DFT 來檢測(cè)信號(hào)中的普遍頻率。我一直試圖圍繞傅立葉變換是什么以及 DFT 算法如何工作，但顯然我還有很長(zhǎng)的路要走。我編寫了一些代碼來生成信號(hào)（因?yàn)槟康氖翘幚硪魳?，所以我生成了一個(gè)主要的 C 和弦，因此是奇怪的頻率值），然后嘗試回到頻率數(shù)。這是我的代碼sr = 44100 # sample ratex = np.linspace(0, 1, sr) # one second of signaltpi = 2 * np.pidata = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)freqs = np.fft.fftfreq(sr)fft = np.fft.fft(data)idx = np.argsort(np.abs(fft))fft = fft[idx]freqs = freqs[idx]print(freqs[-6:] * sr)這給了我[-262. 262. -330. 330. -392. 392.]不同于我編碼的頻率（261.63、329.63 和 392.0）。我做錯(cuò)了什么，我該如何解決？

查看完整描述

3 回答

森欄

TA貢獻(xiàn)1810條經(jīng)驗(yàn) 獲得超5個(gè)贊

實(shí)際上，如果幀持續(xù)T數(shù)秒，則 DFT 的頻率是k/THz，其中 k 是整數(shù)。因此，過采樣不會(huì)提高估計(jì)頻率的準(zhǔn)確性，只要這些頻率被識(shí)別為 DFT 幅度的最大值。相反，考慮到持續(xù) 100 秒的較長(zhǎng)幀會(huì)導(dǎo)致 DFT 頻率之間的間距為 0.01Hz，這可能足以產(chǎn)生預(yù)期的頻率。通過將峰值頻率估計(jì)為其相對(duì)于功率密度的平均頻率，可能會(huì)更好。

圖 1：即使在應(yīng)用了 Tuckey 窗口之后，加窗信號(hào)的 DFT 也不是 Dirac 的總和：在峰值底部仍然存在一些頻譜泄漏。在估計(jì)頻率時(shí)必須考慮到該功率。

另一個(gè)問題是幀的長(zhǎng)度不是信號(hào)周期的倍數(shù)，它可能不是周期性的。盡管如此，DFT 的計(jì)算方式好像信號(hào)是周期性的，但在幀的邊緣是不連續(xù)的。它會(huì)引起被描述為頻譜泄漏的雜散頻率。加窗是處理此類問題和緩解與人為不連續(xù)相關(guān)問題的參考方法。實(shí)際上，窗口的值在靠近幀的邊緣處不斷減小到零。有一個(gè)窗口函數(shù)列表，在scipy.signal 中有很多窗口函數(shù)可用。一個(gè)窗口被應(yīng)用為：

tuckey_window=signal.tukey(len(data),0.5,True)

data=data*tuckey_window

在這一點(diǎn)上，出現(xiàn)最大幅度的頻率仍然是 262、330 和 392。應(yīng)用窗口只會(huì)使峰值更加明顯：加窗信號(hào)的 DFT 具有三個(gè)不同的峰值，每個(gè)都具有中心瓣和旁瓣，具體取決于窗口的 DFT。這些窗口的波瓣是對(duì)稱的：中心頻率因此可以計(jì)算為峰值的平均頻率，與功率密度有關(guān)。

import numpy as np

from scipy import signal

import scipy

sr = 44100 # sample rate

x = np.linspace(0, 1, sr) # one second of signal

tpi = 2 * np.pi

data = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)

#a window...

tuckey_window=signal.tukey(len(data),0.5,True)

data=data*tuckey_window

data -= np.mean(data)

fft = np.fft.rfft(data, norm="ortho")

def abs2(x):

return x.real**2 + x.imag**2

fftmag=abs2(fft)[:1000]

peaks, _= signal.find_peaks(fftmag, height=np.max(fftmag)*0.1)

print "potential frequencies ", peaks

#compute the mean frequency of the peak with respect to power density

powerpeak=np.zeros(len(peaks))

powerpeaktimefrequency=np.zeros(len(peaks))

for i in range(1000):

dist=1000

jnear=0

for j in range(len(peaks)):

if dist>np.abs(i-peaks[j]):

dist=np.abs(i-peaks[j])

jnear=j

powerpeak[jnear]+=fftmag[i]

powerpeaktimefrequency[jnear]+=fftmag[i]*i

powerpeaktimefrequency=np.divide(powerpeaktimefrequency,powerpeak)

print 'corrected frequencies', powerpeaktimefrequency

由此產(chǎn)生的估計(jì)頻率為 261.6359 Hz、329.637Hz 和 392.0088 Hz：它比 262、330 和 392Hz 好得多，并且滿足這種純無噪聲輸入信號(hào)所需的 0.01Hz 精度。

反對(duì) 回復(fù) 2021-11-09

繁星coding

TA貢獻(xiàn)1797條經(jīng)驗(yàn) 獲得超4個(gè)贊

DFT 結(jié)果箱在頻率上由 Fs/N 分隔，其中 N 是 FFT 的長(zhǎng)度。因此，DFT 窗口的持續(xù)時(shí)間限制了 DFT 結(jié)果箱頻率中心間距方面的分辨率。

但是，對(duì)于低噪聲（高 S/N）中分離良好的頻率峰值，您可以通過在 DFT 結(jié)果箱之間插入 DFT 結(jié)果來估計(jì)頻率峰值位置，以達(dá)到更高的分辨率，而不是增加數(shù)據(jù)的持續(xù)時(shí)間。您可以嘗試對(duì)粗略頻率峰值位置估計(jì)進(jìn)行拋物線插值，但加窗 Sinc 插值（本質(zhì)上是香農(nóng)-惠特克重建）將提供更好的頻率估計(jì)精度和分辨率（假設(shè)感興趣的頻率峰值周圍有足夠低的本底噪聲，例如，在您的人工波形情況下沒有附近的正弦曲線）。

反對(duì) 回復(fù) 2021-11-09