Вычислительный фильтр (b, a, x, zi) с использованием БПФ - PullRequest
Новая
       41

Вычислительный фильтр (b, a, x, zi) с использованием БПФ

3 голосов
/ 12 декабря 2010

Я хотел бы попытаться вычислить y=filter(b,a,x,zi) и dy[i]/dx[j], используя FFT, а не во временной области для возможного ускорения в реализации GPU.

Я не уверен, что это возможно, особенно когда zi не равен нулю.Я посмотрел, как реализованы scipy.signal.lfilter в scipy и filter в октаве.Они оба делаются непосредственно во временной области, при этом scipy использует прямую форму 2 и октавную прямую форму 1 (из просмотра кода в DLD-FUNCTIONS/filter.cc).Я нигде не видел реализации FFT, аналогичной fftfilt для FIR-фильтров в MATLAB (то есть 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)

Ответы [ 3 ]

2 голосов
/ 12 декабря 2010

Вы можете уменьшить ошибку в существующей реализации, заменив P = T + MN - 1 на P = T + 2*MN - 1.Это чисто интуитивно понятно, но мне кажется, что для разделения bf и af потребуются 2*MN термины из-за переноса.

CS Burrus довольно кратко описывает, как рассматривать фильтрациюFIR или IIR, в блочно-ориентированном виде, здесь .Я не читал это подробно, но я думаю, что это дает вам уравнения, необходимые для реализации БИХ-фильтрации по свертке, включая промежуточные состояния.

1 голос
/ 12 декабря 2010

Попробуйте scipy.signal.lfiltic(b, a, y, x=None), чтобы получить начальные условия.

Текст документа для lfiltic: 

Given a linear filter (b,a) and initial conditions on the output y
and the input x, return the inital conditions on the state vector zi
which is used by lfilter to generate the output given the input.

If M=len(b)-1 and N=len(a)-1.  Then, the initial conditions are given
in the vectors x and y as

x = {x[-1],x[-2],...,x[-M]}
y = {y[-1],y[-2],...,y[-N]}

If x is not given, its inital conditions are assumed zero.
If either vector is too short, then zeros are added
  to achieve the proper length.

The output vector zi contains

zi = {z_0[-1], z_1[-1], ..., z_K-1[-1]}  where K=max(M,N).
1 голос
/ 12 декабря 2010

Я забыл, что мало что знал о БПФ, но вы можете взглянуть на sedit.py и interval.py на http://jc.unternet.net/src/ и посмотреть, поможет ли что-нибудь там.

...