1 '''
2 FHT
3 John Bryan, 2017
4 Python 2.7.3
5 '''
6
7
8 import numpy as np
9 import time
10 import matplotlib.pyplot as plt
11 import warnings
12 np.set_printoptions(threshold = np.nan, precision = 3, suppress = 1)
13 warnings.filterwarnings("ignore")
14
15
16 def swap(xarray, i, jvar):
17 '''
18 swap
19 '''
20 temp = xarray[i]
21 xarray[i] = xarray[jvar]
22 xarray[jvar] = temp
23 return None
24
25
26 def digitreversal(xarray, radix, log2length, length):
27 '''
28 digitreversal
29 '''
30 if log2length % 2 == 0:
31 n1var = int(np.sqrt(length))
32 else:
33 n1var = int(np.sqrt(int(length/radix)))
34
35 reverse = np.zeros(n1var, dtype = int)
36 reverse[1] = int(length/radix)
37 for jvar in range(1, radix):
38 reverse[jvar] = reverse[jvar-1]+reverse[1]
39 for i in range(1, int(n1var/radix)):
40 reverse[radix*i] = int(reverse[i]/radix)
41 for jvar in range(1, radix):
42 reverse[int(radix*i)+jvar] = reverse[int(radix*i)]+reverse[jvar]
43
44 for i in range(0, n1var-1):
45 for jvar in range(i+1, n1var):
46 uvar = i+reverse[jvar]
47 vvar = jvar+reverse[i]
48 swap(xarray, uvar, vvar)
49 if log2length % 2 == 1:
50 for zvar in range(1, radix):
51 uvar = i+reverse[jvar]+(zvar*n1var)
52 vvar = jvar+reverse[i]+(zvar*n1var)
53 swap(xarray, uvar, vvar)
54 return xarray
55
56
57 def fht(xvector, cosine, sine, length, log2length):
58 '''
59 fast hartley transform
60 '''
61 xarray = np.zeros(length)
62 b_p = length/2
63 n_p = 2
64 tss = length/2
65 for pvar in range(0, log2length):
66 npp = n_p/2
67 baset = 0
68 for bvar in range(0, b_p):
69 baseb = baset+npp
70 basebb = baset+n_p
71 for nvar in range(0, npp):
72 t_f = nvar*tss
73 nmn = (basebb-nvar)%basebb
74 if (pvar%2 == 0):
75 if (nvar == 0):
76 xcs = xvector[baseb+nvar]
77 xarray[baset+nvar] = xvector[baset+nvar]+xcs
78 xarray[baseb+nvar] = xvector[baset+nvar]-xcs
79 else:
80 xcs = (xvector[baseb+nvar]*cosine[t_f])+(xvector[nmn]*sine[t_f])
81 xarray[baset+nvar] = xvector[baset+nvar]+xcs
82 xarray[baseb+nvar] = xvector[baset+nvar]-xcs
83 else:
84 if (nvar == 0):
85 xcs = xarray[baseb+nvar]
86 xvector[baset+nvar] = xarray[baset+nvar]+xcs
87 xvector[baseb+nvar] = xarray[baset+nvar]-xcs
88 else:
89 xcs = (xarray[baseb+nvar]*cosine[t_f])+(xarray[nmn]*sine[t_f])
90 xvector[baset+nvar] = xarray[baset+nvar]+xcs
91 xvector[baseb+nvar] = xarray[baset+nvar]-xcs
92
93 baset = baset+n_p
94
95 b_p = b_p/2
96 n_p = n_p*2
97 tss = tss/2
98
99 if (pvar%2 == 0):
100 return xarray
101 else:
102 return xvector
103
104
105 def fft2 (xarray, twiddle, svar) :
106 '''
107 radix-2 dit
108 '''
109 nvar = xarray.size
110 b_p = nvar/2
111 nvar_p = 2
112 twiddle_step_size = nvar/2
113 for pvar in range(0, svar):
114 nvar_pp = nvar_p/2
115 base_t = 0
116 for bvar in range(0, b_p):
117 base_b = base_t+nvar_pp
118 for nvar in range(0, nvar_pp):
119 if nvar == 0:
120 bot = xarray[base_b+nvar]
121 else:
122 twiddle_factor = nvar*twiddle_step_size
123 bot = xarray[base_b+nvar]*twiddle[twiddle_factor]
124 top = xarray[base_t+nvar]
125 xarray[base_t+nvar] = top+bot
126 xarray[base_b+nvar] = top-bot
127 base_t = base_t+nvar_p
128 b_p = b_p/2
129 nvar_p = nvar_p*2
130 twiddle_step_size = twiddle_step_size/2
131 return xarray
132
133
134 def plot_times(times2, times3):
135 '''
136 plot fht and fft performance
137 '''
138 uvar = np.zeros(10, dtype = int)
139 for i in range(4, 14):
140 uvar[i-4] = np.power(2, i)
141 plt.figure(figsize = (7, 5))
142 plt.rc("font", size = 9)
143 plt.loglog(uvar, times3*1000, '^', ms = 5, markerfacecolor = "None", \
144 markeredgecolor = 'blue', markeredgewidth = 1, \
145 basex = 2, basey = 10, label = 'FFT')
146 plt.loglog(uvar, times2*1000, 'o', ms = 5, markerfacecolor = "None", \
147 markeredgecolor = 'black', markeredgewidth = 1, \
148 basex = 2, basey = 10, label = 'FHT')
149 plt.legend(loc = 2)
150 plt.grid()
151 plt.xlim([9, 16000])
152 plt.ylim([.1, 3000])
153 plt.ylabel("time (milliseconds)")
154 plt.xlabel("sequence length")
155 plt.title("Time vs Length for Python code")
156 plt.savefig('phartley.png', bbox_inches = 'tight')
157 plt.show()
158 return None
159
160
161 def test():
162 '''
163 Test w/ different length random sequences
164 '''
165 flag = 0
166 i = 0
167 times = np.zeros(10)
168 times2 = np.zeros(10)
169 for log2length in range (4, 14):
170 length = np.power(2, log2length)
171 np.random.seed(1)
172 garray = np.random.rand(length)
173 np.random.seed(1)
174 gfy = np.random.rand(length)
175 gpy = np.fft.fft(garray)
176 hpy = gpy.real-gpy.imag
177 gpy = gpy/length
178 hpy = hpy/length
179 sine = np.zeros(length, dtype = np.float)
180 cosine = np.zeros(length, dtype = np.float)
181 for index in range(0, length):
182 sine[index] = np.sin(2*np.pi*index/length)
183 cosine[index] = np.cos(2*np.pi*index/length)
184 kmax = (float(length)/2.)-1
185 kvar = np.linspace(0, kmax, kmax+1)
186 twiddles = np.exp(-2j*np.pi*kvar/length)
187 radix = 2
188 gfy = digitreversal(gfy, radix, log2length, length)
189 time0 = time.time()
190 gfy = fft2(gfy, twiddles, log2length)
191 times2[i] = time.time()-time0
192 gfy = gfy/length
193 garray = digitreversal(garray, radix, log2length, length)
194 time0 = time.time()
195 g2y = fht(garray, cosine, sine, length, log2length)
196 times[i] = time.time()-time0
197 g2y = g2y/length
198 t_f = np.allclose(g2y, hpy)
199 if t_f == 0:
200 flag = 1
201
202 assert(t_f)
203 i = i+1
204 if flag == 0:
205 print ("All results were correct.")
206 plot_times(times, times2)
207
208
209 test()
210
211
212