MATLAB simulation of bilinear transformation

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.

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

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