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

There is this method for to set 0dB gain to be at wanted frequency (fc) (Octave/Matlab example for biquad LPF):

% needed for Octave -------------------------
% -------------------------------------------
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;


You just need to evaluate the transfer function on the unit circle at the frequency of interest:

$$H(e^{j\omega_0})=\frac{\displaystyle\sum_{k=0}^Nb_ke^{-jk\omega_0}}{\displaystyle\sum_{k=0}^Na_ke^{-jk\omega_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{\displaystyle\sum_{k=0}^Nb_k}{\displaystyle\sum_{k=0}^Na_k}\tag{2}$$

and

$$H(e^{j\pi})=\frac{\displaystyle\sum_{k=0}^N(-1)^kb_k}{\displaystyle\sum_{k=0}^N(-1)^ka_k}\tag{3}$$

respectively.

EDIT: If you want a formula that directly expresses the squared magnitude of $$H(e^{j\omega})$$ then use this:

$$\big|H(e^{j\omega})\big|^2=\frac{\displaystyle r_b[0]+2\sum_{k=1}^Nr_b[k]\cos(k\omega)}{\displaystyle r_a[0]+2\sum_{k=1}^Nr_a[k]\cos(k\omega)}\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.

• I know this is right answer for my original question so, had to edit the question a bit... . – Juha P Jul 16 at 12:39
• @JuhaP: Why can't you just compute the complex number $H(e^{j\omega_0})$ and then compute its magnitude? – Matt L. Jul 16 at 12:43
• wouldn't it need some complex math then(?) which I don't comprehend well enough and also don't use in this implementation of weighting filter (using Ruby language). – Juha P Jul 16 at 13:35
• @JuhaP: No, you can just compute the real and imaginary parts separately, square and add them, and you obtain the squared magnitude of the spectrum evaluated at the given frequency. – Matt L. Jul 16 at 13:37
• @JuhaP: I've added the formula for directly computing the squared magnitude of the frequency response. It's a generalization of what you used for $N=2$. – Matt L. Jul 17 at 18:12