# Comb filter - attenuation

I have a signal and from 'standard information' I know the predominant frequency is at $1\textrm{ Hz}$ and then there should be three harmonics at lower amplitudes. My signal however has a lower amplitude at $1\textrm{ Hz}$, higher at the first harmonic and then several harmonics, more than what should be present.

I am trying to design a comb filter to rectify this but I am not sure where to start. Any advice from anyone who has some experience with this filter would be great.

• BTW, there is a different convention in the usage of the word "harmonic" than a similar word "overtone". the "first harmonic" is the sinusoidal component with the same frequency as the fundamental (1 Hz in your case). perhaps the "first overtone" is the same as the "second harmonic" (or 2 Hz in your case). Commented Sep 26, 2014 at 17:32

I'm sorta writing this as a separate answer because it's about a sub-topic (how to get an arbitrary delay with continuous-time precision in a discrete-time system).

The simplest way to think about fractional-sample interpolation, in my opinion, is just to consider the reconstruction half of the Nyquist-Shannon-Whittaker-Kotelnikov-whomever sampling and reconstruction theorem. If a continuous-time signal $$x(t)$$ is sampled sufficiently at sample rate $$f_s=\frac{1}{T}$$ ($$T$$ is the sampling period)

$$x[n] \triangleq x(nT) \quad n \in \mathbb{Z}$$

then $$x(t)$$ is reconstructed from the samples $$x[n]$$ with

\begin{align} x(t) & = \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \frac{\sin\left(\pi \tfrac{t-nT}{T} \right)}{\pi \tfrac{t-nT}{T}}\\ \\ & = \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \operatorname{sinc}\left( \tfrac{t-nT}{T} \right) \\ \\ & = \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \operatorname{sinc}\left( \tfrac{t}{T} - n \right) \\ \end{align}

where

$$\operatorname{sinc}(u) \triangleq \begin{cases} \frac{\sin(\pi u)}{\pi u} \qquad & u \ne 0 \\ 1 & u = 0 \end{cases}$$

Note that $$\operatorname{sinc}(u)$$ is an even-symmetry function

$$\operatorname{sinc}(-u) = \operatorname{sinc}(u)$$

Now, without loss of generality, let's say that $$t=0$$ is the present instant of time (which means that $$x[0]$$ is the most current sample of $$x(t)$$, all past samples, $$x[n]$$, have negative indices, and none of the future samples with positive indices exist) and we want to know what $$x(t)$$ is at some delayed time $$x(-D)$$ where delay $$D>0$$.

$$x(-D) = \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \operatorname{sinc}\left( \tfrac{-D}{T} - n \right)$$

First thing to do is split the normalized delay, $$\frac{D}{T}$$, which is dimensionless, into integer, $$i_D$$, and fractional, $$f_D$$ parts:

$$\tfrac{D}{T} = i_D + f_D$$

where

$$i_D \triangleq \left\lfloor \tfrac{D}{T} \right\rfloor = \operatorname{floor}\left( \tfrac{D}{T} \right)$$

$$f_D = \tfrac{D}{T}- i_D = \tfrac{D}{T} - \left\lfloor \tfrac{D}{T} \right\rfloor$$

Note that $$\tfrac{D}{T} - 1 < i_D \le \tfrac{D}{T} < i_D + 1$$

and $$0 \le f_D < 1$$ .

Then $$x(-D)$$ is reconstructed from the samples $$x[n]$$ with

\begin{align} x(-D) & = \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \operatorname{sinc}\left( \tfrac{-D}{T} - n \right) \\ \\ & = \sum\limits_{n=-\infty}^{+\infty} x[n] \cdot \operatorname{sinc}\left( -(i_D + f_D) - n \right) \\ \end{align}

To make this summation finite, we have to toss out all of the samples of $$x[n]$$ but those closest to $$x[-i_D]$$. This is essentially what windowing is about. If we're keeping the $$L$$ (where $$L$$ is even) closest samples on the left and right of the interpolated value at a time precisely $$\tfrac{D}{T}$$ samples in the past:

\begin{align} x(-D) & \approx \sum\limits_{n=-i_D - L/2}^{-i_D + L/2 - 1} x[n] \cdot \operatorname{sinc}\left( -(i_D + f_D) - n \right) w(n + i_D + f_D) \\ \\ & = \sum\limits_{n= -L/2}^{L/2 - 1} x[n - i_D] \cdot \operatorname{sinc}\left( -f_D - n \right) w(n + f_D) \\ \\ & = \sum\limits_{n= -L/2}^{L/2 - 1} x[n - i_D] \cdot \operatorname{sinc}( f_D + n ) w(n + f_D) \\ \end{align}

where $$w(t)$$ is some window function. If it were a Hamming Window, it would be:

$$w(t) = \begin{cases} 0.54 \ + \ 0.46 \cdot \cos\left(\pi \frac{t}{L/2} \right) \qquad & |t| \le \tfrac{L}{2} \\ 0 & |t| > \tfrac{L}{2} \\ \end{cases}$$

