%
% Demo Example for a Anti-Vuvuzela-Filtering using
% discrete notch and lowpass filter algorithms
%
% R. Brigola, December 2014
%
% We use different implementations of a series of 2nd-order
% IIR-Notch-Filters and a discrete 3rd order Butterworth-Lowpass-Filter
% to damping Vuvuzela noise from a recording of a Soccer Game
% during the FIFA-World Cup 2010 in South Africa.
%
% The notches in the frequency characteristics are chosen at
% 233 Hz, 466 Hz (Octave), 699 Hz (Quint) and 932 Hz (Octave).
%
%
% More detailed:
%
% 1) WE WRITE FILTER-FUNCTIONS AS m-FILES, TEST THEM, OBSERVE THE
% COMPUTATION TIME AND THE PERSISTENT NOISE IN THE FILTERED DATA
%
% At first, we construct by the bilinear transform - based on an
% analog frequency characteristics - a discrete 2nd-order Notch-Filter
% for a single notch frequency f as a Matlab function with usage
% [a,b]=myNotch(f,Q,dc_Gain,Fs).
% Q is a filter design parameter, dc_gain the gain of the filter.
% Fs is the sampling frequency used in the application audio file.
% a=[ao,a1,a2], b=[1,b1,b2] are the filter coefficients of the rational
% transfer function H(z) with polynomials in 1/z and coefficients
% a in the nominator, b in the denominator ordered as in the transfer
% function H(z)=(a0+a1/z+a2/z^2)/(1+b1/z+b2/z^2). That coefficients are
% computed in the m-file mynotch.m
%
% Then we design analogously a discrete 3rd order Butterworth Lowpass
% whose coefficients are again saved in a=[a0,a1,2,a3] and b=[b0,b1,b2,b3].
% Usage: [a,b]=myButter(f,dc_Gain,Fs) with f the cutoff frequency,
% dc_Gain, Fs as above.
%
% 2) WE REPEAT THE FILTERING WITH THE COMPUTED FILTER COEFFICIENTS, BUT
% USE MATLAB'S INTERNAL FILTER FUNCTION "filter", OBSERVING AGAIN THE
% COMPUTATION TIME
%
% 3) WE SUBSTITUE MATLAB'S FILTER FUNCTION BY A SELF-IMPLEMENTED FILTER
% FUNCTION "myfilter" - WRITTEN IN C++ AND COMPILED WITH MEX - AND
% COMPARE AGAIN THE COMPUTATION TIMES.
% Remark:
% In case your mex doen't work with Microsoft Visual C++ 2010 SP1
% install the patch as described at the URL
% http://www.microsoft.com/en-us/download/details.aspx?id=4422
%
%
% Input:
% Audio filenames
% Samples per second
% Number of notch frequencies and their location
% Filter quality factor Q, dc_Gain A0
% Switch for lowpass on/off
% Gain and cutoff frequency for the lowpass
% the notch frequencies
% Notch and Lowpass Filter Order
% Frequency range for the plots of the absolut
% amplitude distortions of the filters
clear all;
tic; % start_time=tic %cputime; % to compute computation time
% for a few processing steps
%%%%%%%%%%%%% Set the input data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
audio_input = 'Vuvuzela.wav'; % Test also with the file Vuvuzela.wav
audio_output= 'Vuvuzela_filtered.wav';
Fs=44100;
% We assume, that the audio files are stereo recordings with 2 channels
NF=4; % Number of notch frequencies
cutoff=2000; % lowpass cutoff frequency in Hz
switchlow=1; % switch=1: lowpass on, switch=0 lowpass off
Lowpass_Gain=2; % lowpass dc-Gain >1 to get back signal energy lost
% by lowpass filtering
f(1:NF)=0; % the vector for the NF notch frequencies
f(1)=233; % first notch frequency
f(2)=466; % the next harmonics, here Octave
f(3)=699; % Quint
f(4)=932; % 2nd Octave
% f(5)=1396; % 2nd Quint, if wanted
% f(6)=1864; % 3rd Octave, if wanted
NotchOrder=2; % Order of the Notch Filter, i.e. we use simple Biquads
LowpassOrder=3;% Order of the Bitterworth Lowpass Filter
Q=5; % filter parameter that controls the decay at the notch
% often called quality factor Q
A0=1; % dc-Gain for all consecutive notch filters
P=6000; % frequency range for the plots, here 6000 Hz
%%%%%%%%%%%%%%% end of input data section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Initialize Data Matrices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(switchlow==1)
NFF=NF+1; % dimension for signal matrix
% according to number of applied filters
else NFF=NF;
end;
% Now read the audio wav-file that should be processed
% It is assumed that we have a stereo file with 2 channels
[x,Fs]=audioread(audio_input);
N=size(x,1); % No. of samples per channel
x(1:5,1:2)=0; % set a few samples to zero as initial condition
nom=zeros(NF,NotchOrder+1); % matrices for the coefficients
denom=zeros(NF,NotchOrder+1); % of NF notch filters
preparation_time=toc % processing computation time
tic; %start_time=cputime;
%%%%%%%%%%%%% Compute the Notch Filter Coefficients %%%%%%%%%%%%%%%%%%%%%%
%%%%% Call myNotch for the filter coefficients of the notch filters %%%%%
for k=1:NF; % here, myNotch computes the coeff. of 2nd. order IIR-Filters
[nom(k,1:NotchOrder+1),denom(k,1:NotchOrder+1)]=myNotch(f(k),Q,A0,Fs);
end;
notch_computation_time=toc %cputime - start_time;
% initialize matrices for the successively filtered audio data
Y=zeros(N,2);
Z=zeros(N,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 1) Successive Notch Filtering with the m-file myBiquad.m %%%%%%%%%%
%%%%%%%%%%%%%%%% Start Filtering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%% 2) Substitute below 'myBiquad' by 'filter', i.e. use Matlab's %%%
%%%% internal filter routine and compare the computation time %%%
%%%%
%%%% 3) Substitute the 'filter' routine now by the routine 'myfilter' %%%
%%%% compiled with mex -setup C++ option from 'myfilter.c' %%%
%%%% and compare again the computation time. The function 'myfilter'%%
%%%% implements a so-called 'direct form II transposed' filter %%
%%%% with a rational transfer function %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic; % start filtering the audio file
% Apply the Notch Filters with the filter function myBiquad
Z=myfilter(nom(1,:),denom(1,:),x); % filter(nom(1,:),denom(1,:),x);
% myBiquad(nom(1,:),denom(1,:),x);
for k=2:NF;
Y=myfilter(nom(k,:),denom(k,:),Z); % filter(nom(1,:),denom(1,:),x);
% myBiquad(nom(k,:),denom(k,:),Z);
Z=Y;
end;
notch_filtering_time=toc % processing computation time
tic; % start time lowpass filtering
% Calculate the 3rd order Butterworth lowpass
% then apply the lowpass, here with fixed order 3
if(switchlow==1)
[nomlow(1:4),denomlow(1:4)]=myButter(cutoff,Lowpass_Gain,Fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lowpass filtering both stereo channel data, if filter switch is set
% 1) Use the filter function 'my3rdOrderIIR' from the corresponding m-file
% 2) Substitute it by Matlab's 'filter' as above
% 3) Substitute it by the function 'myfilter'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y=myfilter(nomlow,denomlow,Z); % my3rdOrderIIR(nomlow,denomlow,Z);
% filter(nomlow,denomlow,Z); myfilter(nomlow,denomlow,Z);
end;
lowpass_filtering_time=toc % processing computation time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Now write the filtered audio data and plot the amplitude
%%%%%% distortions of the filters.
%%%%%% (In that demonstration we do not further work as phase corrections
%%%%%% after the IIR-filtering).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic; % start write process
% Write the Output Audio wav-File
audiowrite(audio_output,Y,Fs,'BitsPerSample',24);
writing_time=toc % time for the write process
% Listen to short sections of the original and of the filtered sound
fprintf('\nPlay 15 seconds of the Original Sound:');
sound(x(1:15*Fs,:),Fs,24);
pause(17);% pause while sound is playing
fprintf(' OK\n');
fprintf('Play 15 seconds of the filtered Sound: ');
sound(Y(1:15*Fs,:),Fs,24);
pause(17);% pause while sound is playing
fprintf('OK\n');
% Plot the Amplitude Distortions of the Filters:
% The Notch filter
tic; % start time neede for the plots
P=P+1;
w(1:NF)=0;
h(1:P)=0;
for m=1:P;
for k=1:NF;
e_vec1=exp_vector(NotchOrder,m-1,Fs);
% that vector contains the exponential terms
% used in the nominator and the denominator
% of the frequency characteristics
w(k)=(nom(k,:)*e_vec1')/(denom(k,:)*e_vec1');
end;
h(m)=prod(w,2); % saves the amplitude distortions of the
% notch filtering for all given frequencies
end;
f=linspace(0,P-1,P);
figure(1);plot(f,abs(h(1:P)));
title('Single-sided Amplitude Distortion of the Notch Filter');
xlabel('Frequency in Hz');
ylabel('|Amplitude|');
% Butterworth Lowpass
if(switchlow==1) % i.e. we used the lowpass
hh(1:P)=0;
for m=1:P;
e_vec2=exp_vector(LowpassOrder,m-1,Fs);
w=(nomlow(1:4)*e_vec2')/(denomlow(1:4)*e_vec2');
hh(m)=w; % amplitude distortion of the lowpass
h(m)=w*h(m); % and of the notch plus lowpass series
end;
figure(2);plot(f,abs(hh(1:P)));% plot the amplitude distortions
% of the lowpass
title('Amplitude Distortion of the Lowpass Filter');
xlabel('Frequency in Hz');
ylabel('|Amplitude|');
% plot analogously the amplitude distortions of notch and lowpass in series
figure(3);plot(f,abs(h(1:P)));
title('Amplitude Distortion of the Notch and Lowpass Filtering');
xlabel('Frequency in Hz');
ylabel('|Amplitude|');
end;
plotting_time=toc % computation time for the plots