Calculation of gain at certain frequency in case of nth order IIR filter

Question

There is this method for to set 0dB gain to be at wanted frequency (fc) (Octave/Matlab example for biquad LPF):
% needed for Octave -------------------------
pkg load signal 
% -------------------------------------------
clf;
% calculate coefficients --------------------
fs = 44100; % sample rate
fc = 700; %Hz
fpi = pi*fc;
wc = 2*fpi;
wc2 = wc*wc;
wc22 = 2*wc2;
k = wc/tan(fpi/fs);
k2 = k*k;
k22 = 2*k2;
wck2 = 2*wc*k;
tmpk = (k2+wc2+wck2);

a0 = 1;
a1 = (-k22+wc22)/tmpk;
a2 = (-wck2+k2+wc2)/tmpk;

b0 = (wc2)/tmpk;
b1 = (wc22)/tmpk;
b2 = (wc2)/tmpk;

b = [b0 b1 b2];
a = [a0 a1 a2];

FLT1 = tf(b, a, 1/fs);

% adjust 0dB @ 1kHz -----------------------------
fc = 1000;  % Hz
w = 2.0*pi*(fc/fs);
num = b0*b0+b1*b1+b2*b2+2.0*(b0*b1+b1*b2)*cos(w)+2.0*b0*b2*cos(2.0*w);
den = 1.0+a1*a1+a2*a2+2.0*(a1+a1*a2)*cos(w)+2.0*a2*cos(2.0*w);
G = sqrt(num/den);

b0 = b0/G;
b1 = b1/G;
b2 = b2/G;

b = [b0 b1 b2]
% ------------------------------------------------

FLT2 = tf(b, a, 1/fs);

% plot
nf = logspace(0, 5, fs/2);
figure(1);
[mag0, pha0] = bode(FLT1,2*pi*nf);
semilogx(nf, 20*log10(abs(mag0)), 'color', 'g', 'linewidth', 2, 'linestyle', '-');
hold on;
[mag, pha] = bode(FLT2,2*pi*nf);
semilogx(nf, 20*log10(abs(mag)), 'color', 'm', 'linewidth', 2, 'linestyle', '-');
legend('LPF', 'LPF 0dB@1kHz', 'location', 'southwest');
xlabel('Hz');ylabel('dB');
axis([1 fs/2 -30 15]);
grid on;

How are formulas for to resolve num and den derived so calculation of G for nth order filter can be done? As for an example for a 4th order filter:
a = [1.00000  -0.61847  -1.09281   0.43519   0.30006];
b = [6.9411e-03   1.1097e-02   5.2508e-03   6.9077e-04  -3.2936e-06];
fc = 1000;  % Hz
w = 2.0*pi*(fc/fs);
num = ...; % ????
den = ...; % ????
G = sqrt(num/den);
b(1) = b(1)/G;
b(2) = b(2)/G;
b(3) = b(3)/G;
b(4) = b(4)/G;
b(5) = b(5)/G;

Matt L. · Answer

You just need to evaluate the transfer function on the unit circle at the frequency of interest:
$$H(e^{jomega_0})=frac{displaystylesum_{k=0}^Nb_ke^{-jkomega_0}}{displaystylesum_{k=0}^Na_ke^{-jkomega_0}}tag{1}$$
and take the magnitude.
For the special values $omega_0=0$ and $omega_0=pi$, Eq. $(1)$ simplifies to
$$H(e^{j0})=frac{displaystylesum_{k=0}^Nb_k}{displaystylesum_{k=0}^Na_k}tag{2}$$
and
$$H(e^{jpi})=frac{displaystylesum_{k=0}^N(-1)^kb_k}{displaystylesum_{k=0}^N(-1)^ka_k}tag{3}$$
respectively.
EDIT: If you want a formula that directly expresses the squared magnitude of $H(e^{jomega})$ then use this:
$$big|H(e^{jomega})big|^2=frac{displaystyle r_b[0]+2sum_{k=1}^Nr_b[k]cos(komega)}{displaystyle r_a[0]+2sum_{k=1}^Nr_a[k]cos(komega)}tag{4}$$
where $r_a[k]$ and $r_b[k]$ are the autocorrelations of the denominator and numerator coefficients, respectively:
$$r_a[k]=a[k]star a[-k]\r_b[k]=b[k]star b[-k]$$
where $star$ denotes convolution.

Calculation of gain at certain frequency in case of nth order IIR filter

One Answer

Add your own answers!

Ask a Question