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
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
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.