Samlping frequency of normalized units in Digital Filter Design block in simulink

I was trying to simulate a digital filter process with Digital Filter Design block in simulink, but I'm confused with the configuration in the block I used. In the Unit tab of Frequency specifications(marked with red line), I chose Nomalized, but I can't determine the sampling frequency. How could digital filter determine the stop and pass frequency without sampling frequncy being set?

[Mathieu NOE]

hi

if you use the normalized frequency , 1 = Nyquist frequency = Fs/2

in this case , Fs does not need to be know as all characteristics are normalized

(ex cut off at 0.1 means 0.1* Fs/2)

depending of the Fs you will choose, the actual cut off will be proportionnal to Fs as explained above

[Raymond Wong]

Hi Mathieu, thanks. I wonder where to set the Fs.

[Mathieu NOE]

Hi Raymond

I assume this filter will be used somewhere in your simulation. So what is your sampling rate for signal processing ?

[Raymond Wong]

I use this filter for a sinusoidal wave filtering with sinusoidal wave noise. I actually didn't set sampling frequency, and it ran well. This is what I'm confused about. It seems I don't need to configure the sampling frequency.

[Mathieu NOE]

Ok

this is because you may not be comfortable with working with normalized frequency .

If you prefer you can do it another way , From the sampling time of your data (recorded or simulated) , and knowing what filter you need, you can directlty generate some commonly used filters. examples below with audio filters examples (digital filters). Otherwise you have the always usefull butterworth (digital) filter (using butter.m function)

% 
%          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.

[Raymond Wong]

Thanks, Mathieu. It's really helpful.

[Mathieu NOE]

hello

you're welcome !! glad it helped

would you please accept my answer (if it was good enough ?)

Answer 1

[Mathieu NOE]

official answer now ...

this is because you may not be comfortable with working with normalized frequency .

If you prefer you can do it another way , From the sampling time of your data (recorded or simulated) , and knowing what filter you need, you can directlty generate some commonly used filters. examples below with audio filters examples (digital filters). Otherwise you have the always usefull butterworth (digital) filter (using butter.m function)

%

% 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.