DSPRelated.com
Forums

MATLAB simulation of bilinear transformation

Started by simonzz 6 years ago2 replieslatest reply 6 years ago2619 views

Hello,

I am trying to simulate an IIR Butterworth filter using low-pass to band-stop and then using bilinear transformation and I am a little bit confused about the results I obtain for the IIR response.

The code I am using is the following:

% ex4a Bilinear IIR
Fs = 44100;
Fc = 4500;
Bw = 2000;

w = (6000 - 3000)/(5500 - 3500);
aMax = 10^6;
aMin = 10^.05;
% compute Butterworth analog filter order.
N = ceil(log((aMax - 1)/(aMin - 1))/(2*log(w)))

% get the zero poles of analog
% Butterworth low pass filter
[z,p,k] = buttap(N);
% get TF
[blp, alp] = zp2tf(z,p,k);
% get Band-stop analog TF
[bs, as] = lp2bs(blp, alp, Fc, Bw);
% plot bode diagram
f = figure();
bode(tf(bs, as));
title('Bode Diagram Butterworth notch');
waitfor(f);

% band-stop zero poles 
[zs, ps, ks] = tf2zp(bs, as);
% compute bilinear transformation
[zd,pd,kd] = bilinear(zs, ps, ks, Fs);
[sos] = zp2sos(zd,pd,kd);
% plot IIR frequency response
freqz(sos, 8192, Fs);
title('Frequency response of IIR Butterworth notch');

The problem I am having with the simulation is that when I look at the IIR frequency response with freqz, the cutoff frequency is centered at approx. 750 Hz.

The freqz plot looks good if I call the bilinear function this way:

[zd,pd,kd] = bilinear(zs, ps, ks, Fs/(2*pi));

However, the documentation for bilinear function specify that Fs is in Hz..

Could someone help to clarify this mismatch ?

Thank you.

Regards.

[ - ]
Reply by neiroberMarch 26, 2018

Simon,

You have a lot of steps, so it's hard to say quickly where the problem is.  You can design a band-reject filter in one step using:

[b,a]= butter(N,[fL fU]*2/fs,'stop')

where fL and fU are the -3 dB frequencies. It uses the bilinear transform with pre-warping.

You may also want to take a look at this Matlab function I wrote that designs band-reject filters:

https://www.dsprelated.com/showarticle/1131.php

regards,

Neil


[ - ]
Reply by simonzzMarch 26, 2018

Hello @neirobert,

You are right. The butter function already uses the bilinear transformation.

https://es.mathworks.com/help/signal/ref/butter.html

I completely missed that info yesterday. My bad.

So all the steps I have done are quiet unuseful. The good thing is that I have another code tested which uses butterord and butter functions wich works pretty fine.

Thank you for your post and for the article, I'll take a read to it. :)

EDIT: I understand that in my code I should apply pre-warping.

Regards,

Simon