# IIR filter that is causal and stable and has more phase lag at lower frequencies

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 or gain, as well as a phase response, which is clearly lagging as the frequency increases.

The aim is to produce an IIR filter in the hardware of a DSP that can create a net loop gain of 1 and a phase shift of (180' + 360k'). I am limited to IIR filters. 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, even if it is a cycle or 2 delayed.

I have measured the transfer function for the microphone speaker combination and calculated the arbitrary magnitude and phase response required to meet my needs. Here is my Matlab code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  This script calculates the filter coefficients, from measured amplitude
%  and phase data.  The intention is to program an IIR filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc                                                     %Clear the creen from any previous runs

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;         2050;];     %Frequency values
Mag_SM      = 0.001*[20;    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;       10;];       %Magnitude gain of speaker microphone
Phi_SM      =       [-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;        -330;];      %Phase lag of speaker microphone
Phi_DSP     =       [3.5;   3.5;      1;          -1;         -1.75;      -2.25;      -3;         -3.75;      -4.5;       -5;         -5.5;       -6.25;      -7;         -7.5;       -8.5;       -9;         -9.5;       -10.25;     -10.75;     -11.5;      -12.25;     -13;        -13.5;      -14;        -15;        -15.5;      -16.25;     -17;        -17.5;      -18;        -18.75;     -19.25;     -20;        -20.75;     -21.25;     -22;        -22.75;     -23.25;     -24;        -24.75;     -25.5;      -26;        -26.5;];      %Phase lag of the DSP ADC and DAC hardware

Phi_Tot     =       Phi_SM + Phi_DSP;                   %Add the phase lag of the DSP to the speaker and microphone, as this is the total phase lag needed to cmopensate in the IIR

Mag_IIR     =       1./Mag_SM;                          %We want the response of the IIR to be the inverse
Phi_IIR     =       -360 - 180 - Phi_Tot;               %Make sure the IIR phase is always negative so as to be causal

[Rows Columns] = size(Frequency);                       %For counting and what not

Data_Array = horzcat(Frequency, Mag_IIR, Phi_IIR);      %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(1j*(Data_Array(n, 3)*2*pi/360));
end

D = fdesign.arbmagnphase('Nb,Na,F,H',31,31,Frequency/2050,Data_Array(:, 4));
Hd = design(D,'iirls');

hfvt = fvtool(Hd, 'Color','w');                         %For plotting the graphs
legend(hfvt,'Equiripple', 'FIR Least-Squares','Frequency Sampling', ...
'Location', 'NorthEast')
hfvt(2) = fvtool(Hd,'Analysis','phase','Color','white');
legend(hfvt(2),'Equiripple', 'FIR Least-Squares','Frequency Sampling')

pidx = find(Frequency>=0);
unwrap(angle(Data_Array(:,pidx:end)))],'k--')

fcfwrite(Hd,'firfilter.txt');                           %Write filter coefficients to a file


The magnitude plot is perfectly acceptable:

And the phase plot is correct to a certain extent:

It's not entirely correct because instead of spanning from 270' to 540', which is phase LEAD and non-causal, it should be spanning from -450' to -180', which is phase LAG and causal. As a result, because Matlab has created a filter that has a positive phase, it is non-causal and unstable:

% Coefficient Format: Decimal

% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure    : Direct-Form II
% Numerator Length    : 32
% Denominator Length  : 32
% Stable              : No
% Linear Phase        : No


I assume if Matlab were to put in the appropriate delays in the IIR, it COULD add another 720' onto the phase lag and produce a stable filter that is realizable real-time? Any ideas how to tell it to include the -720'?

Any help would be appreciated!

Kind regards Charles