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;

Answer 1

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.