#   prog1p1.mpl -- Restoration of Analog Audio Recordings
#   Written November 4, 2011 by Eric Olson for Math 761
#
#   This program was developed and tested using Maple 9.5
#   and should also work with other versions.

restart:
kernelopts(printbytes=false):
with(LinearAlgebra):
with(DiscreteTransforms):
Digits:=15:
UseHardwareFloats:=true:
kfunc:=t->alpha/(sigma*sqrt(2*Pi))*exp(-t^2/(2*sigma^2))
      +beta/(sigma*sqrt(4*Pi))*exp(-(t+0.5)^2/(4*sigma^2)):

printf("prog1p1.mpl -- Restoration of Analog Audio Recordings\n"
       "Written November 4, 2011 by Eric Olson for Math 761\n\n");

alpha:=1.75:
beta:=0.5:
sigma:=0.0004:
epsilon:=0.000001:

printf("Reading signalw.dat...\n");
fd:=fopen("signalw.dat",READ):
signalw:=readdata(fd,float):
fclose(fd):
w:=Vector(signalw,datatype=complex[8]):
N:=nops(signalw):
sN:=sqrt(N):
k:=Vector([seq(evalhf(sN*kfunc(j/44100.0)/44100.0),j=0..N/2-1),
     seq(evalhf(sN*kfunc(j/44100.0)/44100.0),j=-N/2..-1)],
     datatype=complex[8]):

printf("Approximating X...\n");
w:=FourierTransform(w):
k:=FourierTransform(k):
for i from 1 to N
do
    if evalf(abs(k[i]))>epsilon
    then
        w[i]:=w[i]/k[i];
    else
        w[i]:=0.0
    end
end:
w:=InverseFourierTransform(w):

printf("Writing wtox.dat...\n");
buf:=[seq(Re(w[j]),j=1..N)]:
fd:=fopen("wtox.dat",WRITE):
writedata(fd,buf,float):
fclose(fd):

printf("Reading signalx.dat...\n");
fd:=fopen("signalx.dat",READ):
signalx:=readdata(fd,float):
fclose(fd):
x:=Vector(signalx,datatype=complex[8]):

printf("\nThe error is E=%g\n",
	VectorNorm(w-x,2)/VectorNorm(x,2));