Designing IIR Filters using Bilinear Transform

I have written the following code to design an IIR filter with a cut-off frequency of 800 Hz. However, I'm not getting the desired cut-off in the plot. Can someone tell me where I went wrong ?

clc;
clear all;
close all;

N = 2;
Wc = 800;
FS = 8000;
wc = 2*pi*Wc/FS;
T = 1;

[a, b] = butter(N, wc, 's');
[ad, bd] = bilinear(a, b, T);
freqz(ad,bd,512,8000); 
axis([0 4000 -40 1]);
title('Frequency Response of the Filter')

Answer 1

[Mathieu NOE]

so final and last version of the answer :

N = 2;
fc = [500 1000];
FS = 8000;
fcn = 2*fc/FS;
T = 1/FS;
freq = logspace(2,log10(FS/2.56));

% [ad, bd] = butter(N, fcn); % digital version
[a, b] = butter(N, 2*pi*fc, 's'); % analog version
figure(1); % allows log x 
[g,p]  = bode(a,b,2*pi*freq);grid
subplot(2,1,1),semilogx(freq,20*log10(g));
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid

[ad, bd] = bilinear(a, b, FS);

figure(2); freqz(ad,bd,freq,FS);  % allows only lin x 
title('Frequency Response of the Filter')

figure(3); % allows log x 
[g,p]  = dbode(ad,bd,1/FS,2*pi*freq);
subplot(2,1,1),semilogx(freq,20*log10(g));grid
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid
% axis([0 4000 -40 1]);

[ad, bd] = butter(N, fcn); % digital version
figure(4); % allows log x 
[g,p]  = dbode(ad,bd,1/FS,2*pi*freq);
subplot(2,1,1),semilogx(freq,20*log10(g));grid
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid

Answer 2

[Mathieu NOE]

hello

this works - simply wrong computation of normalized cut off frequency

FYI butter will by default generate a digital filter; no need to create an first analog version to discretize with bilinear

also I tested the freqz display (linear frequency axis) vs traditionnal log x bode plot;

clc;
clear all;
close all;

N = 2;
fc = 800;
FS = 8000;
fcn = 2*fc/FS;
T = 1/FS;

[ad, bd] = butter(N, fcn);

freq = logspace(2,log10(FS/2.56));

figure(1); freqz(ad,bd,freq,FS);  % allows only lin x 
title('Frequency Response of the Filter')

figure(2); % allows log x 
[g,p]  = dbode(ad,bd,1/FS,2*pi*freq);
subplot(2,1,1),semilogx(freq,20*log10(g));
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);
% axis([0 4000 -40 1]);

[Vincent Abraham]

I agree with your point but I want to create it using bilinear() and using analog approximation

[Mathieu NOE]

sure _ there is nothing that speaks against that

so here we are :

N = 2;
fc = 800;
FS = 8000;
fcn = 2*fc/FS;
T = 1/FS;
freq = logspace(2,log10(FS/2.56));

% [ad, bd] = butter(N, fcn); % digital version
[a, b] = butter(N, 2*pi*fc, 's'); % analog version
figure(1); % allows log x 
[g,p]  = bode(a,b,2*pi*freq);grid
subplot(2,1,1),semilogx(freq,20*log10(g));
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid

[ad, bd] = bilinear(a, b, FS);

figure(2); freqz(ad,bd,freq,FS);  % allows only lin x 
title('Frequency Response of the Filter')

figure(3); % allows log x 
[g,p]  = dbode(ad,bd,1/FS,2*pi*freq);
subplot(2,1,1),semilogx(freq,20*log10(g));grid
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid
% axis([0 4000 -40 1]);

[Vincent Abraham]

Okay yes, this works. But if i wanted to do the same for a BPF, I'm not getting the desired result.

For example,

clc;
clear all;
close all;

N = 2;
Wc = [1000 1400];
FS = 8000;
wc = 2*pi*Wc;

[a, b] = butter(N, wc, 's');
[ad, bd] = bilinear(a, b, FS);
freqz(ad, bd, 512, 8000); 
title('Frequency Response of the Filter')

[Mathieu NOE]

it works in my code :

N = 2;
fc = [1000 1400];
FS = 8000;
fcn = 2*fc/FS;
T = 1/FS;
freq = logspace(2,log10(FS/2.56));

% [ad, bd] = butter(N, fcn); % digital version
[a, b] = butter(N, 2*pi*fc, 's'); % analog version
figure(1); % allows log x 
[g,p]  = bode(a,b,2*pi*freq);grid
subplot(2,1,1),semilogx(freq,20*log10(g));
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid

[ad, bd] = bilinear(a, b, FS);

figure(2); freqz(ad,bd,freq,FS);  % allows only lin x 
title('Frequency Response of the Filter')

figure(3); % allows log x 
[g,p]  = dbode(ad,bd,1/FS,2*pi*freq);
subplot(2,1,1),semilogx(freq,20*log10(g));grid
title('Frequency Response of the Filter')
subplot(2,1,2),semilogx(freq,p);grid
% axis([0 4000 -40 1]);

[Vincent Abraham]

I'm not getting the cut-off frequency of 1000 and 1400Hz with your code

[Mathieu NOE]

yes they are correct

the cut off points are defined by where the modulus is 3 dB lower compared to max gain (here 0 dB)

so the results are in accordance with the theory

there could be a minor mismatch due to the bilinear approximation of the analog version, but if you use the native digital filter output of butter I can testify that it works

[Vincent Abraham]

Damn okay thanks a lot!

[Mathieu NOE]

you're welcome

I will put the last version of the code in the answer section, if you don't mind accep it ! tx