2
$\begingroup$

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')

ax=get(hfvt(2),'CurrentAxes'); set(ax,'NextPlot','Add');
pidx = find(Frequency>=0);
plot(ax,Frequency,[fliplr(unwrap(angle(Data_Array(:,pidx-1:-1:1)))) ... % Mask
    unwrap(angle(Data_Array(:,pidx:end)))],'k--')

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

The magnitude plot is perfectly acceptable: enter image description here

And the phase plot is correct to a certain extent: enter image description here

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

$\endgroup$
1
$\begingroup$

A pole inside the unit circle will cause the phase angle (representing phase lag), of that pole's contribution, to increase as the angle around the unit circle (representing frequency) increases. A pole outside does the opposite. So your choices for a generally decreasing phase angle around the unit circle are: (1) poles outside, an unstable filter, or (2) more zeros inside than poles + zeros outside, representing a FIR filter (instead of, or in conjunction with an IIR).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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