%  This script solves the example in lecture notes where a sample voice
%  signal is given and the average specturm 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 = 512;  %  FFT window length in samples
nfft = 4*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('nad1.wav');  %  Read in wavefile
[pd,f] = pwelch(y,hamming(fftwin),psdolap,nfft,fs);  %  Examine average specturm
figure(1); plot(f,10*log10((fs*pd/2)))
xlabel('Hz')
ylabel('Magnitude dB')
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; 118.4 0.03343;  236.9  0.05724; ...
      484.5  0.03562; 979.8  3.593e-4; 1615 2.211e-3; 2832  3.389e-5; ...
      3478 3.566e-4; 4425  9.793e-6; 5103 4.475e-5; 5501 5.043e-6; ...
      6761 4.283e-5; 7655 7.787e-6; 8043 2.635e-5; fs/2  1e-6]
hold on
% Compare to average spectrum
plot(fa(:,1),10*log10(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 into own vector (convert power to amplitude)
b = fir2(ford,f1,a); %  Obtain FIR filter

%  Examine filter spectrum
figure(2)
[h,ff] = freqz(b,1,2048,fs);
plot(f,10*log10((fs*pd/2))); hold on; plot(ff,20*log10(abs(h)),'g'); hold off
xlabel('Hz')
ylabel('dB')
title('resulting filter spectrum (green)')

%  Get new voice signal
[yn, fs2] = wavread('nad3.wav');
if fs2 ~= fs  %  Check sampling rate
    disp('Problem!  the sampling frequenies are not equal')
else
[pdn,f] = pwelch(yn,hamming(fftwin),psdolap,nfft,fs);  %  Examine average specturm
figure(3); plot(f,10*log10((fs*pdn/2))); hold on; plot(ff,20*log10(abs(h)),'r'); hold off
xlabel('Hz')
    %  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