Why does notch filter coefficients have 10(5 numerator and 5 denominator) for butterworth order 2

Hi ,

Implemented butterworth stopband(notch) filter, order is 2

Getting 10 coefficients,

Designed in matlab

Can you please explain why I am getting 10 coefficients, Is is correct or something wwrong with my interpretation.

Answer 1

[Mathieu NOE]

hello

your IIR filter size is much bigger than needed for an order 2

see example below (search for the notch filter )

% 
%          Cookbook formulae for audio EQ biquad filter coefficients
% ----------------------------------------------------------------------------
%            by Robert Bristow-Johnson  <[e-mail]>
% 
% 
% All filter transfer functions were derived from analog prototypes (that
% are shown below for each EQ filter type) and had been digitized using the
% Bilinear Transform.  BLT frequency warping has been taken into account for
% both significant frequency relocation (this is the normal "prewarping" that
% is necessary when using the BLT) and for bandwidth readjustment (since the
% bandwidth is compressed when mapped from analog to digital using the BLT).
% 
% First, given a biquad transfer function defined as:
% 
%             b0 + b1*z^-1 + b2*z^-2
%     H(z) = ------------------------                                  (Eq 1)
%             a0 + a1*z^-1 + a2*z^-2
% 
% This shows 6 coefficients instead of 5 so, depending on your architechture,
% you will likely normalize a0 to be 1 and perhaps also b0 to 1 (and collect
% that into an overall gain coefficient).  Then your transfer function would
% look like:
% 
%             (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
%     H(z) = ---------------------------------------                   (Eq 2)
%                1 + (a1/a0)*z^-1 + (a2/a0)*z^-2
% 
% or
% 
%                       1 + (b1/b0)*z^-1 + (b2/b0)*z^-2
%     H(z) = (b0/a0) * ---------------------------------               (Eq 3)
%                       1 + (a1/a0)*z^-1 + (a2/a0)*z^-2
% 
% 
% The most straight forward implementation would be the "Direct Form 1"
% (Eq 2):
% 
%     y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]
%                         - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]            (Eq 4)
% 
% This is probably both the best and the easiest method to implement in the
% 56K and other fixed-point or floating-point architechtures with a double
% wide accumulator.
% 
% 
% 
% Begin with these user defined parameters:
% 
%     Fs (the sampling frequency)
% 
%     f0 ("wherever it's happenin', man."  Center Frequency or
%         Corner Frequency, or shelf midpoint frequency, depending
%         on which filter type.  The "significant frequency".)
% 
%     dBgain (used only for peaking and shelving filters)
% 
%     Q (the EE kind of definition, except for peakingEQ in which A*Q is
%         the classic EE Q.  That adjustment in definition was made so that
%         a boost of N dB followed by a cut of N dB for identical Q and
%         f0/Fs results in a precisely flat unity gain filter or "wire".)
% 
%      _or_ BW, the bandwidth in octaves (between -3 dB frequencies for BPF
%         and notch or between midpoint (dBgain/2) gain frequencies for
%         peaking EQ)
% 
%      _or_ S, a "shelf slope" parameter (for shelving EQ only).  When S = 1,
%         the shelf slope is as steep as it can be and remain monotonically
%         increasing or decreasing gain with frequency.  The shelf slope, in
%         dB/octave, remains proportional to S for all other values for a
%         fixed f0/Fs and dBgain.
% 
% 
% 
% Then compute a few intermediate variables:
% 
%     A  = sqrt( 10^(dBgain/20) )
%        =       10^(dBgain/40)     (for peaking and shelving EQ filters only)
% 
%     w0 = 2*pi*f0/Fs
% 
%     cos(w0)
%     sin(w0)
% 
%     alpha = sin(w0)/(2*Q)                                       (case: Q)
%           = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) )           (case: BW)
%           = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 )         (case: S)
% 
%         FYI: The relationship between bandwidth and Q is
%              1/Q = 2*sinh(ln(2)/2*BW*w0/sin(w0))     (digital filter w BLT)
%         or   1/Q = 2*sinh(ln(2)/2*BW)             (analog filter prototype)
% 
%         The relationship between shelf slope and Q is
%              1/Q = sqrt((A + 1/A)*(1/S - 1) + 2)
% 
%     2*sqrt(A)*alpha  =  sin(w0) * sqrt( (A^2 + 1)*(1/S - 1) + 2*A )
%         is a handy intermediate variable for shelving EQ filters.
% 
% 
% Finally, compute the coefficients for whichever filter type you want:
%    (The analog prototypes, H(s), are shown for each filter
%         type for normalized frequency.)

