# This routine uses the fast fourier transform to compute the DFT. The input is a 
# vector [fh[1], ...fh[n]] of samples and the length, n of this vector. The routine finds 
# the smallest power of 2 larger than n and computes a FFT of that length with zero 
# padding to fill in the missing entries. It returns a vector whose first entry is N the 
# length of the FFT, the next N entries are the FFT itself.
mFFT := proc(fh,n) local i,j,k,a,b,M,N:  
M:=ceil(evalf(log[2](n))):  
N := 2^(M): 
a := hfarray(1..N): b := hfarray(1..N):
for i from 1 to n do a[i] := evalf(Re(fh[i])): b[i] := evalf(Im(fh[i])): od: 
if n< N then for k from 1 to N-n do a[k+n] :=0.0: b[k+n] := 0.0 od: end if: 
evalhf(FFT(M,var(a),var(b))): 
[seq( a[k]+I*b[k], k=1..N)]: 
end:
# This routine inverts the discrete Fourier transform. The input is the vector 
# [[fh[1],...,fh[2^n]] followed by n where 2^n is the length of this vector, the output is 
# the inverse Fourier transform of fh.
miFFT := proc(fh,n) local a,b,c,i,j,k;
a := hfarray(1..2^n); b := hfarray(1..2^n);
for i from 1 to 2^n do a[i] := evalf(Re(fh[i])); b[i] := evalf(Im(fh[i]));  od; 
evalhf(iFFT(n,var(a),var(b))); 
[seq( a[k]+I*b[k], k=1..2^n)];
end:
#This routine applies a filter with discretized transfer function [h[1],...,h[2^n]] to a 
# function whose 2^n point FFT is given in [f[1],...f[2^n]].
fil := proc(f,h,n) local a, b,c,i,j,k;
c := [seq(f[j]*h[j],j=1..2^n)];
a := hfarray(1..2^n); b := hfarray(1..2^n);
for i from 1 to 2^n do a[i] := evalf(Re(c[i])); b[i] := evalf(Im(c[i])); od;
evalhf(iFFT(n,var(a),var(b)));
[seq( a[k]+I*b[k], k=1..2^n)];
end: