1 '''
2 Radix-4 N=64 DIT FFT with reduced twiddle table size of 9
3 John Bryan, 2017
4 Python 2.7.3
5 '''
6
7
8 import numpy as np
9
10
11 def tobinary(dec):
12 '''
13 decimal to binary
14 '''
15 binary = ''
16 dividend = dec
17 while dividend != 0:
18 binary = binary+str(dividend%2)
19 dividend = dividend/2
20 return binary[::-1]
21
22
23 def todec(binary):
24 '''
25 binary to decimal
26 '''
27 i = 0
28 decimal = 0
29 while (i<4):
30 decimal = decimal+(binary[i]*np.power(2, 3-i))
31 i = i+1
32 return decimal
33
34
35 def adder(avar, q0_):
36 '''
37 binary adder
38 '''
39 aavar = np.array([0, 0, 0])
40 cvar = np.zeros(4, dtype=int)
41 bvar = q0_
42 i = 0
43 while (i<3):
44 aavar[i] = avar[i]^q0_
45 i = i+1
46 i = 3
47 while (i>0):
48 cvar[i] = aavar[i-1]^bvar
49 bvar = aavar[i-1]&bvar
50 i = i-1
51 cvar[0] = bvar
52 return todec(cvar)
53
54
55 def digitreversal(xarray, radix, log2length, length):
56 '''
57 Yu-Loh-Miller digit-reversal algorithm
58 inputs: xarray is input array; radix = 2 for single binary bit reversal;
59 log2length = log2(xarray length); length = xarray length.
60 output: digit-reversed array xarray.
61 '''
62 if log2length % 2 == 0:
63 n1var = int(np.sqrt(length))
64 else:
65 n1var = int(np.sqrt(int(length/radix)))
66
67 reverse = np.zeros(n1var, dtype = int)
68 reverse[1] = int(length/radix)
69 for jvar in range(1, radix):
70 reverse[jvar] = reverse[jvar-1]+reverse[1]
71 for i in range(1, int(n1var/radix)):
72 reverse[radix*i] = int(reverse[i]/radix)
73 for jvar in range(1, radix):
74 reverse[int(radix*i)+jvar] = reverse[int(radix*i)]+reverse[jvar]
75
76 for i in range(0, n1var-1):
77 for jvar in range(i+1, n1var):
78 uvar = i+reverse[jvar]
79 vvar = jvar+reverse[i]
80
81 temp = xarray[uvar]
82 xarray[uvar] = xarray[vvar]
83 xarray[vvar] = temp
84 if log2length % 2 == 1:
85 for zvar in range(1, radix):
86 uvar = i+reverse[jvar]+(zvar*n1var)
87 vvar = jvar+reverse[i]+(zvar*n1var)
88
89 temp = xarray[uvar]
90 xarray[uvar] = xarray[vvar]
91 xarray[vvar] = temp
92 return xarray
93
94
95 def twiddlefactor(ind, mw_):
96 '''
97 twiddle factor
98 '''
99 index = tobinary(ind)
100 index = index.rjust(6, '0')
101 q0var = int(index[2])
102 q1var = int(index[1])
103 q2var = int(index[0])
104 rvar = index[3:6]
105 qvar = np.zeros(3, dtype=int)
106 for svar in range(0, 3):
107 qvar[svar] = int(rvar[svar])
108 index0 = adder(qvar, q0var)
109 mvar = mw_[index0]
110 q0q1 = q0var^q1var
111 q0q2 = q0var^q2var
112 q0q1q2 = q0var^q1var^q2var
113 rm_ = np.real(mvar)
114 im_ = np.imag(mvar)
115 if q0q2 == 0:
116 pq0q2 = 1
117 else:
118 pq0q2 = -1
119 if q0q1q2 == 0:
120 pq0q1q2 = 1
121 else:
122 pq0q1q2 = -1
123 pq0q2rm = pq0q2*rm_
124 pq0q2im = pq0q2*im_
125 jpq0q1q2im = 1j*pq0q1q2*im_
126 jpq0q1q2rm = 1j*pq0q1q2*rm_
127 q0q1q2 = q0var^q1var^q2var
128 if q0q1 == 0:
129 wvar = pq0q2rm+jpq0q1q2im
130 else:
131 wvar = pq0q2im+jpq0q1q2rm
132 return wvar
133
134
135 def fft4(xarray, mw_, svar):
136 '''
137 radix-4 dit fft
138 '''
139 nvar = 4
140 tss = np.power(4, svar-1)
141 krange = 1
142 block = int(xarray.size/4)
143 base = 0
144 for wvar in range(0, svar):
145 for zvar in range(0, block):
146 for kvar in range(0, krange):
147
148 offset = nvar/4
149 avar = base+kvar
150 bvar = base+kvar+offset
151 cvar = base+kvar+(2*offset)
152 dvar = base+kvar+(3*offset)
153 if kvar == 0:
154 xbr1 = xarray[bvar]
155 xcr2 = xarray[cvar]
156 xdr3 = xarray[dvar]
157 else:
158 index1 = int(kvar*tss)
159 index2 = int(2*kvar*tss)
160 index3 = int(3*kvar*tss)
161 r1var = twiddlefactor(index1, mw_)
162 r2var = twiddlefactor(index2, mw_)
163 r3var = twiddlefactor(index3, mw_)
164 xbr1 = (xarray[bvar]*r1var)
165 xcr2 = (xarray[cvar]*r2var)
166 xdr3 = (xarray[dvar]*r3var)
167 evar = xarray[avar]+xcr2
168 fvar = xarray[avar]-xcr2
169 gvar = xbr1+xdr3
170 hvar = xbr1-xdr3
171 j_h = 1j*hvar
172 xarray[avar] = evar+gvar
173 xarray[bvar] = fvar-j_h
174 xarray[cvar] = -gvar+evar
175 xarray[dvar] = j_h+fvar
176 base = base+(4*krange)
177 block = block/4
178 nvar = 4*nvar
179 krange = 4*krange
180 base = 0
181 tss = float(tss)/4.
182 return xarray
183
184
185
186 def test():
187 '''
188 test for correctness
189 '''
190 flag = 0
191 log4length = 3
192 xarray = np.random.rand(2*np.power(4, log4length)).view(np.complex128)
193 xpy = np.fft.fft(xarray)
194 radix = 4
195 length = np.power(4, log4length)
196 xarray = digitreversal(xarray, radix, log4length, length)
197 hvar = np.linspace(0, 8./float(length), 9)
198 mw_ = np.exp(-2j*np.pi*hvar)
199 xarray = fft4(xarray, mw_, log4length)
200 t_f = np.allclose(xarray, xpy)
201 if t_f == 0:
202 flag = 1
203 assert(t_f)
204 if flag == 0:
205 print ("All radix-4 results were correct.")
206 return None
207
208
209 test()