How to filter data with given numerator and denominator coefficients of an IIR Filter, in MATLAB Implementation

Question

Coefficients are:

b = [9.02201251375918   0.485095309463300   -0.571759308660462  3.63774928924158    4.42592447788520    1.30093840271037    -0.898592924453941  -9.08965268561448   -12.4939309921467   -3.12442307305275   -0.384648004433882  5.71967440411881    0.143575175018517   -4.16015410667989   -1.39561043024989   3.68556878840504    1.37118920559575    9.05040751155859    5.14436163104499    1.50282327260915];

a = [1  -1.27434851603317   1.17767387410182    -0.0684817529891594 -0.731946865017767  -0.982532182263342  2.12595096907370    3.52705106824119    -4.03469353256380   0.932348168542137   1.53860232891322    2.29668138303561    -2.05888003206196   -2.16687979136895   -0.835820231426690  -1.41930565104868   -1.90147577647205   -2.92891900118810   2.84270807925957    -3.92597765206848];

According to @Hilmar @Matt L. and another kind comment, after I do abs(roots(a)), and find there're actually many values larger than 1, mine filter turned out to be an unstable one. Thank you again for your attention, I'm so appreciated your guidance.

1. Problem description

I'm fairly weak in DSP or filtering thing, but recently I'm trying to filter data with an IIR filter, using Y=filter(b,a,X) in MATLAB. I've got the numerator and denominator coefficients, and with freqz(b,a) I'm sure this is the filter I need. However, after I do filter(b,a,sig), the output seemed to be beyond my expectation. I've worked on it a few days with no progress.

Anyone may help me out here? What's wrong with my code, or what's exactly the right way to do with numerator and denominator coefficients to filter data?

2. About my IIR filter

I design my filter with numerator and denominator length both equal to 20, so the order is 38. With freqz(b,a) or fvtool(b,a), we can see the magnitude response as follows: The filter shall amplify the input signal with no more than 30 dB.

3. Main code

Here is my code, for brevity, I omit some coefficients:

clear;close all;clc;
%%% load data
[wav,Fs]= audioread('G:\Program\MATLAB\wavEQ\Experiment\raw.wav'); 
sig = wav(:,1);                           % data to be filtered
%%% get numerator and denominator coefficients
allCoef = [9.022012513759181  0.48509530946330026  ...];
len_b = ceil(length(allCoef)/2);
b = allCoef(1:len_b);
a(1) = 1;
a(2:length(allCoef)-len_b+1)= allCoef(len_b+1:length(allCoef));
%%% filtering
filt = filter(b, a, sig);                % IIR filtering
%%% sos-type filtering
[SOS,~] = tf2sos(b,a);
filt_sos = sosfilt(SOS, sig);
%%% plot unfiltered and filtered signal
T=(length(sig)-1)/Fs;
t=0:1/Fs:T;
figure(1)
subplot(321)
plot(t,sig);title('orgin signal');
subplot(323)
plot(t,filt,'r');title('filt');
subplot(325)
plot(t,filt_sos,'g');title('filt_sos');
subplot(322)
plot(t,db(sig));title('db(sig)');
subplot(324)
plot(t,db(filt),'r');title('db(filt)');
axis([0,4,-100,1000]);
subplot(326)
plot(t,db(filt_sos),'g');title('db(filt_sos)');
axis([0,4,-100,1000]);
% wav_filtering = [sig1_filtering,wav(:,2)];
% audiowrite('raw_FIR_filtering.wav', wav_filtering, Fs);

4. Filtering result

The filtering result was plotted as follow: From "origin signal" figure, we can see that the magnitude of unfiltered signal keeps at a level of less than 0.1, and it gets an mean value of 1.9337e-07, while the magnitude of the filtered signal seemed to be amplified to 'Nan' level by simply having a glance of mean(filt) or mean(filt_sos) and got 'NaN' as result. Even after plotting the filtered signal again but in dB scale and cut off y-axis, it doesn't work out. The filtered signal's amplitude is so much much higher than the origin signal's, which is quite not consistent with the filter magnitude response evidently.

By the way, the frequency range of the origin signal is between 0 and 10000 Hz, and the sampling rate of it is 51200 Hz, which corresponds to the sampling frequency of the filter.

5. My consideration

Firstly, I thought there must be something wrong with my usage of filter(b,a,X), then I help filter and learn that "The filter is a 'Direct Form II Transposed' implementation of the standard difference equation:". I tried to convert coefficients to SOS type with [sos,~]=tf2sos(b,a), but the filtering result was still far from expectation (using sosfilt(sos)).

I couldn't help believing that I may have trouble with understanding the relationship between transfer function and coefficients of the filters. But what happened with filter(b,a,X) since I could get ideal frequency response using the same "b" & "a" with freqz(b,a). Are "b" and "a" in these two functions represent different things? Any suggestion would be great appreciated.

Answer 1

I'm pretty sure that your filter is unstable. As mentioned in a comment, it would help if you provided the coefficients in an easily readable form to your question. You can try abs(roots(a)) to compute the magnitude of the filter's poles. If any of these values is larger than $1$ then the filter is unstable, and that would be the explanation why your output signal grows without bounds.

Note that even an unstable filter can have a finite magnitude response as long as there are no poles on the unit circle of the complex plane. So the plot of the magnitude response doesn't say much about stability.

It would be interesting to know how you came up with those filter coefficients. By the way, the filter order is

max( length(b), length(a) ) - 1

which is $19$ in your case.

EDIT: Now that I have the filter coefficients, I see that the filter is indeed unstable. The maximum pole radius is $1.4697$. If the phase response doesn't matter, you can replace all poles outside the unit circle by reflecting them inside the circle. Of course there mustn't be any poles on the circle. You can use the function polystab for generating stable denominator coefficients from the ones you have. Then you just have to scale the filter accordingly.

Answer 2

Most likely your filter is unstable, i.e some poles on or outside the unit circle. You can test this by executing

p = roots(a); 
display(abs(p) >= 1);

If you see any "1" in the result, you have bad poles.

This happens primarily through lack of numerical precision. To find the poles from the transfer function you have to calculate the roots of the polynomial which is a numerically tricky problem especially for higher order.

You need to get the original poles and zeros of your filter, once you have converted it to transfer function, there is no easy way back.

$freqz()$ doesn't care since it works directly in the Z domain and ignores any convergence issues.