Problem with plotting filter characteristics using 'freqz()' function in matlab

I want to design a high-pass digital filter at cutoff frequency 0.5 Hz. Sampling frequency of the input signal is 500 Hz.

I tried this code:

Fs = 500; % Hz

fc = 0.5; % Hz

[b,a] = butter(12,fc/(Fs/2), "high");

[h,w] = freqz(b,a);

plot(w/pi*Fs/2, abs(h))

After running the code, I get this plot as the magnitude of the filter:

Please help me to fix this problem. Thanks.

Answer 1

[Mathieu NOE]

hello

I have already had that problem with high order filters designed with butter

butter will have numerical problems and the coefficients are wrong

I usually don't go beyong order 5 (always plot the result to double check)

if you really need a 12th order filter, simply put three 4th order filters in serie or use another filter design tool

Order 4 (ok)

Fs = 500;  % Hz
fc = 0.5;  % Hz
[b,a] = butter(4,fc/(Fs/2), "high"); % order 5 MAX !
N = 5000;
[h,w] = freqz(b,a,N);
f = w/pi*Fs/2;
ind = (f>0);
f = f(ind);
H_dB = 20*log10(abs(h));
H_dB = H_dB(ind); 
semilogx(f, H_dB)
xlabel(' Hz ');
ylabel(' gain (dB) ');
title('HPF plot');

[Mathieu NOE]

FYI

some other ways to design digital filters...

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

[Mathieu NOE]

hello

problem solved ?

[Mathieu NOE]

hello

Problem solved ?