rendered paste body# encoding: utf-8"""by ulrik.sverdrup@gmail.comlic: GPL3"""from __future__ import divisionimport scipy as scimport scipy.linalgsz = 8base = sc.zeros([sz, sz, sz, sz])v = sc.array(xrange(sz))alpha = sc.ones([sz])* sc.sqrt(2/8)alpha[0] = sc.sqrt(1/8)# notice that sc.outer is the wrong one!outer = sc.multiply.outercos = sc.cos( outer(v, sc.pi/8*(v + 1/2)) )supercos = outer(cos,cos)superalpha = outer(alpha, alpha)for x in v: for y in v: base[x,y] = superalpha[x,y] * supercos[x,:,y,:]superalpha = outer(alpha,alpha)ncos = supercosprint superalphaprint supercosprint basem = sc.array([[52 , 55 , 61 , 66 , 70 , 61 , 64 , 73 ],[63 , 59 , 55 , 90 , 109 , 85 , 69 , 72],[62 , 59 , 68 , 113 , 144 , 104 , 66 , 73],[63 , 58 , 71 , 122 , 154 , 106 , 70 , 69],[67 , 61 , 68 , 104 , 126 , 88 , 68 , 70],[79 , 65 , 60 , 70 , 77 , 68 , 58 , 75],[85 , 71 , 64 , 59 , 55 , 61 , 65 , 83],[87 , 79 , 69 , 68 , 65 , 76 , 78 , 94]])m = m - 128print mdef dct(m): res = sc.zeros([sz,sz]) for x in v: for y in v: res[x,y] = sc.sum(sc.multiply(m,base[x,y])) return resprint sc.round_(dct(m))