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_inner_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
58 if SVSHAPE
.invxyz
[0]: x_r
.reverse()
63 # reference (read/write) the in-place data in *reverse-bit-order*
65 ri
= [ri
[reverse_bits(i
, levels
)] for i
in range(n
)]
67 # reference list for not needing to do data-swaps, just swap what
68 # *indices* are referenced (two levels of indirection at the moment)
69 # pre-reverse the data-swap list so that it *ends up* in the order 0123..
71 inplace_mode
= SVSHAPE
.skip
not in [0b10, 0b11]
73 ji
= halfrev2(ji
, True)
75 # start an infinite (wrapping) loop
78 for size
in x_r
: # loop over 3rd order dimension (size)
79 x_end
= size
== x_r
[-1]
80 # y_r schedule depends on size
83 for i
in range(0, n
, size
):
86 if SVSHAPE
.invxyz
[1]: y_r
.reverse()
87 for i
in y_r
: # loop over 2nd order dimension
92 for j
in range(i
, i
+halfsize
):
96 if SVSHAPE
.invxyz
[2]: k_r
.reverse()
97 if SVSHAPE
.invxyz
[2]: j_r
.reverse()
98 hz2
= halfsize
// 2 # zero stops reversing 1-item lists
99 # if you *really* want to do the in-place swapping manually,
100 # this allows you to do it. good luck...
103 for j
in j_r
: # loop over 1st order dimension
105 # now depending on MODE return the index
106 if SVSHAPE
.skip
in [0b00, 0b10]:
107 result
= ri
[ji
[j
]] # lower half
108 elif SVSHAPE
.skip
in [0b01, 0b11]:
109 result
= ri
[ji
[size
-j
-1]] # upper half, reverse order
111 ((y_end
and z_end
)<<1) |
112 ((y_end
and x_end
and z_end
)<<2))
114 yield result
+ SVSHAPE
.offset
, loopends
118 for ci
, (jl
, jh
) in enumerate(zip(j
[:hz2
], jr
[:hz2
])):
120 tmp1
, tmp2
= ji
[jlh
], ji
[jh
]
121 ji
[jlh
], ji
[jh
] = tmp2
, tmp1
124 # totally cool *in-place* DCT algorithm
129 levels
= n
.bit_length() - 1
131 # and pretend we LDed data in half-swapped *and* bit-reversed order as well
132 # TODO: merge these two
133 vec
= halfrev2(vec
, False)
134 vec
= [vec
[ri
[i
]] for i
in range(n
)]
136 # start the inner butterfly
140 tablestep
= n
// size
141 ir
= list(range(0, n
, size
))
143 # two lists of half-range indices, e.g. j 0123, jr 7654
144 j
= list(range(i
, i
+ halfsize
))
145 jr
= list(range(i
+halfsize
, i
+ size
))
147 for ci
, (jl
, jh
) in enumerate(zip(j
, jr
)):
148 vec
[ri
[ji
[jl
]]] = t1
+ t2
149 vec
[ri
[ji
[jh
]]] = (t1
- t2
) * (1/coeff
)
150 hz2
= halfsize
// 2 # can be zero which stops reversing 1-item lists
151 for ci
, (jl
, jh
) in enumerate(zip(j
[:hz2
], jr
[:hz2
])):
153 tmp1
, tmp2
= ji
[jlh
], ji
[jh
]
154 ji
[jlh
], ji
[jh
] = tmp2
, tmp1
157 def dct_outer_butterfly(SVSHAPE
):
162 ir
= list(range(0, halfsize
))
163 print ("itersum", halfsize
, size
, ir
)
165 jr
= list(range(i
+halfsize
, i
+n
-halfsize
, size
))
166 print ("itersum jr", i
+halfsize
, i
+size
, jr
)
168 vec
[jh
] += vec
[jh
+size
]
169 print (" itersum", size
, i
, jh
, jh
+size
)