Forums

Designing an IIR comb (peak) filter

Started by bobby_k 2 weeks ago9 replieslatest reply 2 weeks ago100 views

Hi there,

I'm trying to design an IIR comb filter to keep only 10 Hz harmonics in a noisy signal which is periodic at 10 Hz for noise rejection.  The comb filter should then have peaks at 10 Hz, 20 Hz, 30 Hz ... (probably up to 200 Hz), and low amplitude response for other frequencies. 

This is the code I've been working on in MATLAB:  IIR_comb.m

When I run the code, the figures I get look like complete mess. 

I'm not too concerned about performance trade-offs at the moment - just trying to get a sensible response for now.

Cheers

[ - ]
Reply by neiroberSeptember 3, 2019

bobby_k,

Your sample rate is 15 kHz.  So you can decimate by 4 to 3750 kHz.  Then this code gives the response you are looking for, I think.

regards,

Neil


fs= 3750;

 b= [1 zeros(1,374) 1]/2;

[h,f]= freqz(b,1,8192,fs);

plot(f,abs(h)),grid

axis([0 200 0 1]);

xlabel('Hz')


comb1_45714.png

[ - ]
Reply by napiermSeptember 4, 2019

Nice thing is this can be implemented with 1 adder and a memory.


[ - ]
Reply by bobby_kSeptember 3, 2019

Hi Neil,

Thanks for the help.

Could you please explain a few things with your design process?

1) Why did you downsample by a factor of 4?

2) It would seem that your filter is an FIR filter (since denominator in freqz is 1), how did you select these coefficients?

3) What's the significance of 8192?

Cheers

[ - ]
Reply by neiroberSeptember 3, 2019

bobby_k:

1.  Since you said your band of interest was 0 -200 Hz, it is ok to decimate by 4 to give fs= 3750 Hz.  You don't have to decimate, but if you don't, the comb filter is 1501 taps instead of 376 taps.

To decimate, you would need to use a decimation filter, not just downsample.  For example, a very simple decimation-by-4 filter would have coefficients

 [1 1 1 1]/4.  Then you downsample by 4 after filtering.  As I noted, decimation allows shortening the length of the comb filter by a factor of 4.

2.  Yes, it is FIR.  See this article for the theory:

https://www.dsprelated.com/freebooks/pasp/Comb_Filters.html

3.  8192 is the number of frequency bins computed by the matlab function freqz - freqz just computes half the points of a 16384-point fft, and computes a frequency vector f.

regards,

Neil


[ - ]
Reply by kazSeptember 3, 2019

Without going into details of your code... your poleAngle and poleCoeff is updated inside loop for each value of i and takes in last value. Are you aware of that

[ - ]
Reply by bobby_kSeptember 3, 2019

Hi kaz,

Yes, that's correct. I'm placing a pair of conjugate poles per iteration of the for loop (at increments of 10 Hz) and multiplying the transfer functions each time to yield the overall transfer function. 

cheers

[ - ]
Reply by Rick LyonsSeptember 3, 2019

Hi bobby_k

I tried to run your code but it contains calls to functions 'tf()' and 'tfdata()" that I do not have.


[ - ]
Reply by rbjSeptember 3, 2019

it turns out that *every* digital filter is a comb filter (of sorts).  if the digital filter is a LPF, you will see comb teeth at every integer multiple of the sample rate $n \dcot f_\text{s} = \frac{n}{T_\text{s}}$.  if the digital filter is a HPF, you will see teeth halfway between every integer multiple of the sample rate.

now change your unit sample delay $z^-1 = e^{-s T_\text{s}}$ to a longer delay that i will express only in the continuous $s$-domain: $e^{-sT}$ where $\frac{n}{T}$ are the frequency you want to tune your teeth of the comb to.

so you need to make a precision delay to tune it to exactly the fundamental frequency you want, and then you can design the teeth of the comb filter just like you design any old LPF.

and i still haven't figured out how to do $\LaTeX$ here at dsprelated.

[ - ]
Reply by Tim WescottSeptember 3, 2019

$$ H(z) = \left(1 - d \right)\frac{z^N}{z^N - d} $$

Choose $$ N = f_s / 10\mathrm{Hz} $$ and make d as small as you want your peaks to be pronounced.