# Resonance circuit peak filter with unity stopband in MATLAB

I would like to insert a filter with an resonant filter characteristic of a realistic mechanic oscillator. I need to apply it in the time domain, e.g. with a function such as filter or filtfilt. Such, I need filter coefficients fitting my problem. I have found the iirpeak function in MATLAB and implemented it.

fs=22e3;
dur = 2;
t = 1:1/fs:dur-1/fs;
n = randn(1,length(t)) * 1e-4;
x = 2 * sin(2*pi*850*t) + n;

f_res = 800;
[b, a] = iirpeak(f_res/(fs/2),5/(fs/2));
[H, f_iirpeak] = freqz(b,a,fs/10,fs);

loglog(f_iirpeak,abs(H))


The problem is that it damps the signal in the stopband and has unity at resonance. The second problem can be solved by adding a gain factor to all b coefficients. The other problem I was not able to solve, because a mechanic oscillator is described mainly by its $Q$ factor and resonance frequency, and that it does not alter the signal far away from the resonance.

1. How can I alter this filter to have unity stopband characteristic, i.e. shift the filter on the $y$-axis plus one?
2. Is there any other method to get filter coefficients with the wanted characteristic?

EDIT: The solution by Maximilian suggested to see it as an addition of an allpas and the resonator filter. This does part of the trick, but changes the filter characteristic slightly, whether it is "correct" or not. The change in the code is just to add the filter coefficient a to the filter coefficient b: b = a+b.

[b, a] = iirpeak(f_res/(fs/2),100/(fs/2)); % <- here I changed the bw to increase visibility
[H_improved, f_iirpeak] = freqz(b+a,a,fs/10,fs);


Here is a picture of the differences:

So how can this change be avoided? One approach would be to adapt the bandwidth in the initial iirpeak function.

One can get the same bandwidth, if, in the course of coefficients creation, one multiplies the bandwidth (bw) with sqrt(2), for whatever reason.

[b2, a2] = iirpeak(f_res/(fs/2),bw/(fs/2));


EDIT2: The edit from Maximilian correctly suggested that the change in the form of the filter is due to the phase. This will be neglected in the following. Following the comment from Matt L. implemented a peakingEQ filter (see Audio-EQ-Cookbook, link in comment from Matt L.) and compared. The bandwidth function just returns the 3dB bandwidth.

fs=22e3;
dur = 2;
t = 1:1/fs:dur-1/fs;
n = randn(1,length(t)) * 1e-4;
x = 2 * sin(2*pi*850*t) + n;

f_res               = 800;
Q                   = 10;
bw                  = f_res/Q;
fprintf(['Bandwidth: ' num2str(bw) '.\n\n'])

dBgain_pEQ          = 3;
A_pEQ               = 10^(dBgain_pEQ/20);  % (for peaking and shelving EQ filters only)
w0_pEQ              = 2*pi*f_res/fs;       % get angular frequency
alpha_pEQ           = sin(w0_pEQ)/(2*Q);   % (case: Q) <-- (2*Q/sqrt(2))

b_pEQ(1)            =   1 + alpha_pEQ*A_pEQ;
b_pEQ(2)            =  -2*cos(w0_pEQ);
b_pEQ(3)            =   1 - alpha_pEQ*A_pEQ;
a_pEQ(1)            =   1 + alpha_pEQ/A_pEQ;
a_pEQ(2)            =  -2*cos(w0_pEQ);
a_pEQ(3)            =   1 - alpha_pEQ/A_pEQ;

[H_pEQ, f_pEQ]      = freqz(b_pEQ,a_pEQ,fs/10,fs);

[b_ipk, a_ipk]      = iirpeak(f_res/(fs/2),bw/(fs/2));
b_ipk               = b_ipk + a_ipk;
[H_ipk, f_ipk]      = freqz(b_ipk,a_ipk,fs/10,fs);

[bandwidth_ipk] = find_3dB( H_ipk, f_ipk, 0.01 );
[bandwidth_pEQ] = find_3dB( H_pEQ, f_pEQ, 0.01 );

figure
ax(1) = subplot(211);
semilogx(f_pEQ,abs(H_pEQ));
hold on;
semilogx(f_ipk,abs(H_ipk));
semilogx(f_ipk,ones(length(f_ipk),1)*((max(H_pEQ)-min(H_pEQ))/2+1),'k--');
legend('pEQ','ipk','3dB')
set(gca,'XLIM',[10^2 10^4]);
ax(2) = subplot(212);
semilogx(f_pEQ,angle(H_pEQ));
hold on;
semilogx(f_ipk,angle(H_ipk));
legend('pEQ','ipk')


With the hack from Maximilian they look quite alike except for a difference in bandwidth, which can be solved by adding a term to the alpha_pEQ in the script (see comment in source code). This might be due to a power - amplitude confusion at this point, but is neglected at this point.

