1
$\begingroup$
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:magnitude response 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:enter image description here 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.

$\endgroup$
2
  • $\begingroup$ I think you should add your coefficients to the question numerically and not just as a screenshot. Additionally I would add the spectra of the unfiltered and filtered data. $\endgroup$ Sep 4, 2019 at 13:09
  • $\begingroup$ Thank you very much for your attention. It turned out that my filter is an unstable one, it did do with my coefficients. Plus, your suggestion teaches me how to describe my question more complete, I would remember that. $\endgroup$
    – Jack Wolf
    Sep 4, 2019 at 14:10

2 Answers 2

0
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Sir, you are right. Mine filter is unstable, I used to think it doesn't matter to do with "unstable" as long as the magnitude response is okay. I kinda know what I need to do next. As for the filter coefficients, I use genetic algorithm to search optimal coefficients "b" and "a", under the condition that freq(b,a) approximates the ideal magnitude response. Thank you so much for your attention and correction. $\endgroup$
    – Jack Wolf
    Sep 4, 2019 at 13:37
  • $\begingroup$ SUPER helpful!! After polystab,I find the magnitude response is almost the same as the unstable one, but with a all magnitude increase of 20 dB. So would you please give me a hint about scaling the filter? I guessed you might talk about polyscale(a,alpha), and I DID try to utilize it. But the shape of the magnitude response tends to change a little much, and I find it seems to be its inherent feature——"By reducing the radius of the roots in an autoregressive polynomial, the bandwidth of the spectral peaks in the frequency response is expanded (flattened)". $\endgroup$
    – Jack Wolf
    Sep 5, 2019 at 9:25
  • $\begingroup$ @JackWolf: If a2 are the result of polystab, and a are your original filter coefficients, then you can scale the b coefficients by abs(sum(a2)/sum(a)). This should give you the same magnitude as before, just with a stable filter. $\endgroup$
    – Matt L.
    Sep 5, 2019 at 10:20
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Your point would help me a lot in my next work, thanks a lot sir. $\endgroup$
    – Jack Wolf
    Sep 4, 2019 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.