(*  prog1p1.c -- Restoration of Analog Audio Recordings
    Written November 16, 2011 by Eric Olson for Math 761

    This program makes use of the fast Fourier transform given
    in the file fourier.pas.  The commands

        fpc fourier.pas
        fpc prog1p1.pas

    will compile this program using the FPK Pascal compiler in
    most POSIX/GNU environments.
*)

Program Main;

Uses ucomplex,fourier;

const N=262144;
    alpha=1.75;
    beta=0.5;
    sigma=0.0004;
    epsilon=0.000001;

type nvect=array [1..N] of complex;

(* The filter kfunc specifies the loss of high frequencies and tape
   print through by means of the constants alpha, beta and sigma. *)

function kfunc(t:double):double;
var tm1,tm2:double;
begin
    tm1:=t;
    tm2:=t+0.5;
    kfunc:=alpha*exp(-tm1*tm1/(2*sigma*sigma))/(sigma*sqrt(2*PI))
        +beta*exp(-tm2*tm2/(4*sigma*sigma))/(sigma*sqrt(4*PI));
end;

procedure readsignal(var X:nvect;p:string);
var buf: array [1..N] of single;
    r:longint;
    fd:file;
    ii:longint;
begin
    writeln('Reading ',p,'...');
    {$I-}
    assign(fd,p);
    reset(fd,sizeof(single));
    {$I+}
    if(ioresult<>0) then begin
        writeln(stderr,'Can''t open ',p,' for read!');
        halt(1);
    end;
    blockread(fd,buf,N,r);
    if(r<>N) then begin
        writeln(stderr,'Error reading data from ',p,'!');
        halt(2);
    end;
    close(fd);
    for ii:=1 to N do X[ii]:=buf[ii];
end;

procedure writesignal(var X:nvect;p:string);
var buf: array [1..N] of single;
    r:longint;
    fd:file;
    ii:longint;
begin
    writeln('Writing ',p,'...');
    for ii:=1 to N do buf[ii]:=X[ii].re;
    {$I-}
    assign(fd,p);
    rewrite(fd,sizeof(single));
    {$I+}
    if(ioresult<>0) then begin
        writeln(stderr,'Can''t open ',p,' for write!');
        halt(3);
    end;
    blockwrite(fd,buf,N,r);
    if(r<>N) then begin
        writeln(stderr,'Error writing data to ',p,'!');
        halt(4);
    end;
    close(fd);
end;

(*  Compute relative error how well X approximates x *)

function relerr(var XX,x:nvect):double;
var r1,r2:double;
    ii:longint;
    t:complex;
begin
    r1:=0.0;
    r2:=0.0;
    for ii:=1 to N do begin
        t:=x[ii];
        r1+=t.re*t.re+t.im*t.im;
        t-=XX[ii];
        r2+=t.re*t.re+t.im*t.im;
    end;
    relerr:=sqrt(r2/r1);
end;

var K,W,X:nvect;
    ki:complex;
    t:double;
    ii:longint;
begin
    writeln('prog1p1.pas -- Restoration of Analog Audio Recordings');
    writeln('Written November 16, 2011 by Eric Olson for Math 761');
    writeln;

    readsignal(W,'signalw.raw');

    (*  The function kfilt(t) takes an argument of time in seconds
        whereas the K[i] is sampled at 44.1 khz.

            i=2 corresponds to t=1/44100 of a second and
            i=N corresponds to t=-1/44100 of a second. *)

    for ii:=1 to N do begin
        if(ii<=N div 2) then t:=(ii-1)/44100.0
        else t:=(ii-N-1)/44100.0;
        K[ii]:=kfunc(t)/44100.0;
    end;

    writeln('Approximating X...');

    fft(K,N,1);
    fft(W,N,1);
    for ii:=1 to N do begin

        (*  This is the line that needs to be modified to implement
            the optimal Wiener-Kolmogorov filter. *)

        Ki:=K[ii];
        if(cmod(Ki)>epsilon) then W[ii]/=Ki
        else W[ii]:=0;
    end;
    fft(W,N,-1);

    writesignal(W,'wtox.raw');

    (*  Compute the error in the recovered signal *)
        
    readsignal(X,'signalx.raw');

    writeln;
    writeln('The error is E=',relerr(W,X));
    halt(0);
end.