hello
below some digital filters used in audio signal processing
hope it helps
%
% Cookbook formulae for audio EQ biquad filter coefficients
% ----------------------------------------------------------------------------
% by Robert Bristow-Johnson <rbj@audioimagination.com>
%
%
% 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.
