IIR filter is unstable despite being modeled on real-world causal measurements.

1 view (last 30 days)
Hi Everyone,
I have a physical system consisting of a microphone, which feeds to a DSP, which drives a speaker. The microphone and speaker system have an amplitude response, as well as a phase response, which is clearly lagging. The aim is to produce an IIR filter in the DSP that can create a net loop gain of 1 and a phase shift of (180' + 360k'). The aim is that whatever the microphone reads, will produce a signal from the speaker that is out of phase with the input acoustic signal. I have measured the transfer function for the microphone speaker combination and as a start tried to make an IIR filter that replicates this. The phase at 20Hz lags by 90' and as 2kHz is approached it lags by 310'. The amplitude peaks at about 300Hz, gently rolling off either side of that. So, as this is a causal physical system, I assume that I can at least replicate this response in an IIR filter, before attempting my original goal.
clc
Frequency = [0; 20; 50; 100; 150; 200; 250; 300; 350; 400; 450; 500; 550; 600; 650; 700; 750; 800; 850; 900; 950; 1000; 1050; 1100; 1150; 1200; 1250; 1300; 1350; 1400; 1450; 1500; 1550; 1600; 1650; 1700; 1750; 1800; 1850; 1900; 1950; 2000;]; %Frequency vector in Hz
Mag_D = 0.001*[0; 42.105; 63.158; 78.947; 78.947; 118.421; 118.421; 126.316; 113.158; 107.895; 84.211; 71.053; 48.421; 50.526; 41.053; 40.000; 37.895; 33.684; 29.474; 28.421; 26.316; 24.211; 23.158; 23.158; 22.105; 21.053; 21.053; 20.000; 18.947; 15.789; 16.316; 15.789; 16.842; 15.789; 16.842; 16.842; 16.842; 17.368; 16.316; 15.263; 14.211; 13.158;]; %Magnitude vector in mV
Phi_D = [-90; -94; -140; -173; -135; -173; -194; -216; -232; -256; -262; -270; -269; -272; -290; -282; -281; -288; -269; -282; -287; -284; -295; -285; -286; -294; -297; -285; -301; -292; -287; -292; -285; -288; -285; -294; -290; -311; -306; -301; -298; -310;]; %Phase vector in degrees
Data_Array = horzcat(Frequency, Mag_D, Phi_D); %Put them all into an array for ease of use.
for n = 1:Rows %Make a complex representation off the magnitude and phase data.
Data_Array(n, 4) = Data_Array(n, 2).*exp(j*(Data_Array(n, 3)*2*pi/360)); %Convert to radians remember.
end
D = fdesign.arbmagnphase('Nb,Na,F,H',200,200,Frequency/2000,Data_Array(:, 4)); %Set the filter specification.
disp('---------------------------------------------------------')
designmethods(D) %Display the methods that may be used to design this filter.
disp('---------------------------------------------------------')
Hd = design(D,'iirls'); %Create the IIR filter
HFVT_Mag = fvtool(Hd, 'magnitude'); %We will visualise the filter now.
set(HFVT_Mag, 'MagnitudeDisplay', 'Magnitude (dB)', 'FrequencyScale', 'Log', 'NormalizedFrequency', 'Off');
HFVT_Phi = fvtool(Hd, 'phase');
set(HFVT_Phi, 'PhaseUnits', 'Degrees', 'FrequencyScale', 'Log', 'NormalizedFrequency', 'Off');
fcfwrite(Hd,'firfilter.txt'); %Write filter to a file.
Now, the resulting filter amplitude matches the specification, and so does the phase, except for the wrapping issue:
The filter coefficients generate just fine, except the report in Matlab says this:
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure : Direct-Form II
% Numerator Length : 201
% Denominator Length : 201
% Stable : No
% Linear Phase : No
What I am confused about is why the system is not stable? As far as I can see, the phase is lagging, so there is no non-causal behaviour. What could be wrong?
Any help would be appreciated!
Kind regards Charles
  1 Comment
Mark Knecht
Mark Knecht on 8 Dec 2014
Charles - this thread is quite old but the topic interests me. As a new MatLab user I wanted to run the code but 'Rows' is undefined. Do you remember what you wanted to use there? length(frequency)? Or something else?

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!