1 '''
2 Transform of length 2N real sequence using length N complex FFT.
3 Python 2.7.3
4 John Bryan, 2017
5 '''
6
7 import numpy as np
8 import time
9 import matplotlib.pyplot as plt
10 import warnings
11 np.set_printoptions(threshold=np.nan, precision=3, suppress=1)
12 warnings.filterwarnings("ignore")
13
14
15 def bitreversal(xarray, length, log2length):
16 '''
17 bracewell-buneman bit reversal function
18 inputs: xarray is array; length is array length; log2length=log2(length).
19 output: bit reversed array xarray.
20 '''
21 muplus = int((log2length+1)/2)
22 mvar = 1
23 reverse = np.zeros(length, dtype = int)
24 upper_range = muplus+1
25 for nuvar in np.arange(1, upper_range):
26 for kvar in np.arange(0, mvar):
27 tvar = 2*reverse[kvar]
28 reverse[kvar] = tvar
29 reverse[kvar+mvar] = tvar+1
30 mvar = mvar+mvar
31 if (int(log2length) & 0x01):
32 mvar = mvar/2
33 for qvar in np.arange(1, mvar):
34 nprime = qvar-mvar
35 rprimeprime = reverse[qvar]*mvar
36 for pvar in np.arange(0, reverse[qvar]):
37 nprime = nprime+mvar
38 rprime = rprimeprime+reverse[pvar]
39 temp = xarray[nprime]
40 xarray[nprime] = xarray[rprime]
41 xarray[rprime] = temp
42 return xarray
43
44
45 def fft0 (xarray, twiddle, log2length) :
46 '''
47 radix-2 dit
48 '''
49 nvar = xarray.size
50 b_p = nvar/2
51 nvar_p = 2
52 twiddle_step_size = nvar/2
53 for pvar in range(0, log2length):
54 nvar_pp = nvar_p/2
55 base_t = 0
56 for bvar in range(0, b_p):
57 base_b = base_t+nvar_pp
58 for nvar in range(0, nvar_pp):
59 if nvar == 0:
60 bot = xarray[base_b+nvar]
61 else:
62 twiddle_factor = nvar*twiddle_step_size
63 bot = xarray[base_b+nvar]*twiddle[twiddle_factor]
64 top = xarray[base_t+nvar]
65 xarray[base_t+nvar] = top+bot
66 xarray[base_b+nvar] = top-bot
67 base_t = base_t+nvar_p
68 b_p = b_p/2
69 nvar_p = nvar_p*2
70 twiddle_step_size = twiddle_step_size/2
71 return xarray
72
73
74 def steps(garray, halftwiddles, twiddles, length, log2length):
75 '''
76 steps function carries out the steps to calculate transform.
77 Input garray is real sequence.
78 Output g2array is transform of garray.
79 '''
80 n2var = length/2
81 g2array = np.zeros(length, dtype=np.complex_)
82 x1array = np.zeros(n2var, dtype=np.complex_)
83 x2array = np.zeros(n2var, dtype=np.complex_)
84 y1array = np.zeros(n2var, dtype=np.complex_ )
85 y2array = np.zeros(n2var, dtype=np.complex_ )
86
87 x1array = garray[::2]
88
89 x2array = garray[1::2]
90 x2j = [1j*mvar for mvar in x2array]
91 xarray = [avar+bvar for avar, bvar in zip(x1array, x2j)]
92 x0array = np.array(xarray)
93 x_array = x0array.astype(np.complex_)
94 halflength = int(x_array.size)
95 log2halflength = int(np.log2(halflength))
96 z0array = bitreversal(x_array, halflength, log2halflength)
97 x3array = fft0(z0array, halftwiddles, log2halflength)
98 x3conj = np.conjugate(x3array)
99 for kvar in range(1, n2var):
100 y1array[kvar] = 0.5*(x3array[kvar]+x3conj[n2var-kvar])
101 y2array[kvar] = -0.5*1j*(x3array[kvar]-x3conj[n2var-kvar])
102 y1array[0] = 0.5*(x3array[0]+x3conj[0])
103 y2array[0] = -0.5*1j*(x3array[0]-x3conj[0])
104 for kvar in range(0, n2var):
105 g2array[kvar] = y1array[kvar]+(twiddles[kvar]*y2array[kvar])
106 g2array[n2var] = 0.5*((x3array[0]+x3conj[0])+1j*(x3array[0]-x3conj[0]))
107 for kvar in range (1, n2var):
108 g2array[length-kvar] = np.conjugate(g2array[kvar])
109 return g2array
110
111
112
113 def plot_times(times2, times3):
114 '''
115 plot performance as a function of sequence length
116 '''
117 lengths = np.zeros(15, dtype=int)
118 for i in range(4, 19):
119 lengths[i-4] = np.power(2, i)
120 plt.figure(figsize=(7, 5))
121 plt.rc("font", size=9)
122 plt.loglog(lengths, times2, 'o', ms=5, markerfacecolor="None", markeredgecolor='blue', \
123 markeredgewidth=1, basex=2, basey=10, label='Half-length FFT plus additional steps')
124 plt.loglog(lengths, times3, '^', ms=5, markerfacecolor="None", markeredgecolor='green', \
125 markeredgewidth=1, basex=2, basey=10, label='Full-length complex FFT')
126 plt.legend(loc=2)
127 plt.grid()
128 plt.xlim([9, 530000])
129 plt.ylim([.00001, 8])
130 plt.ylabel("time (seconds)")
131 plt.xlabel("sequence length")
132 plt.title("Performance vs Sequence Length")
133 plt.savefig('t22.png', bbox_inches='tight')
134 plt.show()
135 return None
136
137
138 def test():
139 '''
140 Test w/ different length random sequences
141 '''
142 flag = 0
143 i = 0
144 times = np.zeros(15)
145 times2 = np.zeros(15)
146 for log2length in range (4, 19):
147 length = np.power(2, log2length)
148
149 np.random.seed(0)
150 garray = np.random.rand(length)
151 np.random.seed(0)
152 garray0 = np.random.rand(length)
153 np.random.seed(0)
154 garray00 = np.random.rand(length)
155
156 gpy = np.fft.fft(garray00)
157 n2var = garray.size/2
158 twiddles = np.zeros(n2var, dtype=np.complex)
159 halftwiddles = np.exp(-2j*np.pi*np.arange(0, 0.5, 1./(0.5*garray.size), dtype=np.complex_))
160 twiddles = np.exp(-2j*np.pi*np.arange(0, 0.5, 1./(garray.size), dtype=np.complex_))
161 t0time = time.time()
162 g2array = steps(garray, halftwiddles, twiddles, length, log2length)
163 times[i] = time.time()-t0time
164 t0time = time.time()
165 garray0 = bitreversal(garray0, length, log2length)
166 garray0 = garray0.astype(np.complex_)
167 g3array = fft0(garray0, twiddles, log2length)
168 times2[i] = time.time()-t0time
169 t_f = np.allclose(g3array, g2array)
170 if t_f == 0:
171 flag = 1
172
173 assert(t_f)
174 print t_f
175 i = i+1
176 if flag == 0:
177 print ("All results were correct.")
178 print 'Algorithm A (half-length complex fft with additional steps) fft time is', times
179 print 'Algorithm B (full-length complex fft) time is', times2
180 plot_times(times, times2)
181
182
183 test()