Comb Filter Design and Implementation using MATLAB

Chinmaykshirsagar
8 min readJan 4, 2022

What are Comb Filters?

Comb filters are fundamental components of digital sound effects. A comb filter is used in the acoustic echo modeling. The two main comb-filter types, feedforward and feedback, are discussed in this article, along with a frequency-response analysis.

Types of Comb Filters

There are mainly two major types of comb filtering i.e Feed Forward and Feedback.

Feedforward comb filters

Figure 1 depicts the feedforward comb filter. Around the delay line, the direct signal “feeds forward.” The direct and delayed signals are combined into a linear output.

Fig. 1 Feedforward comb filter

The feedforward comb filter’s difference equation is,

The feedforward comb filter is a form of FIR filter, as we can see. It’s also a unique type of Tapped Delay Line (TDL). By placing b0=1 and bM=g in the feedforward comb filter, the echo simulator may be implemented. It’s a calculation based physical model of a single discrete echo, in other words. This is among the most basic applications of signal processing elements in acoustic modelling. The feedforward comb filter simulates the superposition of a “direct signal” b0x(n) and an attenuated, delayed signal bMx(n-M), where the attenuation (by |bM|<1) is caused by “air absorption” and/or spherical spreading losses, and the lag is caused by acoustic transmission over a range of cMT metres, where T is the sampling period in seconds and c is sound speed. When the simulated transmission latency must be more precise than the closest integer value of samples M, delay-line interpolation must be utilised. When air absorption has to be more precisely replicated, the constant attenuation factor bM can be substituted by a linear, time-invariant filter G(z) that gives a distinct attenuation at each frequency. G(z) has a lowpass quality due to the physics of air absorption.

Feedback comb filters

As demonstrated in Fig.2 (drawn in “direct form 2”), the feedback comb filter employs feedback instead of a feedforward signal.

Fig 2. Feedback comb filters

In “direct form 1”, a difference equation defining the feedback comb filter may be expressed as,

Because there is feedback from the delayed output to the input, the feedback comb filter is a particular instance of an Infinite Impulse Response (IIR) (“recursive”) digital filter. The feedback comb filter is a computational physical model of a succession of echoes that decay exponentially and are spaced equally in time. Take, for example, the unique situation,

is a computational model of an ideal plane wave bouncing back and forth between two parallel walls, where g denotes complete round-trip attenuation (two wall-to-wall traversals, including two reflections).The feedback coefficient a M must be less than 1 in magnitude for stability, i.e. |aM| > 1. If a |aM| > 1 is not set, each echo will be louder than the one before it, resulting in a never-ending, expanding chain of echoes. When the output signal arrives from the end of the delay line rather than the commencement, the difference equation changes to,

This output option simply delays the signal by M samples.

Filtered Feedback Comb Filters

Rather than only a feedback gain, the filtered-feedback comb filter (FFBCF) incorporates filtered feedback. The filtered-feedback comb filter’s transfer function may be represented as Hl(z) where the feedback-filter transfer function is denoted by

The FFBCF may be viewed theoretically a particular instance of the generic all pole transfer function in which the first M-1 denominator coefficients are bound to be zero when Hl(z) is a causal filter:

The FFBCF is computationally more efficient than other, more general-purpose IIR filter architectures due to the “sparseness” of the filter coefficients.

A feedback-comb-physical filter’s understanding is that it simulates a plane-wave bouncing back and forth between two walls. The addition of a lowpass filter to the feedback loop replicates frequency-dependent losses that occur during a propagation round-trip, as they do in actual rooms.

Air absorption and the coefficient of absorption at each wall are the two primary physical drivers of plane-wave attenuation. Scattering causes additional “losses” for plane waves in actual rooms. (When a plane wave collides with something other than a wall, it reflects in a variety of directions.) Textured wall surfaces are a common scatterer in performance venues. Reflections from such walls are commonly approximated using a specular and diffuse component in ray-tracing simulations. Wavelengths larger than the texture grain size reflect specularly (with some attenuation owing to wall motion), but wavelengths on the order of or less than the texture grain size scatter in numerous directions, contributing to the diffuse component of reflection.

In computer music, the filtered-feedback comb filter has several uses. Schroeder, it appears, was the first to propose artificial reverberation, and Moorer was the first to put it into practice. The FFBCF may be thought of as a transfer-function physical-model of a vibrating string when the Karplus-Strong algorithm is interpreted physically. FFBCFs are frequently created as computationally optimised equivalent types depending on a basic waveguide model constructed in terms of bidirectional delay-lines (“digital waveguides”) in digital waveguide modelling of string and wind instruments.

The amplitude-response of the feedback-filter Hl(z) must be smaller than 1 at all frequencies for stability, i.e.,

Applications of comb filtering:

  • Cascaded integrator–comb (CIC) filters are often employed for anti-aliasing in discrete-time systems during interpolation and decimation processes that modify the sample rate.
  • Dot crawl is reduced by 2D and 3D comb filters included in hardware (and rarely software) in Phase Alternating Line (PAL) and National Television System Committee (NTSC) analogue television decoders.
  • Delay, flanging, physical modelling synthesis, and digital waveguide synthesis are all examples of audio signal processing. A comb filter may simulate the impact of acoustic standing waves in a cylindrical cavity or a vibrating string if the delay is adjusted to a few milliseconds.
  • In astronomy, the astro-comb promises to roughly a hundredfold enhance the accuracy of present spectrographs.