%%%%%%%%%%%%%%%%%%%% inputs %%%%%%%%%%%%%%%
Fs = 1e4;
f0 = 100;

%%%%%%%%%%%%%%%%%%%% outputs %%%%%%%%%%%%%%%
w0 = 2*pi*f0/Fs;


%%%%%%%%%%%%%%% simulation %%%%%%%%%%%%%%%
% freq = logspace(1,(log10(Fs/2.5)),300);
freq = logspace(1,3,300);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LPF:        H(s) = 1 / (s^2 + s/Q + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Q = 10;
alpha = sin(w0)/(2*Q);

            b0 =  (1 - cos(w0))/2;
            b1 =   1 - cos(w0);
            b2 =  (1 - cos(w0))/2;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = mydbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(1);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title(' LPF: H(s) = 1 / (s^2 + s/Q + 1)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig            
            


% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HPF:        H(s) = s^2 / (s^2 + s/Q + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = 10;
alpha = sin(w0)/(2*Q);

            b0 =  (1 + cos(w0))/2;
            b1 = -(1 + cos(w0));
            b2 =  (1 + cos(w0))/2;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = mydbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(2);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title(' HPF: H(s) = s^2 / (s^2 + s/Q + 1)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig
            
   
             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BPF:        H(s) = (s/Q) / (s^2 + s/Q + 1)      (constant 0 dB peak gain)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = 10;
alpha = sin(w0)/(2*Q);

            b0 =   alpha;
            b1 =   0;
            b2 =  -alpha;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(3);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title(' BPF: H(s) = (s/Q) / (s^2 + s/Q + 1)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% notch :      H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = 10;
alpha = sin(w0)/(2*Q);


            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
% 
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(4);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APF:        H(s) = (s^2 - s/Q + 1) / (s^2 + s/Q + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = 10;
alpha = sin(w0)/(2*Q);

            b0 =   1 - alpha;
            b1 =  -2*cos(w0);
            b2 =   1 + alpha;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;

% 
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(5);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title(' APF: H(s) = (s^2 - s/Q + 1) / (s^2 + s/Q + 1)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig% % 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% peakingEQ:  H(s) = (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 4;
Q = 3;
alpha = sin(w0)/(2*Q);

            b0 =   1 + alpha*A;
            b1 =  -2*cos(w0);
            b2 =   1 - alpha*A;
            a0 =   1 + alpha/A;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha/A;

num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(6);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title(' peakingEQ:  H(s) = (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig% % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lowShelf: H(s) = A * (s^2 + (sqrt(A)/Q)*s + A)/(A*s^2 + (sqrt(A)/Q)*s + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
A = 4;
Q = 3;
alpha = sin(w0)/(2*Q);

            b0 =    A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha );
            b1 =  2*A*( (A-1) - (A+1)*cos(w0)                   );
            b2 =    A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha );
            a0 =        (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha;
            a1 =   -2*( (A-1) + (A+1)*cos(w0)                   );
            a2 =        (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha;

num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(7);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title('lowShelf: H(s) = A * (s^2 + (sqrt(A)/Q)*s + A)/(A*s^2 + (sqrt(A)/Q)*s + 1)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
fixfig% % % 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% 
A = 4;
Q = 3;
alpha = sin(w0)/(2*Q);

            b0 =    A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha );
            b1 = -2*A*( (A-1) + (A+1)*cos(w0)                   );
            b2 =    A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha );
            a0 =        (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha;
            a1 =    2*( (A-1) - (A+1)*cos(w0)                   );
            a2 =        (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha;

num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);

figure(8);
subplot(2,1,1),semilogx(freq,g1db,'b');grid
title('highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)');
ylabel('dB ');
subplot(2,1,2),semilogx(freq,p1,'b');grid
xlabel('Fréquence (Hz)');
ylabel(' phase (°)')
% % % 
% 
% 
% 
% FYI:   The bilinear transform (with compensation for frequency warping)
% substitutes:
% 
%                                   1         1 - z^-1
%       (normalized)   s  <--  ----------- * ----------
%                               tan(w0/2)     1 + z^-1
% 
%    and makes use of these trig identities:
% 
%                      sin(w0)                               1 - cos(w0)
%       tan(w0/2) = -------------           (tan(w0/2))^2 = -------------
%                    1 + cos(w0)                             1 + cos(w0)
% 
% 
%    resulting in these substitutions:
% 
% 
%                  1 + cos(w0)     1 + 2*z^-1 + z^-2
%       1    <--  ------------- * -------------------
%                  1 + cos(w0)     1 + 2*z^-1 + z^-2
% 
% 
%                  1 + cos(w0)     1 - z^-1
%       s    <--  ------------- * ----------
%                    sin(w0)       1 + z^-1
% 
%                                       1 + cos(w0)     1         -  z^-2
%                                   =  ------------- * -------------------
%                                         sin(w0)       1 + 2*z^-1 + z^-2
% 
% 
%                  1 + cos(w0)     1 - 2*z^-1 + z^-2
%       s^2  <--  ------------- * -------------------
%                  1 - cos(w0)     1 + 2*z^-1 + z^-2
% 
% 
%    The factor:
% 
%                     1 + cos(w0)
%                 -------------------
%                  1 + 2*z^-1 + z^-2
% 
%    is common to all terms in both numerator and denominator, can be factored
%    out, and thus be left out in the substitutions above resulting in:
% 
% 
%                  1 + 2*z^-1 + z^-2
%       1    <--  -------------------
%                     1 + cos(w0)
% 
% 
%                  1         -  z^-2
%       s    <--  -------------------
%                       sin(w0)
% 
% 
%                  1 - 2*z^-1 + z^-2
%       s^2  <--  -------------------
%                     1 - cos(w0)
% 
% 
%    In addition, all terms, numerator and denominator, can be multiplied by a
%    common (sin(w0))^2 factor, finally resulting in these substitutions:
% 
% 
%       1         <--   (1 + 2*z^-1 + z^-2) * (1 - cos(w0))
% 
%       s         <--   (1         -  z^-2) * sin(w0)
% 
%       s^2       <--   (1 - 2*z^-1 + z^-2) * (1 + cos(w0))
% 
%       1 + s^2   <--   2 * (1 - 2*cos(w0)*z^-1 + z^-2)
% 
% 
%    The biquad coefficient formulae above come out after a little
%    simplification.

[632541]

Hello ,

Thank you!

Can you please tell me

How many coefiicients will be generated for 2 order bandpass butterworth digital iir filter(49 to 51)

[Mathieu NOE]

hello again

same amount of coefficients as the order 2 notch

see the demo provided (bandpass is also there)

[632541]

fs= 45e3; %%Sampling frequency 45KHz

order = 2;

t=0:1/fs:1-(1/fs);

%%Input signal

Sig=sin(2*pi*50*t);

[B,A] = butter(order,[(49*2/fs),(51*2/fs)], 'stop'); %%Filter coefficients

Could you please check this..

[Mathieu NOE]

see the doc

you finally generated a 4th order filter because you ignored this :

%   If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an
%   order 2N bandpass filter with passband  W1 < W < W2.

[Mathieu NOE]

corrected code :

fs= 45e3;         %%Sampling frequency 45KHz
order = 2;
t=0:1/fs:1-(1/fs);
%%Input signal 
Sig=sin(2*pi*50*t);
[B,A] = butter(order/2,[(49*2/fs),(51*2/fs)], 'stop');    %%Filter coefficients 
%   If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an
%   order 2N bandpass filter with passband  W1 < W < W2.

[632541]

Thank you so much Mathieu NOE

Understood my mistake in code.

I am using this for notch filter for 50Hz.

[Mathieu NOE]

My pleasure

would you mind accepting my answer ?

tx and good luck for the future