If it were a Kaiser window it would be

$$w(t) = \begin{cases} \frac{1}{I_0(\beta)} \, I_0\left(\beta \sqrt{1 - \left(\frac{t}{L/2}\right)^2 } \right) \qquad & |t| \le \tfrac{L}{2} \\ 0 & |t| > \tfrac{L}{2} \\ \end{cases}$$

where

$$I_0(u) \triangleq \sum\limits_{k=0}^{\infty} \frac{(-1)^k \big( \tfrac{u}{2} \big)^{2k}}{(k!)^2}$$

is the zeroth-order Bessel function of the first kind. With the Kaiser window, knowing the length $$L$$ (like $$L$$ might be 16 samples, or 32 if you want it really good) you can choose $$\beta$$ to tradeoff between the stopband attenuation and transition bandwidth of the reconstruction brickwall filter. A good value of $$\beta$$ is 6 or 7, to get you about 63 dB or 72 dB of stopband attenuation to make your reconstruction brick wall filter have very solid bricks.

Again, assuming for simplicity that $$x[0]$$ is your current sample, if the normalized delay is at least as long as $$\tfrac{L}{2}-1$$ samples, then

\begin{align} y_D(0) & \triangleq x(0-D) \\ & \approx \sum\limits_{n= -L/2}^{L/2 - 1} x[n - i_D] \cdot \operatorname{sinc}\left( n + f_D \right) w(n + f_D) \\ \\ & = \sum\limits_{n= -L}^{-1} x[n + \tfrac{L}{2} - i_D] \cdot \operatorname{sinc}\left(n + f_D + \tfrac{L}{2} \right) w(n + f_D + \tfrac{L}{2}) \\ \\ & = \sum\limits_{n= -L+1}^{0} x[n + \tfrac{L}{2} - 1 - i_D] \cdot \operatorname{sinc}\left(n + f_D + \tfrac{L}{2}-1 \right) w(n + f_D + \tfrac{L}{2}-1) \\ \\ & = \sum\limits_{n= -L+1}^{0} x[n - (i_D-\tfrac{L}{2}+1)] \cdot h_{f_D}[0-n] \\ \\ & = \ \sum\limits_{n=0}^{L-1} \ x[(\tfrac{L}{2}-1-i_D)-n] \cdot h_{f_D}[n] \\ \end{align}

That's a simple dot product of the $$L$$ samples from $$x[-i_D+\tfrac{L}{2}-1]$$ to $$x[-i_D-\tfrac{L}{2}]$$ with a set of FIR tap coefficients $$h_{f_D}[n]$$, (remember, here we're defining the most current sample as $$x[0]$$ and we know $$x[0], x[-1], x[-2], ...$$ ) and where $$i_D$$ is the integer part of the delay and $$i_D \ge \tfrac{L}{2}-1$$.

The FIR tap coefficients:

$$h_{f_D}[n] = \operatorname{sinc}\left(n + 1 - f_D - \tfrac{L}{2} \right) w(n + 1 - f_D - \tfrac{L}{2}) \qquad 0 \le n \le L-1$$

are dependent only on the fractional part $$f_D$$ of the delay.

So this is how you implement a precision delay of $$\frac{D}{T}$$ samples, or $$D$$ seconds of delay (assuming $$T$$ is expressed in seconds). In the continuous-time domain, this would be a delay of $$D$$ seconds with a Laplace transfer function of

$$H_D(s) \triangleq \frac{\mathscr{L}\{y_D(t)\}}{\mathscr{L}\{x(t)\}} = \frac{Y_D(s)}{X(s)} \approx e^{-sD}$$

which in the other partial answer would replace $$z^{-1}$$ in your Direct Form II filter structure to get you a comb filter with an IIR or FIR design of a digital low-pass filter $$H(z)$$.

