1 # DCT "REMAP" scheduler
3 # Modifications made to create an in-place iterative DCT:
4 # Copyright (c) 2021 Luke Kenneth Casson Leighton <lkcl@lkcl.net>
8 # Original fastdctlee.py by Nayuki:
9 # Copyright (c) 2020 Project Nayuki. (MIT License)
10 # https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms
14 # bits of the integer 'val'.
15 def reverse_bits(val
, width
):
17 for _
in range(width
):
18 result
= (result
<< 1) |
(val
& 1)
23 # iterative version of [recursively-applied] half-rev.
24 # relies on the list lengths being power-of-two and the fact
25 # that bit-inversion of a list of binary numbers is the same
26 # as reversing the order of the list
27 # this version is dead easy to implement in hardware.
28 # a big surprise is that the half-reversal can be done with
29 # such a simple XOR. the inverse operation is slightly trickier
30 def halfrev2(vec
, pre_rev
=True):
32 for i
in range(len(vec
)):
34 res
.append(i ^
(i
>>1))
38 for ji
in range(1, bl
):
44 # python "yield" can be iterated. use this to make it clear how
45 # the indices are generated by using natural-looking nested loops
46 def iterate_dct_butterfly_indices(SVSHAPE
):
47 # get indices to iterate over, in the required order
49 # createing lists of indices to iterate over in each dimension
50 # has to be done dynamically, because it depends on the size
51 # first, the size-based loop (which can be done statically)
57 # invert order if requested
64 # reference (read/write) the in-place data in *reverse-bit-order*
66 if SVSHAPE
.mode
== 0b01:
67 levels
= n
.bit_length() - 1
68 ri
= [ri
[reverse_bits(i
, levels
)] for i
in range(n
)]
70 # reference list for not needing to do data-swaps, just swap what
71 # *indices* are referenced (two levels of indirection at the moment)
72 # pre-reverse the data-swap list so that it *ends up* in the order 0123..
74 inplace_mode
= SVSHAPE
.mode
== 0b01 and SVSHAPE
.skip
not in [0b10, 0b11]
76 print ("inplace mode")
77 ji
= halfrev2(ji
, True)
82 # start an infinite (wrapping) loop
85 for size
in x_r
: # loop over 3rd order dimension (size)
86 x_end
= size
== x_r
[-1]
87 # y_r schedule depends on size
90 for i
in range(0, n
, size
):
93 if SVSHAPE
.invxyz
[1]: y_r
.reverse()
94 for i
in y_r
: # loop over 2nd order dimension
96 # two lists of half-range indices, e.g. j 0123, jr 7654
97 j
= list(range(i
, i
+ halfsize
))
98 jr
= list(range(i
+halfsize
, i
+ size
))
100 # invert if requested
101 if SVSHAPE
.invxyz
[2]: k_r
.reverse()
102 if SVSHAPE
.invxyz
[2]: j_r
.reverse()
103 hz2
= halfsize
// 2 # zero stops reversing 1-item lists
104 # if you *really* want to do the in-place swapping manually,
105 # this allows you to do it. good luck...
106 if SVSHAPE
.mode
== 0b01 and not inplace_mode
:
109 print ("xform jr", jr
)
110 for jl
, jh
in zip(j
, jr
): # loop over 1st order dimension
112 # now depending on MODE return the index
113 if SVSHAPE
.skip
in [0b00, 0b10]:
114 result
= ri
[ji
[jl
]] # lower half
115 elif SVSHAPE
.skip
in [0b01, 0b11]:
116 result
= ri
[ji
[jh
]] # upper half, reverse order
118 ((y_end
and z_end
)<<1) |
119 ((y_end
and x_end
and z_end
)<<2))
121 yield result
+ SVSHAPE
.offset
, loopends
124 if SVSHAPE
.mode
== 0b01 and inplace_mode
:
125 for ci
, (jl
, jh
) in enumerate(zip(j
[:hz2
], jr
[:hz2
])):
127 #print ("inplace swap", jh, jlh)
128 tmp1
, tmp2
= ji
[jlh
], ji
[jh
]
129 ji
[jlh
], ji
[jh
] = tmp2
, tmp1
132 def pprint_schedule(schedule
, n
):
137 tablestep
= n
// size
138 print ("size %d halfsize %d tablestep %d" % \
139 (size
, halfsize
, tablestep
))
140 for i
in range(0, n
, size
):
141 prefix
= "i %d\t" % i
142 for j
in range(i
, i
+ halfsize
):
143 (jl
, je
), (jh
, he
) = schedule
[idx
]
144 print (" %-3d\t%s j=%-2d jh=%-2d "
145 "j[jl=%-2d] j[jh=%-2d]" % \
146 (idx
, prefix
, j
, j
+halfsize
,
149 "end", bin(je
)[2:], bin(je
)[2:])
153 # totally cool *in-place* DCT algorithm using yield REMAPs
159 print ("transform2", n
)
160 levels
= n
.bit_length() - 1
162 # reference (read/write) the in-place data in *reverse-bit-order*
164 ri
= [ri
[reverse_bits(i
, levels
)] for i
in range(n
)]
166 # and pretend we LDed data in half-swapped *and* bit-reversed order as well
167 # TODO: merge these two
168 vec
= halfrev2(vec
, False)
169 vec
= [vec
[ri
[i
]] for i
in range(n
)]
171 # create a cos table: not strictly necessary but here for illustrative
172 # purposes, to demonstrate the point that it really *is* iterative.
173 # this table could be cached and used multiple times rather than
174 # computed every time.
180 for i
in range(n
//size
):
181 for ci
in range(halfsize
):
182 ctable
.append((math
.cos((ci
+ 0.5) * math
.pi
/ size
) * 2.0))
198 SVSHAPE0
.lims
= [xdim
, ydim
, zdim
]
199 SVSHAPE0
.order
= [0,1,2] # experiment with different permutations, here
202 SVSHAPE0
.offset
= 0 # experiment with different offset, here
203 SVSHAPE0
.invxyz
= [1,0,0] # inversion if desired
204 # j+halfstep schedule
206 SVSHAPE1
.lims
= [xdim
, ydim
, zdim
]
207 SVSHAPE1
.order
= [0,1,2] # experiment with different permutations, here
210 SVSHAPE1
.offset
= 0 # experiment with different offset, here
211 SVSHAPE1
.invxyz
= [1,0,0] # inversion if desired
213 # enumerate over the iterator function, getting new indices
214 i0
= iterate_dct_butterfly_indices(SVSHAPE0
)
215 i1
= iterate_dct_butterfly_indices(SVSHAPE1
)
216 for k
, ((jl
, jle
), (jh
, jhe
)) in enumerate(zip(i0
, i1
)):
219 t1
, t2
= vec
[jl
], vec
[jh
]
222 vec
[jh
] = (t1
- t2
) * (1/coeff
)
223 print ("coeff", size
, i
, "ci", ci
,
225 "i/n", (ci
+0.5)/size
, coeff
, vec
[jl
],
228 print("transform2 pre-itersum", vec
)
230 # now things are in the right order for the outer butterfly.
235 ir
= list(range(0, halfsize
))
236 print ("itersum", halfsize
, size
, ir
)
238 jr
= list(range(i
+halfsize
, i
+n
-halfsize
, size
))
239 print ("itersum jr", i
+halfsize
, i
+size
, jr
)
241 vec
[jh
] += vec
[jh
+size
]
242 print (" itersum", size
, i
, jh
, jh
+size
)
245 print("transform2 result", vec
)
251 # set the dimension sizes here
253 ydim
= 0 # not needed
254 zdim
= 0 # again, not needed
256 # set total. err don't know how to calculate how many there are...
257 # do it manually for now
263 tablestep
= n
// size
264 for i
in range(0, n
, size
):
265 for j
in range(i
, i
+ halfsize
):
278 SVSHAPE0
.lims
= [xdim
, ydim
, zdim
]
279 SVSHAPE0
.order
= [0,1,2] # experiment with different permutations, here
282 SVSHAPE0
.offset
= 0 # experiment with different offset, here
283 SVSHAPE0
.invxyz
= [0,0,0] # inversion if desired
284 # j+halfstep schedule
286 SVSHAPE1
.lims
= [xdim
, ydim
, zdim
]
287 SVSHAPE1
.order
= [0,1,2] # experiment with different permutations, here
290 SVSHAPE1
.offset
= 0 # experiment with different offset, here
291 SVSHAPE1
.invxyz
= [0,0,0] # inversion if desired
293 # enumerate over the iterator function, getting new indices
295 i0
= iterate_dct_butterfly_indices(SVSHAPE0
)
296 i1
= iterate_dct_butterfly_indices(SVSHAPE1
)
297 for idx
, (jl
, jh
) in enumerate(zip(i0
, i1
)):
300 schedule
.append((jl
, jh
))
302 # ok now pretty-print the results, with some debug output
303 print ("inner butterfly")
304 pprint_schedule(schedule
, n
)
313 SVSHAPE0
.lims
= [xdim
, ydim
, zdim
]
314 SVSHAPE0
.order
= [0,1,2] # experiment with different permutations, here
317 SVSHAPE0
.offset
= 0 # experiment with different offset, here
318 SVSHAPE0
.invxyz
= [1,0,0] # inversion if desired
319 # j+halfstep schedule
321 SVSHAPE1
.lims
= [xdim
, ydim
, zdim
]
322 SVSHAPE1
.order
= [0,1,2] # experiment with different permutations, here
325 SVSHAPE1
.offset
= 0 # experiment with different offset, here
326 SVSHAPE1
.invxyz
= [1,0,0] # inversion if desired
328 # enumerate over the iterator function, getting new indices
330 i0
= iterate_dct_butterfly_indices(SVSHAPE0
)
331 i1
= iterate_dct_butterfly_indices(SVSHAPE1
)
332 for idx
, (jl
, jh
) in enumerate(zip(i0
, i1
)):
335 schedule
.append((jl
, jh
))
337 # ok now pretty-print the results, with some debug output
338 print ("outer butterfly")
339 pprint_schedule(schedule
, n
)
342 if __name__
== '__main__':