标签：filter octave matlab python signal-processing
我想尝试使用FFT而不是在时域中计算y = filter(b,a,x,zi)和dy [i] / dx [j]来实现GPU实现中的加速.
我不确定是否有可能,尤其是当zi为非零时.我研究了如何实现scipy中的scipy.signal.lfilter和八度中的filter.它们都在时域中直接完成,使用直接格式2和八度直接格式1 scipy(通过在DLD-FUNCTIONS / filter.cc中的代码查找).我在MATLAB的FIR滤波器中都没有看到类似于fftfilt的FFT实现(即a = [1.]).
我尝试做y = ifft(fft(b)/ fft(a)* fft(x)),但这在概念上似乎是错误的.另外,我不确定如何处理初始瞬态zi.任何引用,指向现有实现的指针,将不胜感激.
示例代码,
import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt
# create an IRR lowpass filter
N = 5
b, a = sg.butter(N, .4)
MN = max(len(a), len(b))
# create a random signal to be filtered
T = 100
P = T + MN - 1
x = np.random.randn(T)
zi = np.zeros(MN-1)
# time domain filter
ylf, zo = sg.lfilter(b, a, x, zi=zi)
# frequency domain filter
af = sg.fft(a, P)
bf = sg.fft(b, P)
xf = sg.fft(x, P)
yfft = np.real(sg.ifft(bf/af * xf))[:T]
# error
print np.linalg.norm(yfft - ylf)
# plot, note error is larger at beginning and with larger N
plt.figure(1)
plt.clf()
plt.plot(ylf)
plt.plot(yfft)
解决方法:
您可以通过将P = T MN-1替换为P = T 2 * MN-1来减少现有实现中的错误.这纯粹是直观的,但在我看来bf和af的除法将需要2 * MN项,由于需要环绕.
CS Burrus有一篇非常简洁的文章,介绍了如何以面向块的方式考虑滤波,无论是FIR还是IIR,为here.我没有详细阅读它,但是我认为它为您提供了实现IIR滤波所需的方程式卷积,包括中间状态.
标签：filter,octave,matlab,python,signal-processing 来源： https://codeday.me/bug/20191105/1996342.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。