src/openpower/decoder/isa/fastdctlee.py

   1 #
   2 # Fast discrete cosine transform algorithms (Python)
   3 #
   4 # Copyright (c) 2020 Project Nayuki. (MIT License)
   5 # https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms
   6 #
   7 # Permission is hereby granted, free of charge, to any person obtaining a copy of
   8 # this software and associated documentation files (the "Software"), to deal in
   9 # the Software without restriction, including without limitation the rights to
  10 # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
  11 # the Software, and to permit persons to whom the Software is furnished to do so,
  12 # subject to the following conditions:
  13 # - The above copyright notice and this permission notice shall be included in
  14 #   all copies or substantial portions of the Software.
  15 # - The Software is provided "as is", without warranty of any kind, express or
  16 #   implied, including but not limited to the warranties of merchantability,
  17 #   fitness for a particular purpose and noninfringement. In no event shall the
  18 #   authors or copyright holders be liable for any claim, damages or other
  19 #   liability, whether in an action of contract, tort or otherwise, arising from,
  20 #   out of or in connection with the Software or the use or other dealings in the
  21 #   Software.
  22 #
  23
  24 import math
  25 from copy import deepcopy
  26
  27
  28 # DCT type II, unscaled. Algorithm by Byeong Gi Lee, 1984.
  29 # See: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.118.3056&rep=rep1&type=pdf#page=34
  30 def transform(vector, indent=0):
  31     idt = "   " * indent
  32     n = len(vector)
  33     if n == 1:
  34         return list(vector)
  35     elif n == 0 or n % 2 != 0:
  36         raise ValueError()
  37     else:
  38         half = n // 2
  39         alpha = [(vector[i] + vector[-(i + 1)]) for i in range(half)]
  40         beta  = [(vector[i] - vector[-(i + 1)]) /
  41                  (math.cos((i + 0.5) * math.pi / n) * 2.0)
  42                     for i in range(half)]
  43         alpha = transform(alpha)
  44         beta  = transform(beta )
  45         result = []
  46         for i in range(half - 1):
  47             result.append(alpha[i])
  48             result.append(beta[i] + beta[i + 1])
  49         result.append(alpha[-1])
  50         result.append(beta [-1])
  51         return result
  52
  53
  54 def transform(vector, indent=0):
  55     idt = "   " * indent
  56     n = len(vector)
  57     if n == 1:
  58         return list(vector)
  59     elif n == 0 or n % 2 != 0:
  60         raise ValueError()
  61     else:
  62         half = n // 2
  63         alpha = [0] * half
  64         beta = [0] * half
  65         print (idt, "xf", vector)
  66         print (idt, "coeff", n, "->", end=" ")
  67         for i in range(half):
  68             t1, t2 = vector[i], vector[n-i-1]
  69             k = (math.cos((i + 0.5) * math.pi / n) * 2.0)
  70             print (i, n-i-1, "i/n", (i+0.5)/n, ":", k, end= " ")
  71             alpha[i] = t1 + t2
  72             beta[i] = (t1 - t2) * (1/k)
  73         print ()
  74         print (idt, "n", n, "alpha", end=" ")
  75         for i in range(0, n, 2):
  76             print (i, i//2, alpha[i//2], end=" ")
  77         print()
  78         print (idt, "n", n, "beta", end=" ")
  79         for i in range(0, n, 2):
  80             print (i, beta[i//2], end=" ")
  81         print()
  82         alpha = transform(alpha, indent+1)
  83         beta  = transform(beta , indent+1)
  84         result = [0] * n
  85         for i in range(half):
  86             result[i*2] = alpha[i]
  87             result[i*2+1] = beta[i]
  88         print(idt, "merge", result)
  89         for i in range(half - 1):
  90             result[i*2+1] += result[i*2+3]
  91         print(idt, "result", result)
  92         return result
  93
  94
  95 def transform_itersum(vector, indent=0):
  96     idt = "   " * indent
  97     n = len(vector)
  98     if n == 1:
  99         return list(vector)
 100     elif n == 0 or n % 2 != 0:
 101         raise ValueError()
 102     else:
 103         half = n // 2
 104         alpha = [0] * half
 105         beta = [0] * half
 106         for i in range(half):
 107             t1, t2 = vector[i], vector[i+half]
 108             alpha[i] = t1
 109             beta[i] = t2
 110         alpha = transform_itersum(alpha, indent+1)
 111         beta  = transform_itersum(beta , indent+1)
 112         result = [0] * n
 113         for i in range(half):
 114             result[i*2] = alpha[i]
 115             result[i*2+1] = beta[i]
 116         print(idt, "iter-merge", result)
 117         for i in range(half - 1):
 118             result[i*2+1] += result[i*2+3]
 119         print(idt, "iter-result", result)
 120         return result
 121
 122
 123
 124 def transform2(vec, reverse=True):
 125
 126     vec = deepcopy(vec)
 127     # bits of the integer 'val'.
 128     def reverse_bits(val, width):
 129         result = 0
 130         for _ in range(width):
 131             result = (result << 1) | (val & 1)
 132             val >>= 1
 133         return result
 134
 135     # Initialization
 136     n = len(vec)
 137     print ("transform2", n)
 138     levels = n.bit_length() - 1
 139
 140     size = n
 141     while size >= 2:
 142         halfsize = size // 2
 143         tablestep = n // size
 144         ir = list(range(0, n, size))
 145         print ("  xform", size, ir)
 146         for i in ir:
 147             k = 0
 148             j = list(range(i, i + halfsize))
 149             jr = list(range(i+halfsize, i + size))
 150             jr.reverse()
 151             print ("  xform jr", j, jr)
 152             vec2 = deepcopy(vec)
 153             for ci, (jl, jh) in enumerate(zip(j, jr)):
 154                 # exact same actual computation, just embedded in
 155                 # triple-nested for-loops
 156                 t1, t2 = vec[jl], vec[jh]
 157                 coeff = (math.cos((ci + 0.5) * math.pi / size) * 2.0)
 158                 vec2[jl] = t1 + t2
 159                 vec2[jl+halfsize] = (t1 - t2) * (1/coeff)
 160                 print ("coeff", size, i, k, "jl", jl, "jh", jh,
 161                        "i/n", (k+0.5)/size, coeff, vec[jl], vec[jh])
 162                 k += tablestep
 163             vec = vec2
 164         size //= 2
 165
 166     #vec.reverse()
 167     print("transform2 pre-itersum", vec)
 168     # Copy with bit-reversed permutation
 169
 170
 171     vec = transform_itersum(vec)
 172     print("transform2 intermediate", vec)
 173     return vec
 174
 175     size *= 2
 176     while size <= n:
 177         halfsize = size // 2
 178         tablestep = n // size
 179         ir = list(range(0, n, size))
 180         #ir.reverse()
 181         print ("itersum", size, ir)
 182         for i in ir:
 183             jr = list(range(i+halfsize, i + size-1))
 184             #jr.reverse()
 185             print ("itersum    jr", jr)
 186             for jh in jr:
 187                 # exact same actual computation, just embedded in
 188                 # triple-nested for-loops
 189                 #jh = i+halfsize-j
 190                 vec[jh] += vec[jh+1]
 191                 print ("    itersum", size, i, j, jh, jh+1)
 192         size *= 2
 193
 194     if reverse:
 195         vec = [vec[reverse_bits(i, levels)] for i in range(n)]
 196
 197     print("transform2 result", vec)
 198
 199     return vec
 200
 201
 202 # DCT type III, unscaled. Algorithm by Byeong Gi Lee, 1984.
 203 # See: https://www.nayuki.io/res/fast-discrete-cosine-transform-algorithms/lee-new-algo-discrete-cosine-transform.pdf
 204 def inverse_transform(vector, root=True, indent=0):
 205     idt = "   " * indent
 206     if root:
 207         vector = list(vector)
 208         vector[0] /= 2
 209     n = len(vector)
 210     if n == 1:
 211         return vector, [0]
 212     elif n == 0 or n % 2 != 0:
 213         raise ValueError()
 214     else:
 215         half = n // 2
 216         alpha = [vector[0]]
 217         beta  = [vector[1]]
 218         for i in range(2, n, 2):
 219             alpha.append(vector[i])
 220             beta.append(vector[i - 1] + vector[i + 1])
 221         print (idt, "n", n, "alpha 0", end=" ")
 222         for i in range(2, n, 2):
 223             print (i, end=" ")
 224         print ("beta 1", end=" ")
 225         for i in range(2, n, 2):
 226             print ("%d+%d" % (i-1, i+1), end=" ")
 227         print()
 228         inverse_transform(alpha, False, indent+1)
 229         inverse_transform(beta , False, indent+1)
 230         for i in range(half):
 231             x = alpha[i]
 232             y = beta[i] / (math.cos((i + 0.5) * math.pi / n) * 2)
 233             vector[i] = x + y
 234             vector[-(i + 1)] = x - y
 235             print (idt, " v[%d] = alpha[%d]+beta[%d]" % (i, i, i))
 236             print (idt, " v[%d] = alpha[%d]-beta[%d]" % (n-i-1, i, i))
 237         return vector
 238
 239
 240 def inverse_transform2(vector, root=True):
 241     n = len(vector)
 242     if root:
 243         vector = list(vector)
 244     if n == 1:
 245         return vector
 246     elif n == 0 or n % 2 != 0:
 247         raise ValueError()
 248     else:
 249         half = n // 2
 250         alpha = [0]
 251         beta  = [1]
 252         for i in range(2, n, 2):
 253             alpha.append(i)
 254             beta.append(("add", i - 1, i + 1))
 255         inverse_transform2(alpha, False)
 256         inverse_transform2(beta , False)
 257         for i in range(half):
 258             x = alpha[i]
 259             y = ("cos", beta[i], i)
 260             vector[i] = ("add", x, y)
 261             vector[-(i + 1)] = ("sub", x, y)
 262         return vector
 263
 264
 265 def failllll_transform2(block):
 266     N = len(block)
 267     cos = [0.0] * (N>>1)
 268
 269     front = deepcopy(block)
 270     back = deepcopy(block)
 271
 272     step = 1
 273     j = N *2
 274     half_N = N
 275     prev_half_N = N
 276
 277     while j > 1: #// Cycle of iterations Input Butterfly
 278         half_N = half_N >> 1
 279         current_PI_half_By_N = (math.pi / 2) / prev_half_N
 280         current_PI_By_N = 0.0
 281         step_Phase = current_PI_half_By_N * 2.0
 282         print ("n", N, "cos", end=" ")
 283         for i in range(half_N):
 284             #Precompute Cosine's coefficients
 285             a = current_PI_By_N + current_PI_half_By_N
 286             print (i, a / (math.pi), math.cos(a) * 2, end=" ")
 287             cos[i] = 0.5 / math.cos(a)
 288             current_PI_By_N += step_Phase
 289         print()
 290         k = 0
 291         for x in range(step):
 292             for i in range(half_N):
 293                 shift = k + prev_half_N - i - 1
 294                 back[k + i]          = front[k + i] + front[shift]
 295                 back[k + half_N + i] = (front[k + i] - front[shift]) * cos[i]
 296                 print ("xf coeff", N, j, i, x, "k/kh", k+i, k+half_N+i)
 297             k += prev_half_N
 298         temp = front
 299         front = back
 300         back = temp
 301         j = j >> 1
 302         step = step << 1
 303         prev_half_N = half_N
 304
 305     half_N = 2
 306     prev_half_N = 2
 307     j = 2
 308
 309     print("xform intermediate", front)
 310
 311     while j < N: # Cycle of Out ButterFly
 312         k = 0
 313         print ("out", j, N, step, half_N)
 314         for x in range(step):
 315             for i in range(half_N - 1):
 316                 back[k + (i << 1)] = front[k + i]
 317                 back[k + (i << 1) + 1] = (front[k + half_N + i] +
 318                                           front[k + half_N + i + 1])
 319                 print ("  out", j, x, i, "k", k,
 320                                 "k+i<<1", k+(i<<1), "hh1", k+half_N+i)
 321             back[k + ((half_N - 1) << 1)] = front[k + half_N - 1]
 322             back[k + (half_N << 1) - 1] = front[k + (half_N << 1) - 1]
 323             k += prev_half_N
 324
 325         temp = front
 326         front = back
 327         back = temp
 328         j = j << 1
 329         step = step >> 1
 330         half_N = prev_half_N
 331         prev_half_N = prev_half_N << 1
 332
 333     for i in range(N):
 334         block[i] = front[i] #// Unload DCT coefficients
 335     dN = 2.0
 336     #block[0] = block[0] / dN #// Compute DC.
 337
 338     print("transform2 result", block)
 339     return block
 340
 341
 342 if __name__ == '__main__':
 343     vector = range(8)
 344     ops = inverse_transform(vector)
 345     print (ops)