• I feel silly for bumping this, but I'm curious about something mentioned in Robert's answer: "whatever frequency response you design for a low-pass filter (either FIR or IIR) can be used to design the shape of the teeth of the comb filter". May I please ask how this can be done? Any recommended reading? I found this thread as I was searching for comb filter design information, and this is precisely what I need. (I apologize if this is an inappropriate form of a comment, I'm new to this website).
– Lior
Commented Dec 8, 2019 at 17:42
• @Lior, the description of how the frequency response of a multitap comb filter is determined is in the other answer above. this answer is how to get precisely a delay of amount $D$. Commented Dec 9, 2019 at 15:44
• I think I was previously off by 1 in the finite limits to some of the summations. I think I fixed it correctly. I wouldn't mind if someone else double-checks this bookkeeping. Remember the delay $D>0$ and I think $i_D \ge \frac{L}{2} - 1$ so that $$\tfrac{L}{2} - 1 \le i_D + f_D = \tfrac{D}{T}$$ Commented Apr 29, 2023 at 21:56
• It turns out that three approaches to fractional sample delay using FIR filters are equivalent: Lagrangian interpolation, approximation with a maximally flat error, and windowing the ${\tt sinc}$ function using a binomial window. See here for details.
– Peter K.
Commented Feb 12 at 12:22
• @PeterK. that paper looks scary to me. (my brain is flabby because I'm old.) I have done polynomial interpolation but I hadn't bother for more than 3rd order and interpolating in the middle of 4 adjacent samples. When it's using 32 adjacent samples, I've done it by designing a LPF with cutoff of Nyquist/R where R is the upsample rate. I found that R=512 is plenty good enough if you linearly interpolate between adjacent upsampled samples. So that's 512 sets of 32 coefficients each, so we design an FIR filter of length 32 $\times$ 512 - 1 (and padding it with one zero to be an even 16K. Commented Feb 12 at 22:44

Just something to add, there are two kinds of canonical comb filters. there are those made with only feedforward paths and those that have feedback. In fact, you can have additional long delay elements of length $$\frac{1}{f_0}$$ where $$f_0$$ is the frequency of the spacing between the teeth of the comb.

In fact, a comb filter can be designed as generally as an FIR or IIR filter, but with the unit delay elements replaced by delay elements with delay $$\frac{1}{f_0}$$. Whatever frequency response you design for a low-pass filter (either FIR or IIR) can be used to design the shape of the teeth of the comb filter.

If you use a delay of half that above, a delay of $$\frac{1}{2 f_0}$$, your comb filter will be tuned to only the even-numbered harmonics. In addition, you can instead design for a prototype high-pass filter (again with delay elements of half length $$\frac{1}{2 f_0}$$) and your comb filter will pass only the odd-numbered harmonics.

In continuous-time, using the Direct Form II, it's

$$v(t) = x(t) - \sum\limits_{n=1}^N a_n v(t - n/f_0)$$

$$y(t) = \sum\limits_{n=0}^N b_n v(t - n/f_0)$$

The transfer function is

$$H(s) = \frac{Y(s)}{X(s)} = \frac{\sum\limits_{n=0}^N b_n e^{-n s/f_0}}{1+ \sum\limits_{n=1}^N a_n e^{-n s/f_0}}$$

or frequency response

$$H(j2 \pi f) = \frac{Y(j2 \pi f)}{X(j2 \pi f)} = \frac{\sum\limits_{n=0}^N b_n e^{-j n 2 \pi f/f_0}}{1+ \sum\limits_{n=1}^N a_n e^{-j n 2 \pi f/f_0}}$$

We can discuss how you might design your LPF prototype that depicts the shape of the teeth. And also, maybe we should discuss how to do a precision delay that has sub-sample precision so that $$f_0$$ can be precisely tuned.

In the other answer, the precision delay $$D$$ would be integer multiples of $$\frac{1}{f_0}$$ in seconds. and the normalized (to the sample period) precision delay (with integer and fractional part) is

$$\frac{D}{T} = f_\text{s} D = i_D + f_D = \frac{n f_\text{s}}{f_0} \qquad n=1,2,... N$$

(Don't conflate or confuse the symbols $$f_\text{s}$$, the sample rate, with $$f_0$$, the fundamental frequency or "pitch", or the fractional delay $$f_D$$.)

Comb filters involve adding a delayed version of a signal to itself, they are generally used for musical effects processing and are unlikely to be of much help to you here. You may have been slightly misled by the frequency response which appears as a series of peaks.

To increase the amplitude around a particular frequency, for example your dominant 1Hz, you can consider a peaking filter. To eliminate the additional harmonics you may consider band stop filters. Bear in mind you can apply multiple independent filters to the signal in sequence to achieve the changes in frequency response that you desire.

You probably could achieve the desired result with a single filter but you have a set of conditions which I suspect would lead to a complicated and expensive filter, so multiple passes will be an easier approach.

• Could you offer some advice re using the peak filter? I am trying to use a peak filter to amplify my peak but for some reason I get an error message ??? No appropriate method, property, or field FFTsignal2 for class dfilt.df2sos. I want to amplify a peak at 1Hz (or 1/30 if the frequency is normalised) and so I have used the following; N = 8; F0 = (1/30); BW = 0.1; f = fdesign.peak('N,F0,BW',8,(1/30),.05,80); Then I use: (tbh I am not sure what to use here because ideally I would like a nice tight peak) Hd = design(f,'butter'); then to apply I use FFTfilteredSignal = filter(Hd.FFTsignal2); Commented Aug 28, 2014 at 0:37
• how are you doing your peak filter? how are you doing this whole thing? are you doing this entirely in the frequency domain with the FFT? i had thought originally you were looking for time-domain comb filter design. if you're doing this in the frequency domain with the FFT, you might have to worry about overlap-add (OLA) or overlap-save (OLS) and the shape of the FIR filter. you cannot just choose any arbitrary comb filter shape you want, but it has to correspond to an FIR filter of known length so that you can use OLA or OLS. it's easier implementing a comb filter in the time-domain. Commented Sep 26, 2014 at 20:17