How to create a Butterworth filter from scratch?

I'm trying to write the code for a bandpass Butterworth filter wihtout using the functon butter(), I started creating the low pass filter but I'm having trouble calculating the transfer function

This is what I'm trying to do for the denominator, assuming it's a 4 order filter

syms s;

N=4;

for k=1:N

%Real

w=((2*k+N-1)/(2*N))*pi;

Skr(k) = cos(w);

%Imaginary

Ski(k) = sin(w);

end

Sk=complex(Skr,Ski);

Mag=abs(Sk);

Ang=angle(Sk);

syms s

y1=1;

for k=1:N

y=(s-(Sk(k)));

z=expand(y*y1);

y1=z;

end

[Star Strider]

If you show what you have already done, we might be able to help.

Are you using a specific textbook? If so, which one?

[Marian Vides]

I added some more information!

[Star Strider]

Depending on what you want to do with it, you may need to create a discrete (sampled) version of it to use with sampled data. Use the freqs or freqz functions to detemine if your design does what you want.

[Marian Vides]

I know it’s not doing it cause I’m not obtaining the transfer function with the values of design table for butterworth filters:/

[Rik]

Why did you remove your question? That is very rude.

http://web.archive.org/web/20201028181520/https://webcache.googleusercontent.com/search?q=cache%3AV5bz2IjzDukJ%3Ahttps%3A%2F%2Fuk.mathworks.com%2Fmatlabcentral%2Fanswers%2F625783-how-to-create-a-butterworth-filter-from-scratch+&cd=1&hl=en&ct=clnk&gl=nl

[Rena Berman]

(Answers Dev) Restored edit

[Janisha J]

Derivative of state '1' in block 'ideal_switch_matrixconverter/Unbalanced control block/Computation of phase angle/LPF4' at time 0.00010049999999999999 is not finite. The simulation will be stopped. There may be a singularity in the solution. If not, try reducing the step size (either by reducing the fixed step size or by tightening the error tolerances)

Answer 1

[Mathieu NOE]

hello

please consider this audio digital filter cook book

you should be able to find the filter you need.

all the best

% 
%          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] = mydbode(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] = mydbode(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] = mydbode(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] = mydbode(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] = mydbode(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] = mydbode(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 (°)')
fixfig% % % 
% 
% 
% 
% 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.

[Marian Vides]

Thank you so much! Do you know what kind of filter are they worrying with?

[Mathieu NOE]

sorry I don't understand your question ...

[Marian Vides]

Oh, I mean, do you know if this is a chebyshev filter?

[Mathieu NOE]

no , those formulaes are derivative from analog filters.

so the bandpass filter is butterworth

for chebyshev you can look at the matlab buit in functions

https://www.mathworks.com/help/signal/ref/cheby1.html

https://www.mathworks.com/help/signal/ref/cheby2.html