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
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.
- How can I alter this filter to have unity stopband characteristic, i.e. shift the filter on the $y$-axis plus one?
- Is there any other method to get filter coefficients with the wanted characteristic?
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 = 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
One can get the same bandwidth, if, in the course of coefficients creation, one multiplies the bandwidth (
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') linkaxes(ax(:),'x');
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.