我目前正在解决通过Z变换计算的频率响应函数与numpy的FFT算法之间的一些差异。它是一个由脉冲响应代表的简单回声:
h[n] = δ[n] + αδ[n-L]
其中α是回波的衰减系数,L是样本中的延迟量。相应的传递函数由下式给出:
H(f) = ( e^(j2πfΔL) + α ) / e^(j2πfΔL)
其中Δ
是抽样期。
当直接绘制上面的传递函数幅度和使用numpy的fft算法时,我似乎使用相同数量的频率区来获得不同的结果。特别是,FFT幅度似乎在整个频谱周围形成一个包络 - 我相信我应该得到一个简单的梳状滤波器,就像传递函数方法一样:imgur
任何人都可以澄清为什么会发生这种情况以及我是否有可能忽视任何事情?这是由于所使用的DFT算法中的错误吗?
欣赏你的时间,欢呼!
import matplotlib.pyplot as pyplt
import numpy as np
fs = 48000 # Sample rate
T = 1/fs # Sample period
L = 3000 # Delay
a = 0.5 # Attenuation factor
# h[n] = dirac[n] + a * dirac[n-L]
h = np.zeros(L)
h[0] = 1
h[L-1] = a
# Transfer function H
freqs = np.arange(0, fs, fs/(2*L))
e = np.exp(1j*freqs*2*np.pi*L*T)
H = (e + a)/(e)
# Transfer function H via FFT - same # of bins
H_FFT = np.fft.fft(h, 2*L)
pyplt.figure()
# Correct comb filter
pyplt.plot(np.abs(H))
# Runing FFT gives a form of strange envelope error
pyplt.plot(np.abs(H_FFT))
pyplt.legend()
你的代码几乎没问题。你需要改变什么:
这是一个修改后的代码,用于显示其工作原理,为清晰起见绘制在两个不同的图中:
import matplotlib.pyplot as plt
import numpy as np
fs = 48000 # Sample rate
T = 1/fs # Sample period
L = 3 # Delay
a = 0.5 # Attenuation factor
# h[n] = dirac[n] + a * dirac[n-L]
h = np.zeros(fs)
h[0] = 1
h[L] = a
# Transfer function H
freqs = np.arange(0, fs)
e = np.exp(1j*freqs*2*np.pi*L*T)
H = (e + a)/(e)
# Transfer function H via FFT - same # of bins
H_FFT = np.fft.fft(h)
# Correct comb filter
plt.figure()
plt.plot(np.abs(H))
# Runing FFT gives a form of strange envelope error
plt.figure()
plt.plot(np.abs(H_FFT))