2
$\begingroup$

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;

enter image description here

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;
$\endgroup$
3
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ I know this is right answer for my original question so, had to edit the question a bit... . $\endgroup$ – Juha P Jul 16 at 12:39
  • $\begingroup$ @JuhaP: Why can't you just compute the complex number $H(e^{j\omega_0})$ and then compute its magnitude? $\endgroup$ – Matt L. Jul 16 at 12:43
  • $\begingroup$ 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). $\endgroup$ – Juha P Jul 16 at 13:35
  • $\begingroup$ @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. $\endgroup$ – Matt L. Jul 16 at 13:37
  • $\begingroup$ @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$. $\endgroup$ – Matt L. Jul 17 at 18:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.