How to find inverse of filter such that product of filter and its inverse results in a flat spectrum
57 views (last 30 days)
Show older comments
Hello All,
Suppose I have a Low pass filter which would generate a frequency spectrum of two peaks at +x and -x KHz.
Now, I have to create an opposite filter which generates dip at those frequencies so that when I pass the signal through these two filters it gives me a nearly flat spectrum.
I would want to know what is the easiest way to generate/implement this filter in MATLAB.
Thanks for your time in advance.
Answers (2)
Bruno Luong
on 26 Jul 2022
Edited: Bruno Luong
on 26 Jul 2022
If you are working if IIR filter faminy then the answer is trivial yes, since the inverse of the franfert function if a fraction of two polynomials
H(s) = Y(s)/X(s)
is just the
Hi(s) := 1/H(s) = X(s)/Ys(s)
When you apply both filters succesive you'll get the identity. So write in differnce equation recusive form you simply need to swap a and b, the coefficients of the original filter.
The problem is you might get unstable filter. But that is life you cannot do whatever you want.
2 Comments
Star Strider
on 26 Jul 2022
The easiest way to find out is to do the experiment —
Fs = 1E+5;
t = linspace(0, 1E+5, 1E+6)/Fs;
fc = [1 2 4]*1E+4;
s = sum(sin(fc(:)*2*pi*t));
Fn = Fs/2;
L = numel(t);
NFFT = 2^nextpow2(L);
FTs = fft(s,NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTs(Iv))*2)
grid
xlim([0 5E+3])
xlabel('Frequency')
ylabel('Amplitude')
title('Original Signal Spectrum')
fcomb = [[-9 -5 5 9]+1000, [-9 -5 5 9]+2000, [-9 -5 5 9]+4000]; % Design Multiband FIR Notch Filter
mags = [[1 0 1], [0 1], [0 1]];
dev = [[0.5 0.1 0.5], [0.1 0.5], [0.1 0.5]];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
figure
freqz(hh, 1, 2^20, Fs)
set(subplot(2,1,1), 'XLim',[0 5000])
set(subplot(2,1,2), 'XLim',[0 5000])
sgtitle('Filter Bode Plot')
sf = filtfilt(hh,1,s);
figure
subplot(2,1,1)
plot(t,s)
grid
xlim([0 0.0005])
ylim([-2.5 2.5])
ylabel('Amplitude')
title('Original Signal')
subplot(2,1,2)
plot(t,sf)
grid
xlim([0 0.0005])
ylim([-2.5 2.5])
xlabel('Time')
ylabel('Amplitude')
title('Filtered Signal')
Attenuation_dB = 10*log10(sf.^2 / s.^2) % Measure Overall Result
The filter meets the design requirements reasonably well, and produces an acceptable result. This filter uses the default values for a number of otherwise unspecified variables. A more specific design, for example using more options such as provided by designfilt could do even better. (I do not have the DSP System Toolbox, so I am not familiar with its functions. It has other functions and options as well.)
It is likely not possible to completely eliminate the frequencies-of-interest.
.
0 Comments
See Also
Categories
Find more on Digital Filter Design in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!