Is there a function like "y = filter(b,a,x,zi) uses initial conditions zi for the filter delays" in fftfilt? If not, what's the most efficient way to implement this?

Is there anything like "y = filter(b,a,x,zi) uses initial conditions zi for the filter delays" in fftfilt? If not, what's the most efficient way to implement this by using fftfilt.

Thanks a lot! Every answer helps.

Answer 1

[Paul]

Filtering is a linear operation so the iniital condition response and the input response can be computed separately and then added together. filter() uses a Direct-Form II Tranposed realization. The initial condition response of the filter can be determined from the system and output matrices of that realization.

For example:

Generate an IIR filtr

rng(100)
b = rand(1,5);
a = [1 rand(1,4)]

a = 1×5

1.0000 0.1216 0.6707 0.8259 0.1367

The A-matrix and C-matrix of the DFIIT realization is

A = diag(ones(1,3),1);
A(:,1) = -a(2:end).';
C = zeros(1,size(A,1)); C(1) = 1;

Verify the characteristic polynomial

a

a = 1×5

1.0000 0.1216 0.6707 0.8259 0.1367

poly(A)

ans = 1×5

1.0000 0.1216 0.6707 0.8259 0.1367

Initial conditions for the taps

xi = 1:4;

Generate the 10 samples of the IC response using filter()

y1 = filter(b,a,zeros(1,10),xi);

Generate the IC response from the realization (could probalby be done more efficiently)

for n = 0:9
    y2(n+1) = C*(A^n)*xi(:);
end

Compare

[y1;y2]

ans = 2×10

1.0000 1.8784 2.1009 1.6588 -3.2988 -2.7034 0.8842 4.2034 1.5795 -3.3721 1.0000 1.8784 2.1009 1.6588 -3.2988 -2.7034 0.8842 4.2034 1.5795 -3.3721

However, the Question asks about using fftfilt(), which implies the underlying filter is a FIR filter. In this case, the initial condition response DFIIT filter is simply xi followed by zeros, which can be seen by inspection of the DFIIT realization.

y3 = filter(b,1,zeros(1,10),xi)

y3 = 1×10

1 2 3 4 0 0 0 0 0 0

So if we want to use fftfiltt() instead of filter() with initial conditions for a DFIIT realization:

u = rand(1,10); % random intput
y4 = filter(b,1,u,xi);
y5 = fftfilt(b,u) + [xi zeros(1,numel(u)-numel(xi))];
[y4;y5]

ans = 2×10

1.3125 2.6444 3.6059 5.0232 0.9550 0.4092 0.7965 0.8992 0.9209 1.6637 1.3125 2.6444 3.6059 5.0232 0.9550 0.4092 0.7965 0.8992 0.9209 1.6637

y4 - y5

ans = 1×10

1.0e-15 * -0.2220 -0.4441 -0.4441 0 -0.1110 -0.1110 -0.2220 -0.2220 -0.1110 -0.4441

[祥宇 崔]

Thanks! This is the answer I want!