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?

15 views (last 30 days)

Edited: Paul on 17 Apr 2022
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.

Paul on 16 Apr 2022
Edited: Paul on 17 Apr 2022
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!

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!