Comb filter - attenuation

Question

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.

Answer 1

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)$.

Answer 2

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$.)

Answer 3

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.