The second solution with the peaking EQ is a little handier, since it includes a direct gain control.

• Regarding your edit: Certainly, these two functions are not equal to your previously shown code. The peaks in for your code are much narrower. Also, in the frequency domain, there should not be such a change, since it is $X(f)=1+H(f)$ where $X(f)$ has numerator $b+a$ and $H(f)$ has numerator $b$ only. – Maximilian Matthé Feb 8 '17 at 13:56
• You are correct: I forgot to mention that I increased the bandwidth for the picture to 100 Hz to achieve better visibility of the differences. The change though is as it is. The transformation also seems to affect the bandwidth, which is lower in the right, "transformed", picture. – Irenaius Feb 8 '17 at 14:25
• I've added some comments on why this happens to my answer. Though, to solve this issue is not very trivial, as it requires appropriate design of the allpass filter phase. – Maximilian Matthé Feb 8 '17 at 15:22
• You need a so-called peaking EQ filter. You can find corresponding design formulas in the audio-EQ-cookbook. – Matt L. Feb 8 '17 at 15:57

You have designed a bandpass filter, i.e. it does not allow any frequency far way from the resonance to pass. As far as I understand, you want a filter, that does not change the signal in the stopband, but has a resonance at some frequency.

So, you actually want a parallel connection of your bandpass filter and an allpass filter. In Z-Domain, let the bandpass filter be $H(z)$. The allpass filter can be done as $G(z)=1$. Then, the overall filter X(z) is given by

$$X(z)=G(z)+H(z)=1+\frac{b(z)}{a(z)}=\frac{b(z)+a(z)}{a(z)}$$

where $b(z)$ and $a(z)$ are the polynomials of the bandpass filter. So, all you need to do is change your filter coefficients according to:

fs=22e3;
dur = 2;
t = 1:1/fs:dur-1/fs;
n = randn(1,length(t)) * 1e-4;
x = 2 * sin(2*pi*850*t) + n;

f_res = 800;
[b, a] = iirpeak(f_res/(fs/2),5/(fs/2));
[H, f_iirpeak] = freqz(b,a,fs/10,fs);
[H2, f_iirpeak2] = freqz(a+b,a,fs/10,fs);  % here's the change

loglog(f_iirpeak,abs([H H2]))


Here's the program output. The blue curve is the original filter, the green curve is the new filter. You can modify the peak strength by multiplying the b-component of the numerator of the filter, to get different peak heights.

edit: However, as mentioned in the comments, the magnitude response of the modified filter changes slightly compared to the original filter. This is best shown by a filter with a wider passband:

fs=22e3;
dur = 2;
t = 1:1/fs:dur-1/fs;
n = randn(1,length(t)) * 1e-4;
x = 2 * sin(2*pi*850*t) + n;

f_res = 800;
[b, a] = iirpeak(f_res/(fs/2),100/(fs/2));
[b2, a2] = iirpeak(f_res/(fs/2),100/(fs/2));
[H, f_iirpeak] = freqz(b,a,fs/10,fs);
[H2, f_iirpeak2] = freqz(a2+b2,a2,fs/10,fs);

subplot(121);
plot(f_iirpeak,([abs(H) abs(H2)-1]))
ylim([0,1]); xlim([0, 3000]); grid on;
legend({'original', 'modified-1'});

subplot(122);
plot(f_iirpeak, angle(H));
xlim([0, 3000]); grid on;


The green and blue curve do not overlap. The right hand-side figure is the phase of the original filter. The difference in magnitude is due to the phase of the original filter: The overall frequency response of the filter is indeed

$$X(f)=1+H(f),$$

however, this does not mean that

$$|X(f)|=1+|H(f)|$$

since $H(f)$ is a complex value. In order to fix this, one would need to design an allpass filter that has the same / similar phase response as the original filter, and then perform the addition in the Z-Domain. However, then the equations get more complicated.

To prove this idea, consider the following:

[H, f_iirpeak] = freqz(b,a,fs/10,fs);
[H2, f_iirpeak2] = freqz(1j*a2+b2,a2,fs/10,fs);


Here, I designed the allpass filter to have phase $\pi/2$ all the time, i.e. it (roughly) matches the phase for frequencies < 800Hz. Looking at the resulting figure, you see that the filters have similar magnitude response for the lower frequencies. However, it messes up for the higher frequencies, because then the allpass and the bandpass partially cancel out, since the bandpass introduces $-\pi/2$ phase (i.e. -j).

It is up to the designer, if a very accurate response is needed or if a simple "+1" is sufficient.

• This does part of the trick. Thanks. If I am not mistaken it slightly changes the filter characteristic though, better seen in linear scale. – Irenaius Feb 8 '17 at 13:38