%   prog1p1.m -- Restoration of Analog Audio Recordings
%   Written October 29, 2011 by Eric Olson for Math 761
%
%   This program was developed and tested using GNU Octave.
%   It should also work with MATLAB.

clear all
more off

N=262144;
alphaz=1.75;
betaz=0.5;
sigmaz=0.0004;
epsilonz=0.000001;
kfunc=@(t) alphaz*exp(-t.^2/(2*sigmaz^2))/(sigmaz*sqrt(2*pi)) ...
    +betaz*exp(-(t+0.5).^2/(4*sigmaz^2))/(sigmaz*sqrt(4*pi));

disp('prog1p1.m -- Restoration of Analog Audio Recordings');
disp('Written October 29, 2011 by Eric Olson for Math 761');
disp('');

disp('Reading signalw.dat...');
load signalw.dat;

w=signalw;
k=(kfunc([0:N/2-1,-N/2:-1]/44100.0)/44100.0)';

disp('Approximating X...');
k=fft(k);
w=fft(w);
for i=1:N
    if abs(k(i))>epsilonz
        w(i)=w(i)/k(i);
    else
        w(i)=0;
    end
end
w=ifft(w);

disp('Writing wtox.dat...');
buf=real(w);
save -ascii wtox.dat buf

disp('Reading signalx.dat...');
load signalx.dat;
x=signalx;

disp(sprintf('\nThe error is E=%g',norm(w-x)/norm(x)));
