out ccon;
%off nat;
on rounded;

procedure delta(ii,jj);
begin 
if ii eq jj then return 1 else return 0; 
end;

n:=15;m:=6;
operator a,k,r,f;

for all ii such that ii<0 or ii>m-1 let k(ii)=0;
for all ii such that ii<0 or ii>n-1 let a(ii)=0;
for all ii such that ii<0 or ii>n-1 let r(ii)=0;
for all ii,jj let f(ii)*f(jj)=f(ii+jj);
for all ii let f(ii)*f(ii)=f(ii+ii);
f(0):=1;

for ii:=-10:20 do write(k(ii));
for ii:=-10:20 do write(a(ii));
for ii:=-10:20 do write(r(ii));

%-------------A-----------
for ii:=0:(n-1) do a(ii):=0; a(7):=100;
%-------------centered psf kernel-----------
for ii:=0:(m-1) do k(ii):=0; k(0):=0.5; k(1):=0.3; k(2):=0.2; 

%-------------convolution procedure-----------
%psf kernel dim = m
%------------------
procedure cconv();
begin

for tt:=0:(n-1) do <<

if tt > (m - 1) then write( r(tt) := for ii:=0:(m-1) sum a(tt-ii)*k(ii) )
		else write(r(tt):=0);
>>;
end;
%-----------------------------------
%-------------convolution procedure-----------
%psf kernel dim = n
%------------------
procedure cconvn();
begin

for tt:=0:(n-1) do <<

if tt > (m - 1) then write( r(tt) := for ii:=0:(n-1) sum a(tt-ii)*k(ii) )
		else write(r(tt):=0);
>>;
end;
%-----------------------------------

%-------------fft procedure-----------
%psf kernel dim = n
%------------------
procedure ffft(y,x);
begin

for uu:=0:(n-1) do <<

write( y(uu) := for tt:=0:(n-1) sum x(tt)*f(tt*uu) )

>>;
end;
%-----------------------------------


cconv();
cconvn();

for ii:=0:(m-1) do write("k(",ii,"):=",k(ii) );
for ii:=0:(n-1) do write("a(",ii,"):=",a(ii) );
for ii:=0:(n-1) do write("r(",ii,"):=",r(ii) );


%-----------------------------------

ffft(fa,a);
ffft(fk,k);
ffft(fr,r);

%-----------------------------------

for uu:=0:(n-1) do write(errr(uu):=fr(uu)-fa(uu)*fk(uu));


shut ccon;
bye;
end;