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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
c
c-----------------------------------------------------------------------
c main program: time-efficient radix-4 fast fourier transform
c author: l. robert morris
c department of systems engineering and computing science
c carleton university, ottawa, canada k1s 5b6
c
c input: the array "a" contains the data to be transformed
c-----------------------------------------------------------------------
c
c test program for autogen radix-4 fft
c
dimension a(2048),b(2048)
common /aa/a
c
ioutd=i1mach(2)
c
c compute dft and idft for n = 64, 256, and 1024 complex points
c
do 1 mm=3,5
n=4**mm
do 2 j=1,n
a(2*j-1)=uni(0)
a(2*j )=uni(0)
b(2*j-1)=a(2*j-1)
2 b(2*j )=a(2*j)
c
c forward dft
c
call radix4(mm,1,-1)
c
if(mm.ne.3) go to 5
c
c list dft input, output for n = 64 only
c
write(ioutd,98)
write(ioutd,100)
do 3 j=1,n
write(ioutd,96) b(2*j-1),b(2*j),a(2*j-1),a(2*j)
3 continue
c
c inverse dft
c
5 call radix4(mm,0, 1)
c
c list dft input, idft output for n = 64 only
c
if(mm.ne.3) go to 7
c
write(ioutd,99)
write(ioutd,100)
do 6 j=1,n
write(ioutd,96) b(2*j-1),b(2*j),a(2*j-1),a(2*j)
6 continue
c
c calculate rms error
c
7 err=0.0
do 8 j=1,n
8 err=err+(a(2*j-1)-b(2*j-1))**2+(a(2*j)-b(2*j))**2
err=sqrt(err/float(n))
write(ioutd,97) mm,err
1 continue
c
96 format(1x,4(f10.6,2x))
97 format(1x,20h rms error for m =,i2,4h is ,e14.6/)
98 format(1x,43h dft input dft output/)
99 format(1x,43h dft input idft output/)
100 format(1x,44h real imag real imag/)
stop
end
|
