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 ji
= halfrev2(ji
, True)
78 # start an infinite (wrapping) loop
81 for size
in x_r
: # loop over 3rd order dimension (size)
82 x_end
= size
== x_r
[-1]
83 # y_r schedule depends on size
86 for i
in range(0, n
, size
):
89 if SVSHAPE
.invxyz
[1]: y_r
.reverse()
90 for i
in y_r
: # loop over 2nd order dimension
95 for j
in range(i
, i
+halfsize
):
99 if SVSHAPE
.invxyz
[2]: k_r
.reverse()
100 if SVSHAPE
.invxyz
[2]: j_r
.reverse()
101 hz2
= halfsize
// 2 # zero stops reversing 1-item lists
102 # if you *really* want to do the in-place swapping manually,
103 # this allows you to do it. good luck...
104 if SVSHAPE
.mode
== 0b01 and not inplace_mode
:
106 for j
in j_r
: # loop over 1st order dimension
108 # now depending on MODE return the index
109 if SVSHAPE
.skip
in [0b00, 0b10]:
110 result
= ri
[ji
[j
]] # lower half
111 elif SVSHAPE
.skip
in [0b01, 0b11]:
112 result
= ri
[ji
[size
-j
-1]] # upper half, reverse order
114 ((y_end
and z_end
)<<1) |
115 ((y_end
and x_end
and z_end
)<<2))
117 yield result
+ SVSHAPE
.offset
, loopends
120 if SVSHAPE
.mode
== 0b01 and inplace_mode
:
121 for ci
, jl
in enumerate(j_r
[:hz2
]):
122 jh
= size
-jl
-1 # upper half, reverse order
124 tmp1
, tmp2
= ji
[jlh
], ji
[jh
]
125 ji
[jlh
], ji
[jh
] = tmp2
, tmp1
128 def pprint_schedule(schedule
, n
):
133 tablestep
= n
// size
134 print ("size %d halfsize %d tablestep %d" % \
135 (size
, halfsize
, tablestep
))
136 for i
in range(0, n
, size
):
137 prefix
= "i %d\t" % i
138 for j
in range(i
, i
+ halfsize
):
139 (jl
, je
), (jh
, he
) = schedule
[idx
]
140 print (" %-3d\t%s j=%-2d jh=%-2d "
141 "j[jl=%-2d] j[jh=%-2d]" % \
142 (idx
, prefix
, j
, j
+halfsize
,
145 "end", bin(je
)[2:], bin(je
)[2:])
151 # set the dimension sizes here
153 ydim
= 0 # not needed
154 zdim
= 0 # again, not needed
156 # set total. err don't know how to calculate how many there are...
157 # do it manually for now
163 tablestep
= n
// size
164 for i
in range(0, n
, size
):
165 for j
in range(i
, i
+ halfsize
):
178 SVSHAPE0
.lims
= [xdim
, ydim
, zdim
]
179 SVSHAPE0
.order
= [0,1,2] # experiment with different permutations, here
182 SVSHAPE0
.offset
= 0 # experiment with different offset, here
183 SVSHAPE0
.invxyz
= [0,0,0] # inversion if desired
184 # j+halfstep schedule
186 SVSHAPE1
.lims
= [xdim
, ydim
, zdim
]
187 SVSHAPE1
.order
= [0,1,2] # experiment with different permutations, here
190 SVSHAPE1
.offset
= 0 # experiment with different offset, here
191 SVSHAPE1
.invxyz
= [0,0,0] # inversion if desired
193 # enumerate over the iterator function, getting new indices
195 i0
= iterate_dct_butterfly_indices(SVSHAPE0
)
196 i1
= iterate_dct_butterfly_indices(SVSHAPE1
)
197 for idx
, (jl
, jh
) in enumerate(zip(i0
, i1
)):
200 schedule
.append((jl
, jh
))
202 # ok now pretty-print the results, with some debug output
203 print ("inner butterfly")
204 pprint_schedule(schedule
, n
)
213 SVSHAPE0
.lims
= [xdim
, ydim
, zdim
]
214 SVSHAPE0
.order
= [0,1,2] # experiment with different permutations, here
217 SVSHAPE0
.offset
= 0 # experiment with different offset, here
218 SVSHAPE0
.invxyz
= [1,0,0] # inversion if desired
219 # j+halfstep schedule
221 SVSHAPE1
.lims
= [xdim
, ydim
, zdim
]
222 SVSHAPE1
.order
= [0,1,2] # experiment with different permutations, here
225 SVSHAPE1
.offset
= 0 # experiment with different offset, here
226 SVSHAPE1
.invxyz
= [1,0,0] # inversion if desired
228 # enumerate over the iterator function, getting new indices
230 i0
= iterate_dct_butterfly_indices(SVSHAPE0
)
231 i1
= iterate_dct_butterfly_indices(SVSHAPE1
)
232 for idx
, (jl
, jh
) in enumerate(zip(i0
, i1
)):
235 schedule
.append((jl
, jh
))
237 # ok now pretty-print the results, with some debug output
238 print ("outer butterfly")
239 pprint_schedule(schedule
, n
)
242 if __name__
== '__main__':