% code modified from: http://www.aquaphoenix.com/lectures/matlab10/page4.html
% Workshop  	UN/ESA/NASA/UAE Workshop, Lecture at UAE University, Al Ain, UAE, November 21-22 2005

%It's important to capture the sampling frequency at which the sound was recorded, otherwise the speed of playback and results of further processing is not guaranteed to be correct:

[funky, f] = wavread('funky.wav');

%To play a wave file at sampling frequency f:

soundsc(funky, f);

%To view the waveform, plot the wave. Since audio is represented with many thousand samples per second, it may be required to plot small portions of the waveform at a time.

subplot(2,1,1), plot(funky), title('Entire waveform');
smallRange = 100000:100000+floor(f/100);
subplot(2,1,2), plot(smallRange, funky(smallRange)), title('100 milliseconds');

% Filtering

% We will examine audio filtering in the sense of specific frequency suppression and extraction. There are many different types of filters available for the construction of filters. We will specifically use the Butterworth filter.

%Matlab includes function butter for building Butterworth filters of three sorts:
%
%    * 'low' : Low-pass filters, which remove frequencies greater than some specified value.
%    * 'high' : High-pass filters, which remove frequencies lower than some specified value.
%    * 'stop' : Stop-band filters, which remove frequencies in a given range of values. 

%Frequencies values are specified in normalized terms between 0.0 and 1.0, where 1.0 corresponds to half the sampling frequency: f/2. A given frequency is thus expressed in terms of this value, for example, 1000Hz = 1000/(f/2).

%Filters are described in terms of 2 vectors ([b, a] = [numerator, denominator]).

%To apply a filter to a 1-D audio waveform, Matlab provides function filtfilt , which takes as arguments the result [b, a] from butter, the waveform, and a value denoting the order (number of coefficients) of the filter.

%A filter's frequency response can be plotted using function freqz . Magnitude values at zero dB are unaffected by the filter. Magnitude values below 0 dB are suppressed.

% Low-pass filter
% We design a 10th order low-pass filter to supress frequencies higher than 200Hz.

clf;
subplot(4,1,1); plot(funky(1:1000));

fNorm = 200 / (f/2);
[b,a] = butter(10, fNorm, 'low');
funkyLow = filtfilt(b, a, funky);

subplot(4,1,2); plot(funkyLow(1:1000));

soundsc(funkyLow,f);

% High-pass filter
% We design a 10th order high-pass filter to supress frequencies below 5kHz.
fNorm = 5000 / (f/2);
[b, a] = butter(10, fNorm, 'high');
funkyHigh = filtfilt(b, a, funky);

subplot(4,1,3); plot(funkyHigh(1:1000));


soundsc(funkyHigh,f);

% Stop-band filter - allows only frequencies between 500 Hz and 2500 Hz
fNorm = [500/(f/2), 2500/(f/2)];
[b, a] = butter(10, fNorm, 'stop');
funkyBand = filtfilt(b, a, funky);

subplot(4,1,4); plot(funkyBand(1:1000));

soundsc(funkyBand,f);