%  This script solves the example in lecture notes where a sample voice
%  signal is given (a repeated vowel sound) and the average specturm
%  of this sound is used to design a filter to separate the voice signal
%  from white noise.
% 
%    Written by Kevin D. Donohue (donohue@engr.uky.edu) March 2004 
%

%  Filter parameters
ford = 500;  %  FIR filter order
fftwin = 256;  %  FFT window length in samples
nfft = 2*fftwin;  %  Number of FFT points
psdolap = fix(fftwin/2);  %  numver of overlapping samples in window

%  Open sample voiced signal and show it average spectrum
[y,fs] = wavread('myvowel.wav');  %  Read in wavefile
[pd,f] = pwelch(y,hamming(fftwin),psdolap,nfft,fs);  %  Examine average specturm
figure(1); semilogy(f,(fs*pd/2))
xlabel('Hz')
ylabel('Power')
title('Used this figure to determine shape of spectral filter')
pause
%  Create an outline of the spectrum with frequency and amplitude points
fa = [0 0; 65 0.05; 100  0.11; 150  0.07; 230  0.22; ...
      300  0.1; 375  0.36; 600  0.01; 820  0.06; 1100  3.506e-3; ...
      1280  .02; 1500 .005; 1850 1.5e-3; 2400 .01; 3000 3.074e-3; 3600  0.0005; ...
      4000 4e-4; fs/2  1e-6]
hold on
% Compare to average spectrum
semilogy(fa(:,1),fa(:,2),'r')
hold off
title('Compare filter specification points (red) to data spectrum')

f1 = 2*fa(:,1)/fs;  %  Normalize frequencies wrt Nyquist
a = sqrt(fa(:,2));  %  Extract amplitudes convert from power
b = fir2(ford,f1,a); %  Obtain FIR filter

%  Examine filter spectrum
figure(2)
[h,ff] = freqz(b,1,2048,fs);
plot(ff,abs(h))
xlabel('Hz')
ylabel('Magnitude')
title('resulting filter spectrum')

%  Get new voice signal
[yn, fs2] = wavread('mytalk.wav');
if fs2 ~= fs  %  Check sampling rate
    disp('Problem!  the sampling frequenies are not equal')
else
%  Add white noise 
nos = randn(size(yn));
%  Find RMS power in mytalk and noise
rmstalk = std(yn);
rmsnos = std(nos);
%  Computed 6 dB ratio
sigovernoise = 10^(6/20);  %  Convert dB to gain of signal avbove noise
scnos = rmstalk/sigovernoise;  %  Compute scaling factor for normalized rms noise
signos = yn + scnos*nos/rmsnos;
ftalk = filter(b,1,signos);  %  Filter new voice signal
soundsc(signos,fs)           %  play old sound
pause(length(signos)/fs)     %  Let it finish playing
soundsc(ftalk,fs)            %  Play filtered sound
end