C++ Reverb – DSP Final Project

Overview

This past year (Fall 2018) I took a DSP (digital signal processing) class at Cedarville University. At the end of the class, as a final project we had to implement something “cool” using the techniques we learned in the rest of class on a Raspberry Pi. This was supposed to be done using the Raspberry Pi package for Simulink Coder.

For my project, I decided to implement a digital reverb unit; as I thought it would be interesting. I also thought it would be fun to add reverb to my music library. The reverb unit was first implemented in Matlab and Simulink for testing on a desktop, then the filters were ported to C++ so that they could be directly implemented on the Pi without using Simulink coder.

A sample input and output (along with all the code needed for the C++ implementation are available on github). The files Invention_no1.mp3 and Invention_no1_reverb.mp3 are the input and output files respectively.

This post does not go into the math behind digital signal processing, mostly because I did not feel like writing about it. I might eventually, but I probably won’t.

30 Second DSP Review

Here is just a quick review of some basic digital signal processing concepts that might make the rest of the post make more sense (in my dreams).

Basically, digital signal processing is what happens when you take an analog signal and sample it (this puts it into the discrete time world). When you do this, some interesting things happen. Namely, you can perfectly recreate the signal up to 1/2 the sampling rate used to sample the signal (thank you Mr. Nyquist). Once the signal has been sampled, filters are needed to do things to the signal. These filters are implemented using a difference equation (the discrete time cousin of the differential equation). An example of one is shown below.

y[n] + b_1*y[n-1] + b_2*y[n-2] = a_1*x[n] + a_2*x[n-1] + ... + a_n*x[n-M]

These differences aproximate derivatives (think of the definition of a derivative, then remove the limit and the x – Δx from the denominator)

derivative1derivative2

Now if you let f be discrete (only valid on integer inputs) let and let each “index” of the function be incremented by Δx, you can see how the first derivative can be aproximated by f(x) – f(x-1).

