1
$\begingroup$

I have tried to make a function to implement Low pass filter using butter worth technique,and i have taken code from different places in article 8.3 of proakis book,"DSP Using Matlab 3rd Ed" I have made a single function combining all those code sections ,my function is running and give correct output [C,B,A] matrix as per book,but my plot has reverse direction(My plot appears to be of High pass) while i intended in my design low pass

My code is below:

function buttr_low_pass
%main function body
clc;clear all;close all
Wp = 0.2*pi; Ws = 0.3*pi; Rp = 7; As = 16;
% Analog filter design:
[b,a] = afd_butt(Wp,Ws,Rp,As)
freqz(b,a)
[C,B,A] = sdir2cas(b,a)

%definitions of called/used functions
function [C,B,A] = sdir2cas(b,a);

% DIRECT-form to CASCADE-form conversion in s-plane

% -------------------------------------------------

% [C,B,A] = sdir2cas(b,a)

%  C = gain coefficient

%  B = K by 3 matrix of real coefficients containing bk's

%  A = K by 3 matrix of real coefficients containing ak's

%  b = numerator polynomial coefficients of DIRECT form

%  a = denominator polynomial coefficients of DIRECT form

%

Na = length(a)-1; Nb = length(b)-1;



% compute gain coefficient C

b0 = b(1); b = b/b0;

a0 = a(1); a = a/a0;

 C = b0/a0;

%

% Denominator second-order sections:

p= cplxpair(roots(a)); K = floor(Na/2);

if K*2 == Na     % Computation when Na is even

   A = zeros(K,3);

   for n=1:2:Na

       Arow = p(n:1:n+1,:);

       Arow = poly(Arow);

       A(fix((n+1)/2),:) = real(Arow);

   end



elseif Na == 1   % Computation when Na = 1

       A = [0 real(poly(p))];



else             % Computation when Na is odd and > 1

   A = zeros(K+1,3);

   for n=1:2:2*K

       Arow = p(n:1:n+1,:);

       Arow = poly(Arow);

       A(fix((n+1)/2),:) = real(Arow);

       end

       A(K+1,:) = [0 real(poly(p(Na)))];

end



% Numerator second-order sections:

z = cplxpair(roots(b)); K = floor(Nb/2);

if Nb == 0           % Computation when Nb = 0

   B = [0 0 poly(z)];



elseif K*2 == Nb     % Computation when Nb is even

   B = zeros(K,3);

   for n=1:2:Nb

       Brow = z(n:1:n+1,:);

       Brow = poly(Brow);

       B(fix((n+1)/2),:) = real(Brow);

   end



elseif Nb == 1       % Computation when Nb = 1

       B = [0 real(poly(z))];



else                 % Computation when Nb is odd and > 1

   B = zeros(K+1,3);

   for n=1:2:2*K

       Brow = z(n:1:n+1,:);

       Brow = poly(Brow);

       B(fix((n+1)/2),:) = real(Brow);

   end

   B(K+1,:) = [0 real(poly(z(Nb)))];

end





function [b,a] = afd_butt(Wp,Ws,Rp,As);

% Analog Lowpass Filter Design: Butterworth

% -----------------------------------------

% [b,a] = afd_butt(Wp,Ws,Rp,As);

%  b = Numerator coefficients of Ha(s)

%  a = Denominator coefficients of Ha(s)

% Wp = Passband edge frequency in rad/sec; Wp > 0

% Ws = Stopband edge frequency in rad/sec; Ws > Wp > 0

% Rp = Passband ripple in +dB; (Rp > 0)

% As = Stopband attenuation in +dB; (As > 0)

%

if Wp <= 0

        error('Passband edge must be larger than 0')

end

if Ws <= Wp

        error('Stopband edge must be larger than Passband edge')

end

if (Rp <= 0) | (As < 0)

        error('PB ripple and/or SB attenuation ust be larger than 0')

end



N = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));

fprintf('\n*** Butterworth Filter Order = %2.0f \n',N)

OmegaC = Wp/((10^(Rp/10)-1)^(1/(2*N)));

[b,a]=u_buttap(N,OmegaC);
function [b,a] = u_buttap(N,Omegac);

% Unnormalized Butterworth Analog Lowpass Filter Prototype

% --------------------------------------------------------

% [b,a] = u_buttap(N,Omegac);

%      b = numerator polynomial coefficients of Ha(s)

%      a = denominator polynomial coefficients of Ha(s)

%      N = Order of the Butterworth Filter

% Omegac = Cutoff frequency in radians/sec

%

[z,p,k] = buttap(N);

      p = p*Omegac;

      k = k*Omegac^N;

      B = real(poly(z));

      b0 = k;

      b = k*B;

      a = real(poly(p));

Please save it in your PC as "buttr_low_pass.m" and then run in your Matlab

I have also attached snap of my output plot,which appears to be high passenter image description here

$\endgroup$
  • $\begingroup$ could you please remove all the superfluous empty lines from your code to make it sensibly readable? $\endgroup$ – Marcus Müller Jul 15 at 13:55
  • $\begingroup$ Can you add the plots to your question, too, so that we don't have to run your code? $\endgroup$ – Marcus Müller Jul 15 at 13:58
3
$\begingroup$

You should use freqs instead of freqz.

| improve this answer | |
$\endgroup$
  • $\begingroup$ This Chp 8, is about IIR filter design,which are part of digital filters category, so should not we use "freqz" command because IIR filters are digital not analog $\endgroup$ – engr Jul 15 at 17:34
2
$\begingroup$

Your problem is here, as Yiftah also pointed out:

% Analog filter design:
[b,a] = afd_butt(Wp,Ws,Rp,As)
freqz(b,a)

The a and b are not vectors of IIR filter coefficients or z-domain polynomials, they are s-domain polynomials.

The code above uses the analog prototypes to calculate the s-domain polynomials.

The name of the section also hints at it:

8.3. Characteristics of Prototype Analog Filters

You have to wait until the next section to do the s-to-z domain transfers.

8.4. Analog-to-Digital Filter Transformations

| improve this answer | |
$\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.