Today I’d like to introduce a guest blogger, Charu, who is a Technical Marketing Manager at The MathWorks. Charu regularly presents webinars on signal processing topics, and most recently presented a webinar highlighting R2010a capabilities in the area of MATLAB signal processing.
- Set up input and output streams
- Extracting the noise spectrum
- Identifying vuvuzela frequencies
- Designing the parametric equalizers
- Visualizing magnitude response of the filters
- Implementing the parametric equalizers
- Filtering the input in a stream processing loop
- Visualizing and listening to output results
- Things to Try
Hi there - soccer fever is in full swing, and the incessant sounds of the vuvuzela are a big part of the soccer viewing experience. This South African instrument has become so popular that it now has a dedicated Facebook page! Well, 2 weeks ago, Loren featured a very interesting MATLAB Central submission by Vincent Choqueuse on vuvuzela denoising using the spectral subtraction method. This cool demo inspired me to try out other filtering techniques for the vuvuzela denoising problem.
A quick search on “vuvuzela filtering” suggested that notch filters and adaptive filters were the most common methods, so I decided to try a method that uses notch parametric equalizers. Parametric equalizers allow the user to control the bandwidth, gain and center frequency of the filter. Since the game commentary is a streaming broadcast, I also wanted to see if I could apply stream processing techniques using new System objects in MATLAB.
This example demonstrates how you can design parametric equalizer filters in MATLAB and use them to filter out vuvuzela sounds from an audio stream. This demo uses Signal Processing Toolbox, Filter Design Toolbox and Signal Processing Blockset.
First, let us look at how we can create streaming inputs and outputs in MATLAB. To get the data, download 'vuvuzela.wav' from this MATLAB Central entry.
We need an input audio stream that can read data from a WAV file. To do this, we can use the MultimediaFileReader System object.
Create a System object to read an audio stream from the WAV file
hSound = signalblks.MultimediaFileReader(... 'Vuvuzela.wav',... 'SamplesPerAudioFrame', 1024, ... 'AudioOutputDataType', 'double');
Get the input sampling frequency
hInfo = info(hSound); Fs = hInfo.AudioSampleRate;
Create a System object to play back audio to the sound card
hPlayer = signalblks.AudioPlayer('SampleRate', Fs);
Create a System object to write the demo output to a WAV file.
hOutputFile = signalblks.MultimediaFileWriter(... 'with_then_without_vuvuzela.wav',... 'SampleRate', Fs);
Create System objects for logging pure noise and filtered output signals.
hLogNoise = signalblks.SignalLogger; hLogOutput = signalblks.SignalLogger;
The first 5 seconds of the input do not contain any significant data, ie, no noise even, so let’s advance the input by 5 seconds, i.e., about 100 frames (@Fs = 22050Hz, 1024 samples per frame).
for frame_index=1:100 step(hSound); end
To obtain our “pure noise” signal, we use a 1-second portion of the signal that has only the vuvuzela sound, but no speech.
for noise_frames = 1:20 noise = step(hSound); step(hLogNoise, noise); end
Plot the noise spectrum.
pwelch(hLogNoise.Buffer, , , 8192, Fs);
Let’s zoom into the plot to identify noise peaks.
The fundamental and the first 4 harmonic frequencies of the vuvuzela are at: 235Hz, 476Hz, 735Hz, 940Hz, 1180Hz.
Our next task is to design notch parametric equalizers that will selectively filter out these frequencies. Since the notch filters may also remove components of the speech signal at the above frequencies, we will apply a peak filter to boost the speech signal at the output of the notch filters.
To design parametric equalizers, we need to specify filter order, center frequency, and the bandwidth or Q-factor of the filter. To determine the bandwidth or Q-factor of the filter, let’s take a closer look at the noise spectrum. At the fundamental and 1st harmonic frequeuncies (235Hz, 465Hz), we see that the peaks are sharp. Beyond the 1st harmonic, the noise peaks are spread across a wider band. So, for the notch filters at the lower frequencies, we specify a narrow stopband, i.e., a high Q, and for higher harmonics, we specify a wider stopband, i.e., a low Q. After some trial and error, I arrived at the following parameters:
Here are the specifications for notch parametric equalizers.
F01 = 235; F02 = 476; F03 = 735; F04 = 940; F05 = 1180; N = 2; %Filter order Gref = 0; %Reference gain G0 = -20; %Attenuate by 20dB Qa = 10; %Higher Q-factor for fundamental and 1st harmonic Qb = 5; %Lower Q-factor for higher harmonics f1 = fdesign.parameq('N,F0,Qa,Gref,G0',N,F01,Qa,Gref,G0,Fs); f2 = fdesign.parameq('N,F0,Qa,Gref,G0',N,F02,Qa,Gref,G0,Fs); f3 = fdesign.parameq('N,F0,Qa,Gref,G0',N,F03,Qb,Gref,G0,Fs); f4 = fdesign.parameq('N,F0,Qa,Gref,G0',N,F04,Qb,Gref,G0,Fs); f5 = fdesign.parameq('N,F0,Qa,Gref,G0',N,F05,Qb,Gref,G0,Fs);
Here are specifications for peak parametric equalizer.
F06 = 480; %Center frequency to boost N_peak = 8; %Use higher filter order for peak filter Qc = 0.5; %Use very low Q-factor for boost equalizer filter G1 = 9; %Boost gain by 9dB f6 = fdesign.parameq('N,F0,Qa,Gref,G0',N_peak,F06,Qc,Gref,G1,Fs);
Hp1 = design(f1, 'butter'); Hp2 = design(f2, 'butter'); Hp3 = design(f3, 'butter'); Hp4 = design(f4, 'butter'); Hp5 = design(f5, 'butter'); Hp6 = design(f6, 'butter');
hFV = fvtool([Hp1 Hp2 Hp3 Hp4 Hp5 Hp6]); legend(hFV,'NotchEQ #1', 'NotchEQ #2', 'NotchEQ #3',... 'NotchEQ #4','NotchEQ #5','Peak EQ');
To implement each filter, we use a System object that realizes an SOS (second-order-sections) structure. Using a System object to implement the filter enables us to apply the filter in a stream processing loop, without having to worry about updating filter states between each iteration.
hEQ1 = signalblks.BiquadFilter('SOSMatrix',Hp1.sosMatrix,'ScaleValues',Hp1.ScaleValues); hEQ2 = signalblks.BiquadFilter('SOSMatrix',Hp2.sosMatrix,'ScaleValues',Hp2.ScaleValues); hEQ3 = signalblks.BiquadFilter('SOSMatrix',Hp3.sosMatrix,'ScaleValues',Hp3.ScaleValues); hEQ4 = signalblks.BiquadFilter('SOSMatrix',Hp4.sosMatrix,'ScaleValues',Hp4.ScaleValues); hEQ5 = signalblks.BiquadFilter('SOSMatrix',Hp5.sosMatrix,'ScaleValues',Hp5.ScaleValues); hEQ6 = signalblks.BiquadFilter('SOSMatrix',Hp6.sosMatrix,'ScaleValues',Hp6.ScaleValues);
To filter out the vuvuzela noise, we use a stream processing loop that reads in one frame of input data from the WAV file, applies the series of 6 parametric equalizers, and then plays back the filtered output through the sound card.
for frames = 1:400 %Step through input WAV file one frame at a time input = step(hSound); %Apply notch EQ filters first, then peak EQ filter out1 = step(hEQ1, input); out2 = step(hEQ2, out1); out3 = step(hEQ3, out2); out4 = step(hEQ4, out3); out5 = step(hEQ5, out4); filt_output = step(hEQ6, out5); %Play 200 frames of input signal, then 200 frames of filtered output %to compare original and filtered signals if frames < 200 step(hPlayer, input); else step(hPlayer, filt_output); end %Log filtered output to buffer, WAV file step(hLogOutput, filt_output); step(hOutputFile, filt_output); end
Listen to the noisy and filtered signals here.
Visualize the spectrum of the filtered signal. If you zoom into the plot, you can see that the vuvuzela frequencies have been attenuated by our notch parametric equalizers.
pwelch(hLogOutput.Buffer, , , 8192, Fs);
Close input and output System objects to ensure that we leave all objects in their proper states.
close(hSound); close(hPlayer); close(hOutputFile);
This demo was a lot of fun to work on and I hope you’ll find it interesting. Feel free to download the code and try your own equalizer parameters. Have fun!