Is there any other way to find Transfer function H(z) for an IIR Butterworth Filter?

close all;
clc;clear all;


%ButterWorth IIR LowPass Filter
Rp=input('Enter PassBand Ripple:');
As=input('Enter Stopband Ripple:');
wp=input('Enter Digital Passband Frequency:');
ws=input ('Enter Digital Stopband Frequency:');
T = input('Enter the Sampling Time Period in Seconds:');

%PreWarping these frequencies into digital domain
omegap = (2/T)*tan(wp/2);
omegas = (2/T)*tan(ws/2); 

%Calculating order
num=log10((10^(0.1*Rp)-1)/((10^(0.1*As)-1)));
den=2*log10(omegap/omegas);
N=ceil(num/den);

%Calculating Cut off Frequencies
x1=(10^(0.1*Rp)-1);
x2=(10^(0.1*As)-1);
pow=1/(2*N);
omegac1=omegap/(x1^(pow));
omegac2=omegas/(x2^pow);
omegac=(omegac1+omegac2)/2;

%Calculating Poles
limit=((2*N)-1);
p1=zeros(1,limit+1);
for k=0:limit
    pow1=(1j*pi)/(2*N);
    e=exp(pow1*((2*k)+N+1));
    p1(k+1)=omegac*e;
end
p=zeros(1,N);

%Selection of left Half plane Poles
for k=0:limit
    if(real(p1(k+1))<0)
        p(k+1)=p1(k+1);
    end
end

%Determining Analog Prototype Filter Transfer Function Ha(s)
syms s z
den1=1;
num1=omegac^N;
lt=length(p);
for g=1:lt
    den1=den1*(s-p(g));
end
Ha=num1/den1;


%Bilinear Transformation
H=subs(Ha,s,(2/T)*(z-1)/(z+1));

H gives the tranfer function , but the z terms are just as it is. The desired result is that it should be in terms of z^2, z^3 ... etc

Is there any other way to obtain transfer function H(z)? Also How can i find difference equation from H(z)?

[Paul]

Hi Aryan,

It will be easier to help if you provide example values of Rp, As, wp, ws, and T so the result in question can be recreated.

[Walter Roberson]

expand perhaps

Or numden and then coeffs the pieces ?

[Star Strider]

One missing item are the units of the inputs. I partially corrected for them here, although the magnitude (checked using the freqz function) is still not correct. I leave those to you.

% close all;
% clc;clear all;

% % % SPECIFICATIONS: Rp = 1 Db,  As = 50 dB,  wp = 10 Hz,  ws = 20 Hz,  Fs = 100 Hz

%ButterWorth IIR LowPass Filter
% Rp=input('Enter PassBand Ripple:');
% As=input('Enter Stopband Ripple:');
% wp=input('Enter Digital Passband Frequency:');
% ws=input ('Enter Digital Stopband Frequency:');
% T = input('Enter the Sampling Time Period in Seconds:');

Rp = db2mag(-1);                                                % Convert From dB To Magnitude
As = db2mag(-50);                                               % Convert From dB To Magnitude
T = 0.01;                                                       % Sampling Interval (1/Fs)
wp = 10*pi*T*2;                                                 % Converted From Hz To rad/s & Corrected For Nyquist Frequency
ws = 20*pi*T*2;                                                 % Converted From Hz To rad/s & Corrected For Nyquist Frequency

%PreWarping these frequencies into digital domain
omegap = (2/T)*tan(wp/2);
omegas = (2/T)*tan(ws/2); 

%Calculating order
num=log10((10^(0.1*Rp)-1)/((10^(0.1*As)-1)))

num = 2.4952

num = abs(num)

num = 2.4952

den=2*log10(omegap/omegas)

den = -0.6990

den = abs(den)

den = 0.6990

N=ceil(num/den)

N = 4

%Calculating Cut off Frequencies
x1=(10^(0.1*Rp)-1);
x2=(10^(0.1*As)-1);
pow=1/(2*N);
omegac1=omegap/(x1^(pow));
omegac2=omegas/(x2^pow);
omegac=(omegac1+omegac2)/2;

%Calculating Poles
limit=((2*N)-1);
p1=zeros(1,limit+1);
for k=0:limit
    pow1=(1j*pi)/(2*N);
    e=exp(pow1*((2*k)+N+1));
    p1(k+1)=omegac*e;
end
p=zeros(1,N);

%Selection of left Half plane Poles
for k=0:limit
    if(real(p1(k+1))<0)
        p(k+1)=p1(k+1);
    end
end

%Determining Analog Prototype Filter Transfer Function Ha(s)
syms s z
den1=1;
num1=omegac^N;
lt=length(p);
for g=1:lt
    den1=den1*(s-p(g));
end
Ha=num1/den1;


%Bilinear Transformation
H=subs(Ha,s,(2/T)*(z-1)/(z+1))

H = 

[Hn,Hd] = numden(H)

Hn = 

Hd = 

Hnn = sym2poly(Hn)

Hnn = 1×5

1.0e+64 * 0.6966 2.7865 4.1797 2.7865 0.6966

Hnd = sym2poly(Hd)

Hnd = 
   1.0e+64 *

   6.2490 - 0.0000i   1.3617 - 0.0000i   3.1310 - 0.0000i   0.2881 - 0.0000i   0.1161 - 0.0000i

figure
freqz(Hnn, Hnd, 2^16, 1/double(T))
set(subplot(2,1,1), 'YLim',[-10 0])                                 % Zoom Y-Axis Of Magnitude Plot (Optional)

.

Answer 1

[Sai]

I understand that you are trying to get your final transfer function in powers of z. With small changes in code without modifying your logic, you can get it as shown in below code snippet.

close all; 
clc;clear all; 
%ButterWorth IIR LowPass Filter 
Rp=input('Enter PassBand Ripple:'); 
As=input('Enter Stopband Ripple:'); 
wp=input('Enter Digital Passband Frequency:'); 
ws=input ('Enter Digital Stopband Frequency:'); 
T = input('Enter the Sampling Time Period in Seconds:'); 
%PreWarping these frequencies into digital domain 
omegap = (2/T)*tan(wp/2); 
omegas = (2/T)*tan(ws/2);  
%Calculating order 
num=log10((10^(0.1*Rp)-1)/((10^(0.1*As)-1))); 
den=2*log10(omegap/omegas); 
N=abs(ceil(num/den)); 
%Calculating Cut off Frequencies 
x1=(10^(0.1*Rp)-1); 
x2=(10^(0.1*As)-1); 
pow=1/(2*N); 
omegac1=omegap/(x1^(pow)); 
omegac2=omegas/(x2^pow); 
omegac=(omegac1+omegac2)/2; 
%Calculating Poles 
limit=((2*N)-1); 
p1=zeros(1,limit+1); 
for k=0:limit 
    pow1=(1j*pi)/(2*N); 
    e=exp(pow1*((2*k)+N+1)); 
    p1(k+1)=omegac*e; 
end 
p=zeros(1,N); 
%Selection of left Half plane Poles 
for k=0:limit 
    if(real(p1(k+1))<0) 
        p(k+1)=p1(k+1); 
    end 
end 
num1=omegac^N; 
den1 = p; 
%Bilinear Transformation 
[zd,pd] = bilinear(num1,den1,1,T)
H = (tf(zd,pd,T)) 

There are direct MATLAB functions like buttord (for determining filter order and cutoff frequency), butter(for obtaining filter coefficients), bilinear (for converting analog to digital) by which you can design Butterworth filters. Refer to the below documentation pages for more information

• https://in.mathworks.com/help/signal/ref/buttord.html
• https://in.mathworks.com/help/signal/ref/butter.html#buct21r
• https://in.mathworks.com/help/signal/ref/bilinear.html