1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
|
if not ?fboundp (fft) then load (fft);
"The fft of a list:"$
(N:32,
g:makelist(if i<N/2 then 1.0 else 0.0,i,0,N-1),
plot2d([discrete,makelist([i,g[i]],i,1,N)],['y,-0.1,1.1]),
gt:fft(g));
"The inverse_fft of a list:"$
gti:inverse_fft(gt);
"The norm-difference between list G and output of inverse_fft"$
apply(max,map(cabs,gti-g));
"The fft of an array:"$
(N:8,
g:g+%i*g,
ga:make_array(any,N), fillarray(ga,g),
gta:fft(ga));
"The inverse_fft of an array:"$
gtia:inverse_fft(gta);
"The norm-difference between array GA and output of inverse_fft"$
(mx:0, for i:0 thru N-1 do mx:max(mx,cabs(gtia[i]-ga[i])), mx);
"Example from info pages:"$
L : [1, 2, 3, 4, 5, 6, 7, 8];
(n : length (L),
x : make_array (any, n),
fillarray (x, L) );
y : fft (x) ;
(a : make_array (any, n/2 + 1) ,
b : make_array (any, n/2 + 1) ,
a[0] : realpart (y[0]) ,
b[0] : 0 ,
for k : 1 thru n/2 - 1 do
(a[k] : realpart (y[k] + y[n - k]),
b[k] : imagpart (y[n - k] - y[k])),
a[n/2] : y[n/2] ,
b[n/2] : 0 )$
listarray (a);
listarray (b);
f(j) := sum (a[k]*cos(-2*%pi*j*k/n) + b[k]*sin(-2*%pi*j*k/n), k, 0, n/2) $
"The following should equal the list L"$
makelist (float(f (j)), j, 0, n - 1);
|