These difference This can be moved into the discrete frequency domain using some math. (Something something discrete fourier transform. Something something Z-transform). Once in the discrete frequency domain, the frequency response of the filter can be evaluated. This is pretty much equivelent to how it would look in the continous domain (except that the filter response repeats every multiple of the sampling rate.

Basically if you really care, look up an online DSP course, or become an electrical engineer 🙂 (MIT open-courseware)

Reverb Unit

Because I do not have a large knowledge of acoustics, I had to do some research on how the effect is implemented digitally. I primarily used a page about the Schroeder digital reverb unit from Stanford, and a page about digital effects from Spin Semiconductors as my sources of information.

An algorithmic digital reverb unit (as opposed to sampling a room) is based on passing a signal through a series of all-pass filters. These filters do not change the magnitude of the signal passed through them, but only change its phase in a frequency dependent way (see the figure below). The all-pass filter models the effect of a sound wave bouncing around a room from an instrument back to the ear, the relative amplitudes of each of the frequencies probably will not change much, but their relative delays do.

Allpass
All-pass filter frequency and phase response

The phase of the signal can be thought of as how much the signal has been delayed. The large spike in frequency between -pi and pi in the phase plot is due to the fact that Matlab (i.e. the arctan function) cannot tell the difference between positive and negative 180 degrees (pi radians), so anything more than -180 degrees is made to be positive. You can also see that the phase shift on the signal is frequency dependent, so a high frequency signal gets delayed more than a low frequency one.

Allpass-Imp1
All-pass filter block diagram

The image above shows an implementation for an all-pass filter from the Spin Semiconductors website. The circles with pluses are adders and with X’s are multipliers. By changing the gains on the multipliers, the shape of the frequency response can be changed.

Reference Designs

The Stanford page shows an implementation of a Schroeder reverb unit (shown below).

Schroeder_Reverb
Schroeder reverb unit

The blocks labeled AP are all pass filters with a the top number being the gain on the feedback/feed-forward taps (the multiplier blocks in the second image), and the bottom number being the delay of the filter. The blocks labeled FBCF are feedback comb filters, the numbers refer to same thing as in the all-pass filters. The outputs from these blocks are fed into a mixing matrix to produce multi-channel audio outputs.

The article describes the method used to pick the delay values for the all-pass filters
N_i = 100ms/(3^i) * F_s
Where ‘i’ is a positive integer, N_i is the delay for the i’th all-pass filter, and F_s is the sample frequency of the system. (Sorry I don’t have a good way to put equations into wordpress).

The Spin Semiconductors article describes a different reverb topology that chains multiple all-pass filters into a feedback loop (shown below).

SpinSemi_Reverb
Spin Semiconductors reverb unit

The figure shows eight all-pass filters in series with each other, with taps after every two filters to feed the output network. The article also mentions that there should be approximately 200ms of delay in the feedback path (in addition to the all-pass filter delays) in order for the system to remain causal. (Because the all-pass filter output contains some component of the input that has not been delayed, feeding the all-pass filter output directly into its input requires prior knowledge of the output before it has been calculated. This can pose some problems.) Additionally, filters could be added into the loop to make the sound more natural.

I started by implementing something close to the Schroeder reverb unit in Simulink, but found that Simulink had a hard time running all of the parallel comb filters at the same time. Also due to the structure of the reverb unit, it was difficult to change the length of the reverb without making the order of the digital filters absurdly high.

Therefore, I made my design more closely match the reverb unit layout by Spin Semiconductors. This feedback allows for the length of the reverb to be easily adjusted by changing how much output signal is fed back into the unit. The figure below shows the impulse response of the system both with and without feedback. However, when feedback was added to the system, a noticeable high-pitched ringing was generated. In order to combat this, I added a low-pass FIR filter to the filter chain. This almost completely eliminated the ringing (it was also required to have a FIR filter for the project, so it was a double win).

Reverb impulse response with and without feedback
Reverb impulse response with, and without feedback

The FIR filter was designed to have a cutoff frequency of 3.5KHz. I used Matlab’s IIR filter generation tools to produce a frequency response, then used frequency sampling to find the FIR filter taps (A 64th order filter was used). The frequency sampling method takes the desired frequency response of the filter, then samples it evenly from 0 – 2/π radians/sample. The phase response of the filter can be set to anything; however, for the filter to be causal, the filter should have linear phase. This is ensured by making the phase of the filter φ(ω) = -jω (N – 1)/2 where N is the length of the filter.

The following figure shows the frequency, phase, and impulse response of the FIR filter. The code listed below was used to produce the filter plotted in the figure.

FIRFilter

The following code builds a FIR filter from a magnitude response

function [imp_resp, mag_resp, phase_resp] = fir_sellicott(mag_vec)
% fir_sellicott Create a FIR filter from the magnitude response vector
% the input is assumed to be equidistantly spaced between 0 and 2*pi

   N = length(mag_vec);
   resp = zeros(1, N+1);

   for k = 0:floor((N-1)/2)
       resp(k+1) = mag_vec(k+1)*exp(-1j*2*pi/N*k*(N-1)/2);
       resp(N-k+1) = conj(resp(k+1));
   end

   if mod(N, 2) == 0
       resp(N/2+1) = mag_vec(N/2+1);
   end

   mag_resp = abs(resp(1:N));
   phase_resp = angle(resp(1:N));

   imp_resp = ifft(resp, N);
end

The following function produces a magnitude vector from filter coefficients (numerator and denominator).

function H = buildtf(num_coeffs,den_coeffs,s)
  %build a transfer function from two vectors of coefficients
  num = 0; den = 0;
  num_ord = length(num_coeffs);
  den_ord = length(den_coeffs);

  for i = 1:num_ord
      num = num + num_coeffs(i)*s.^(num_ord - i);
  end

  for i = 1:den_ord;
      den = den + den_coeffs(i)*s.^(den_ord - i);
  end
  H = num./den;
end

Finally the following code snippet generates a filter using the two previous functions.

clear; close all; clc;
pkg load communications;

% produce FIR lowpass filter to get rid of HF in the reverb/just
% make it sound better.
Fs = 44100;  % audio sample rate
Ts = 1/Fs;

% filter settings
F_p = 3000;
f_p = F_p/(Fs/2);
N = 7;

% generate a digital filter in matlab
[b, a] = butter(N, f_p);

M = 31; % FIR filter order
N2 = 512; % number of points to make a detailed frequency response of the filter

% one omega for generation, another for plotting
omega = 2*pi*(0:1/M:(1-1/M));
omega2 = 2*pi* ( 0:1/N2:(1-1/N2) );

% one z for generation, another for plotting
z = exp(1j*omega);
z2 = exp(1j*omega2);

% one transfer function for getting sample points
% the other for plotting frequency response
H_sim = buildtf(b, a, z);
H_fin = buildtf(b, a, z2);

freq_hz = omega/(2*pi)*Fs;
freq_hz2 = omega2/(2*pi)*Fs;

% produce an FIR filter from the generated IIR filter
[h, H, phi] = fir_sellicott(abs(H_sim));

% get a detailed plot of the FIR frequency response
H_l = fft(h, N2);

% display one of the IIR all-pass filters I am using to produce the reverb effect
% these filters have the difference equation
% y[n] = a*x[n] + x[n-d] - a*y[n-d]
% H(z) = (a + z^(-d) ) / (1 + a z^(-d))
a = 0.5;
d = 1;
all_pass_H = (a + z2.^(-d) ) ./ (1 + a*z2.^(-d));

% plot the stuff
figure(1);

% Plot the magnitude of the various filters
subplot(3, 1, 1);
scatter(freq_hz, mag2db(H), 'b', 'filled'); hold on
plot(freq_hz2, mag2db(abs(H_fin)), 'r--', 'linewidth', 2);
plot(freq_hz2, mag2db(abs(H_l)), 'linewidth', 2, 'color', 'b'); hold off
xlabel('Frequency (Hz)');
ylabel('|H(f)|');
xlim([0, max(freq_hz2)]);
ylim([-80, 1]);
title('FIR filter 1: Magnitude');
legend('FIR sample points', 'IIR response', 'FIR response');

% Plot the phase of the filters
subplot(3, 1, 2);
scatter(freq_hz, phi, 'b','filled'); hold on
plot(freq_hz2, angle(H_l), 'linewidth', 2, 'color', 'b'); hold off
xlabel('Frequency (Hz)');
ylabel('angle(H)');
xlim([0, max(freq_hz2)]);
title('FIR filter 1: Phase');

% Plot the impulse response
subplot(3, 1, 3);
stem(0:length(h)-1, h, 'linewidth', 2, 'color', 'b', 'markerfacecolor', 'b');
xlabel('Time n (samples)');
ylabel('h[n]');
title('FIR Filter 1: Impulse Response');

figure(2);

% Plot the magnitude of the various filters
subplot(2, 1, 1);
plot(omega2, abs(all_pass_H), 'linewidth', 2, 'color', 'b'); hold off
xlabel('Frequency (Hz)');
ylabel('|H(f)|');
xlim([0, max(omega2)]);
title('All pass filter 1: Magnitude');

% Plot the phase of the filters
subplot(2, 1, 2);
plot(omega2, angle(all_pass_H), 'linewidth', 2, 'color', 'b'); hold off
xlabel('Frequency (Hz)');
ylabel('angle(H)');
xlim([0, max(omega2)]);
title('All pass filter: Phase');

% output taps
fprintf('{ ');
for i=1:(length(h)-1)
  fprintf('%f,', h(i));
end
fprintf('%f }\n', h(length(h)));

In order for the system to be causal, a delay needed to be added to the feedback path. The Spin Semiconductors page recommended a delay of at least 200ms in order to minimize “flutter” in the output signal. I used a delay of 8820 samples based on my sample rate of 44100Hz. Another recommendation was to add a taped delay filter to the reverb unit to smooth the outgoing signal (this is really just another FIR filter). This was added by breaking up the 8820 sample delay into three segments and use exponentially decaying coefficients at each break (1, 0.5, 0.25, 0.125). This made the reverb sound better to me, so I left it in the design. The following figure shows the final design in Simulink for the reverb unit.

ReverbUnit

C++

Now that the reverb unit was working in Simulink, the design needed to be ported to C++, to implement on the Raspberry Pi. I started by modifying some code I found on github which demonstrated mp3 decoding and playback using standard Linux libraries.

The project is split into three parts an audio decoding program, a playback program, and the actual filters. The first two parts are almost directly lifted from the previously mentioned code. The decoding program outputs the raw samples from the audio file over stdout where they are taken in by the filter program which filters them (well, duh). The filter program then outputs the filtered samples over stdout. Last, the playback program takes the samples and does magic on them so the Linux audio system can play it back.

One interesting thing to note is that audio files are typically stereo, where left and right samples are interleaved. This means that each buffer should be twice as long as designed, then every other index is used for each channel (e.g. all of the odd delays are left samples). This has to be taken into account when implementing the filters, so each buffer is twice as long as the filter design calls for, and the same action is performed on each audio channel.

While the code to handle the Linux audio interface was written in C, I wrote a few C++ classes to handle the filtering. The reverb unit class definition (named FilterProject) is listed below (and is avalible on github). This code shows the four all-pass filter objects (described later), the FIR filter object, and the buffer delay-line used by the reverb unit.

Reverb Unit Class

The following two code listings show the implementation of the FilterProject class and the initialization of all of the filter objects. The all-pass filters are initialized by passing in the number of samples they need to delay, and the feedback and feed-forward gains. The FIR filter is initialized by passing in an array of filter coefficients. These coefficients were simply copy-pasted from Matlab into the C++ project. The constructor for the reverb unit takes one parameter, the feedback gain of the reverb unit. This will set how long the reverb lasts.

Reverb class header

class FilterProject{
private:
  using outType = int16_t;
  using deque = std::deque;
  using buffPtr = std::unique_ptr;

public:
  //constructor
  FilterProject(float feedbackGain_ = 0.25);

  // the decoder and playback programs need the take the samples in the form of 8-bit integers
  uint8_t* get_samples(uint8_t* samples, size_t num_samples);

private:
  //define the names of the filters I will be using
  AllpassFilter allpass1;
  AllpassFilter allpass2;
  AllpassFilter allpass3;
  AllpassFilter allpass4;
  FIRFilter firFilter;
  buffPtr delay;
  outType do_filtering(outType new_x);

  float feedbackGain;
  const size_t sample_rate = 44100;
};

Reverb class construction

FilterProject::FilterProject(float feedbackGain_) :
  //initilize filters
  allpass1(4410, 0.6),
  allpass2(1470, -0.6),
  allpass3(490, 0.6),
  allpass4(163, -0.6),

  // pass in FIR coefficients to the FIR filter class
  firFilter({ 0.003369,0.002810,0.001758,0.000340,-0.001255,-0.002793,-0.004014,
  -0.004659,-0.004516,-0.003464,-0.001514,0.001148,0.004157,0.006986,0.009003,
  0.009571,0.008173,0.004560,-0.001120,-0.008222,-0.015581,-0.021579,-0.024323,
  -0.021933,-0.012904,0.003500,0.026890,0.055537,0.086377,0.115331,0.137960,
  0.150407,0.150407,0.137960,0.115331,0.086377,0.055537,0.026890,0.003500,
  -0.012904,-0.021933,-0.024323,-0.021579,-0.015581,-0.008222,-0.001120,
  0.004560,0.008173,0.009571,0.009003,0.006986,0.004157,0.001148,-0.001514,
  -0.003464,-0.004516,-0.004659,-0.004014,-0.002793,-0.001255,0.000340,
  0.001758,0.002810,0.003369 })
{
  delay = std::make_unique(3*2*2940, 0.0);

  feedbackGain = 0.25;
  if (feedbackGain_ > 0 && feedbackGain_ < 1) {
    feedbackGain = feedbackGain_;
  }
}

Now that the filters are initialized, they must be connected together to form reverb. The code listing below shows how the various filters are connected together in order to perform the reverb operation. First, the input sample is combined with the feedback signal, then the signal is passed through a series of filters (the four all-pass filters and the FIR filter). Note that the output of each filter is passed to the input of the next filter in the chain. The output from the last filter is placed into the delay line, where it goes through the final FIR filter.

auto &d = *delay.get();

// the coefficient on the d.back() sets how long the reverb
// will sustain: larger = longer
auto x = 0.7*new_x + feedbackGain*d.back();

// run through the all pass filters
// chain the outputs end to end
auto y1 = allpass1.do_filtering(x);
auto y2 = allpass2.do_filtering(y1);
auto y3 = allpass3.do_filtering(y2);
auto y4 = allpass4.do_filtering(y3);
auto y5 = firFilter.do_filtering(y4);

d.pop_back();
d.push_front(y5);

//add a bit of an FIR filter here, smooth the output
auto reverb_signal = y5 + 0.5*d[2*2940] + 0.25*d[2*2*2940] + 0.125*d[3*2*2940];
auto y = 0.6*reverb_signal + new_x;
return y;

For all of my filters a std::deque is used for the delay line. This container is used over a vector, as elements can be pushed and popped from both ends, not just the back. This allows me to turn the deque into a FIFO buffer quite easily. My use of the deque in this manner can be seen in line 15 and 16 from the prior code.

FIR Filter Class

The following listing shows the header file for the FIR filter class, the code for the all-pass filter is very similar. This code shows the data contained in the FIR filter, which is a vector that contains the filter coefficients, and a buffer which contains the audio samples being filtered. This code is listed here mostly for the sake of completeness, all of the work is done in the implementation file.

class FIRFilter {
using outType = int32_t;
using deque = std::deque;
using buffPtr = std::unique_ptr;

public:
  FIRFilter(std::vector taps_);
  ~FIRFilter() = default;

  outType do_filtering(outType new_x);

private:

  size_t delay;
  std::vector taps;
  buffPtr input_samples;
};

The next listing shows the implementation of the FIR filter class. It should be noted that the size of the delay buffer is set to be twice that of the order of the filter (line 5). This is due to the fact that stereo data is run through the filter, so at each time step there are actually two samples. Each time the do_filtering function is calls it runs on either the left or right channel of audio, so it has to be called twice in order to get one full sample. Also, in the constructor for the class, a vector of the coefficients is loaded into the internal class vector. These coefficients are taken from the impulse response of the filter, which are also the coefficients of the zeros of the filter’s transfer function.

FIRFilter::FIRFilter(std::vector taps_) :
  taps(taps_)
{
  // multiply by 2 to make room for the L/R audio channels
  delay = taps_.size()*2;

  // initilize buffer to zeros
  input_samples = std::make_unique(delay, 0.0);
}

FIRFilter::outType FIRFilter::do_filtering(outType new_x) {
  auto &amp;amp;x = *input_samples.get();

  // replace the oldest value for x
  x.pop_back();
  x.push_front(new_x);

  double new_val = 0;
  // implement the filter.
  for (auto i = 0; i < delay/2; ++i) {
    //multiply by 2 to account for LR channels
    new_val += taps[i]*x[2*i];
  }

  //return the newest value for y
  return new_val;
}

In the do_filtering function, each output sample is determined by the weighted sum of the delay buffer where the weights are the filter coefficients. This function simply implements the difference equation for an FIR filter. In general the equation is

y[n] = a_1*x[n] + a_2*x[n-1] + ... + a_n*x[n-M]

where the a_i’s are the filter coefficients and M is the order of the filter. Note that there is no feedback from previous output values in the filter, hence the impulse response of the filter is finite.

The equation is implemented by a for loop that sums the even valued delay elements (either all the left, or right audio samples) and weights them by their appropriate coefficient. This is a case where a speed improvement could probably be made by using a std::array over a std::vector, as this would allow the compiler to unroll the loop (as if i had hard coded all of the gains and delays in by hand); however I used the std::vector, as it allowed the filter to easily be generalized so that I could just copy in the filter coefficients from Matlab into the constructor.

All-Pass Filter Class

The all-pass filters are implemented in a similar way to the FIR filters. The same note about left and right audio channels applies in this filter as well. The following code listing shows the implementation of the filter in c++.

AllpassFilter::AllpassFilter(size_t delay_, float gain_) {

  // multiply by 2 to make room for the L/R audio channels
  delay = delay_ * 2;
  gain = gain_;
  // initilize the buffer to zeros
  delay_buff = std::make_unique<deque>(delay, 0.0);
}

AllpassFilter::outType AllpassFilter::do_filtering(outType new_x) {
  // grab the delay line
  auto &g = *delay_buff.get();

  //get the latest value to come through the delay line
  auto g_out = g.back();
  g.pop_back();

  // do the thing! follows the standard allpass filter structure
  auto g_in = new_x + gain*g_out;
  g.push_front(g_in);
  auto y = -gain*g_in + g_out;

  //return the newest value for y
  return y;
} 

From the figure showing the block diagram of an all-pass filter, it can be seen that the difference equations for the filter (in direct form 2) are:

g[n] = x[n] + a*g[n-M]

where g[n] is the input to the delay line and M is the number of delay elements in the filter. The output for the filter is:

y[n] = -a*g[n] + g[n-M]

Conclusions

The reverb unit turned out far better than I expected it to, and it forms a nice little demonstration project for what DSP can do. After completing the project, I think the best description of my feelings were best stated by He-Man. I may have gone a little reverb crazy, adding it to every sound file I had on my computer.

He-Man
…To make this sound effect

The next logical step for this project would be to implement it on a microcontroller so that audio could be passed through it to have reverb applied. But that is a project for another time.

Sam Ellicott
Soli Deo Gloria

Advertisements

2 thoughts on “C++ Reverb – DSP Final Project

  1. Hi Sam,

    Many congratulations, this was great to read! Have you ever come across the workflow that allows generating audio plugins directly from MATLAB code? If that sparks any interest, then take a look at the generateAudioPlugin function in the online MATLAB documentation. Since release R2019a that can also generate C++ JUCE projects. It feels like that would be another great way to port your orginal MATLAB design to a real-time implementation.
    Other thought: the MATLAB Plugin AES Student Competition from the Audio Engineering Society, which I am helping to organize, would be a perfect fit to this project in case you still qualify as student. More info at https://www.mathworks.com/academia/student-competitions/aes-matlab.html. The submission deadline is August 16.

    I would be happy to provide more details ditrectly if needed. Feel free to reach out.

    Thanks!
    Gabriele.

    Like

    1. Thanks for the information, I have not seen that before. It seems like whenever I am using MATLAB, I implement a thing and then I find a built-in function that does that thing better. 🙂

      Unfortunately, I will not have any extra time to work on a competition project as I am starting a new job this week. I wish I had the time, it looks like fun.

      Thanks,
      Sam Ellicott
      Soli Deo Gloria

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s