% Script for Wiener filter
load handel   % y is the signal stored in handel.mat
soundsample = y(1:2^16);   % to have the same as in dftfilter.m
n = length(soundsample);   % n is the length of signal
p = 1.3;                              % parameter for cutoff
noise = std(soundsample)*.50;                   % 50 percent noise level
noisecorrupt = noise*randn(n,1);                 % pure noise to be added to signal
noisysample = soundsample + noisecorrupt;                              % noisy signal
FTsample = fft(noisysample);
modFTsample = conj(FTsample).*FTsample;         % take fft of signal, and
approxmodnoise = max(modFTsample-n*(p*noise)^2,0); % apply cutoff
Filtermatrix = approxmodnoise./modFTsample;                 % define phi as derived
filteredsample = real(ifft(Filtermatrix.*FTsample));           % invert the fft
RMSerr_noisy = norm(noisysample-soundsample)/sqrt(n)
RMSerr_filtered = norm(filteredsample-soundsample)/sqrt(n)
disp('Hit key to play original sound')
pause
sound(soundsample)
disp('Hit key to play noisy sound')
pause
sound(noisysample)
disp('Hit key to play filtered sound')
pause
sound(filteredsample)
