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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
c
c-----------------------------------------------------------------------
c main program: test program to exercise the wfta subroutine
c the test waveform is a complex exponential a**i whose
c transform is known analytically to be (1 - a**n)/(1 - a*w**k).
c
c authors:
c james h. mcclellan and hamid nawab
c department of electrical engineering and computer science
c massachusetts of technology
c cambridge, mass. 02139
c
c inputs:
c n-- transform length. it must be formed as the product of
c relatively prime integers from the set:
c 2,3,4,5,7,8,9,16
c invrs is the flag for forward or inverse transform.
c invrs = 1 yields inverse transform
c invrs .ne. 1 gives forward transform
c rad and phi are the magnitude and angle (as a fraction of
c 2*pi/n) of the complex exponential test signal.
c suggestion: rad = 0.98, phi = 0.5.
c-----------------------------------------------------------------------
c
double precision pi2,pin,xn,xj,xt
dimension xr(1260),xi(1260)
complex cone,ca,can,cnum,cden
c
c output will be punched
c
iout=i1mach(3)
input=i1mach(1)
cone=cmplx(1.0,0.0)
pi2=8.0d0*datan(1.0d0)
50 continue
read(input,130)n
130 format(i5)
write(iout,150) n
150 format(10h length = ,i5)
if(n.le.0 .or. n.gt.1260) stop
c
c enter a 1 to perform the inverse
c
c read(input,130) invrs
invrs = 0
c
c enter magnitude and angle (in fraction of 2*pi/n)
c avoid multiples of n for the angle if the radius is
c close to one. suggestion: rad = 0.98, phi = 0.5.
c
c read(input,160) rad,phi
rad = 0.98
phi = 0.5
160 format(2f15.10)
xn=float(n)
pin=phi
pin=pin*pi2/xn
c
c generate z**j
c
init=0
do 200 j=1,n
an=rad**(j-1)
xj=j-1
xj=xj*pin
xt=dcos(xj)
xr(j)=xt
xr(j)=xr(j)*an
xt=dsin(xj)
xi(j)=xt
xi(j)=xi(j)*an
200 continue
can=cmplx(xr(n),xi(n))
ca=cmplx(xr(2),xi(2))
can=can*ca
c
c print first 50 values of input sequence
c
max=50
if(n.lt.50)max=n
write(iout,300)(j,xr(j),xi(j),j=1,max)
c
c call the winograd fourier transform algorithm
c
call wfta(xr,xi,n,invrs,init,ierr)
c
c check for error return
c
if(ierr.lt.0) write(iout,250) ierr
250 format(1x,5herror,i5)
if(ierr.lt.0) go to 50
c
c print first 50 values of the transformed sequence
c
write(iout,300)(j,xr(j),xi(j),j=1,max)
300 format(1x,3hj =,i3,6hreal =,e20.12,6himag =,e20.12)
c
c calculate absolute and relative deviations
c
devabs=0.0
devrel=0.0
cnum=cone-can
pin=pi2/xn
do 350 j=1,n
xj=j-1
xj=-xj*pin
if(invrs.eq.1) xj=-xj
tr=dcos(xj)
ti=dsin(xj)
can=cmplx(tr,ti)
cden=cone-ca*can
cden=cnum/cden
c
c true value of the transform (1. - a**n)/(1. - a*w**k),
c where a = rad*exp(j*phi*(2*pi/n)), w = exp(-j*2*pi/n).
c for the inverse transform the complex exponential w
c is conjugated.
c
tr=real(cden)
ti=aimag(cden)
if(invrs.ne.1) go to 330
c
c scale inverse transform by 1/n
c
tr=tr/float(n)
ti=ti/float(n)
330 tr=xr(j)-tr
ti=xi(j)-ti
devabs=sqrt(tr*tr+ti*ti)
xmag=sqrt(xr(j)*xr(j)+xi(j)*xi(j))
devrel=100.0*devabs/xmag
if(devabs.le.devmx1)go to 340
devmx1=devabs
labs=j-1
340 if(devrel.le.devmx2)go to 350
devmx2=devrel
lrel=j-1
350 continue
c
c print the absolute and relative deviations together
c with their locations.
c
write(iout,380) devabs,labs,devrel,lrel
380 format(1x,21habsolute deviation = ,e20.12,9h at index,i5/
1 1x,21hrelative deviation = ,f11.7,8h percent,1x,9h at index,i5)
go to 50
end
|