MATLAB Implementation:

A. Forward Feed Comb Filter Amplitude and Phase Response

  • All the output responses as we can see look like comb.
  • This type adds echo to the input signal
  • The Z domain equation used for the feedforward comb filter is as follows:
  1. alpha=0.8 and N0=11
w=-pi:0.01:pi;
num=[1 zeros(1,10) 0.8];
den=[1];
h=freqz(num,den,w);
subplot(2,1,1);
plot(w/pi,abs(h),'b');
ylabel('Magnitude');
subplot(2,1,2);
plot(w/pi,angle(h));
ylabel('Phase');

Freqz generates the frequency response of the transfer function given by num, den over the interval w. abs(h) and angle(h) represents the amplitude and phase response of the Comb filter respectively. =0.8 and N0=11

2. alpha=1.2 and N0=9

w=-pi:0.01:pi;
num=[1 zeros(1,8) 1.2];
den=[1];
h=freqz(num,den,w);
subplot(2,1,1);
plot(w/pi,abs(h),'b');
ylabel('Magnitude');
subplot(2,1,2);
plot(w/pi,angle(h));
ylabel('Phase');

The above snippet is more or less same as the above, the only difference is the value of and N0. The value of is increased and that of N0 is reduced. The phase response becomes sharper.

3. alpha=0.5 and N0=11

num=[1 zeros(1,10) 0.5];
den=[1];
h=freqz(num,den,w);

In this we have reduced and N0 remains the same. Rest part of the code i.e the plotting of Amplitude and Phase spectrum remains the same as above. Here w is also referred to as normalized frequency.

B. Feedback Comb Filter Amplitude and Phase Response

  • This type removes echo present in the input signal
  • The Z domain equation used for the feedforward comb filter is as follows:
  1. alpha=0.8 and N0=11
w=-pi:0.01:pi;
den=[1 zeros(1,10) 0.8];
num=[1];
h=freqz(num,den,w);
subplot(2,1,1);
plot(w/pi,abs(h),'b');
ylabel('Magnitude');
subplot(2,1,2);
plot(w/pi,angle(h));
ylabel('Phase');

The snippet is more or less the same as above except the fact that the transfer function in this case is reciprocal of that in feedforward. In the plot function w is divided by pi in order to normalize the output between 0–1.

2. alpha=1.2 and N0=11

den=[1 zeros(1,10) 1.2];
num=[1];
h=freqz(num,den,w);

In this we have increased and N0 remains the same. The function abs(h) basically plots the amplitude response of h in a normalized range of w.

3. alpha=0.5 and N0=11

den=[1 zeros(1,10) 0.5];
num=[1];
h=freqz(num,den,w);

In this we have reduced and N0 remains the same. Freqz function returns a vector with transfer function and w as interval.

C. Impulse response of FeedBack comb filter:

N=1000;
den=[1,zeros(1,10),0.9];
num=[1];
n=0:N-1;
x=(n==0);
y=filter(num, den, x);
stem(n,y);

H(z) = 1/(1+ 0.9z-11) is the transfer function for a comb filter whose impulse response is generated. Filter function in MATLAB is used to filter out x from the transfer function given by num and den.

D. Impulse response of FeedForward Comb Filter

N=1000;
num=[1,zeros(1,10),0.9];
den=[1];
n=0:N-1;
x=(n==0);
y=filter(num, den, x);
stem(n,y);

H(z) =1+ 0.9z-11 is the transfer function for a comb filter whose impulse response is generated.

E. Echo Removal using comb filter

Step1: Add White noise to a signal

Step2: Remove it using filter

N=5000;
x=randn(1,N);
sound(x,8000)
w=-pi:0.01:pi;
den=[1,zeros(1,10),0.9];
num=[1];
h=freqz(num,den,w);
subplot(2,1,1);
plot(w/pi,abs(h),'b');
ylabel('Magnitude');
subplot(2,1,2);
plot(w/pi,angle(h));
ylabel('Phase');

y=filter(num,den,x);
sound(y,8000)

Here we have introduced random noise in the input signal and removed it using a comb filter. sound() function is used to give aquostic output.

F. Pole Zero plot of Feed Forward and FeedBack comb filter

w=-pi:0.01:pi;
num=[1 zeros(1,10) 0.8];
den=[1];
h=freqz(num,den,w);

[num,den] = eqtflength(num,den);
[z,p,k] = tf2zp(num,den);

fvtool(b,a,'polezero');
text(real(z)+.1,imag(z),'Zero');
text(real(p)+.1,imag(p),'Pole');

First the transfer function is defined for the comb filter and then we just call the tf2zp function to plot pole-zero plot. Fvtool is used to visualize the plot.

References:

  1. https://ccrma.stanford.edu/~jos/pasp/Comb_Filters.html
  2. https://youtu.be/yEVoMzHqcoQ
  3. https://www.mathworks.com/help/dsp/ref/dsp.allpassfilter.fvtool.html
  4. https://www.mathworks.com/help/matlab/ref/filter.html
  5. https://www.mathworks.com/help/matlab/ref/stem.html
  6. https://www.mathworks.com/help/signal/ref/tf2zp.html
  7. https://www.mathworks.com/help/signal/ref/freqz.html
  8. https://en.wikipedia.org/wiki/Comb_filter

Authors:

Neeraja Khire

Mahesh Kinge

Chinmay Kshirsagar

Parth Kudal

--

--