论DCT和IDCT的重要性，汇编SIMD版第一，此贴第二，就是这么狂 :-)

输入: [1.000 2.000 3.000]
输出: [ 3.464 -1.414 0.000]
重建: [1.000 2.000 3.000]

[0] cos(0.0*π/3)*sqrt(1/N)*1.0 + cos(0.0*π/3)*sqrt(1/N)*2.0 + cos(0.0*π/3)*sqrt(1/N)*3.0 = 3.464
[1] cos(0.5*π/3)*sqrt(2/N)*1.0 + cos(1.5*π/3)*sqrt(2/N)*2.0 + cos(2.5*π/3)*sqrt(2/N)*3.0 = -1.414
[2] cos(1.0*π/3)*sqrt(2/N)*1.0 + cos(3.0*π/3)*sqrt(2/N)*2.0 + cos(5.0*π/3)*sqrt(2/N)*3.0 = 0.000

恭喜，你（基本）明白DCT了！可M为啥是正交矩阵？！

[0,0] cos(0.0*π/3)*sqrt(1/N)*cos(0.0*π/3)*sqrt(1/N) + cos(0.0*π/3)*sqrt(1/N)*cos(0.0*π/3)*sqrt(1/N) + cos(0.0*π/3)*sqrt(1/N)*cos(0.0*π/3)*sqrt(1/N)
[0,1] cos(0.0*π/3)*sqrt(1/N)*cos(0.5*π/3)*sqrt(2/N) + cos(0.0*π/3)*sqrt(1/N)*cos(1.5*π/3)*sqrt(2/N) + cos(0.0*π/3)*sqrt(1/N)*cos(2.5*π/3)*sqrt(2/N)
[0,2] cos(0.0*π/3)*sqrt(1/N)*cos(1.0*π/3)*sqrt(2/N) + cos(0.0*π/3)*sqrt(1/N)*cos(3.0*π/3)*sqrt(2/N) + cos(0.0*π/3)*sqrt(1/N)*cos(5.0*π/3)*sqrt(2/N)
[1,0] cos(0.5*π/3)*sqrt(2/N)*cos(0.0*π/3)*sqrt(1/N) + cos(1.5*π/3)*sqrt(2/N)*cos(0.0*π/3)*sqrt(1/N) + cos(2.5*π/3)*sqrt(2/N)*cos(0.0*π/3)*sqrt(1/N)
[1,1] cos(0.5*π/3)*sqrt(2/N)*cos(0.5*π/3)*sqrt(2/N) + cos(1.5*π/3)*sqrt(2/N)*cos(1.5*π/3)*sqrt(2/N) + cos(2.5*π/3)*sqrt(2/N)*cos(2.5*π/3)*sqrt(2/N)
[1,2] cos(0.5*π/3)*sqrt(2/N)*cos(1.0*π/3)*sqrt(2/N) + cos(1.5*π/3)*sqrt(2/N)*cos(3.0*π/3)*sqrt(2/N) + cos(2.5*π/3)*sqrt(2/N)*cos(5.0*π/3)*sqrt(2/N)
[2,0] cos(1.0*π/3)*sqrt(2/N)*cos(0.0*π/3)*sqrt(1/N) + cos(3.0*π/3)*sqrt(2/N)*cos(0.0*π/3)*sqrt(1/N) + cos(5.0*π/3)*sqrt(2/N)*cos(0.0*π/3)*sqrt(1/N)
[2,1] cos(1.0*π/3)*sqrt(2/N)*cos(0.5*π/3)*sqrt(2/N) + cos(3.0*π/3)*sqrt(2/N)*cos(1.5*π/3)*sqrt(2/N) + cos(5.0*π/3)*sqrt(2/N)*cos(2.5*π/3)*sqrt(2/N)
[2,2] cos(1.0*π/3)*sqrt(2/N)*cos(1.0*π/3)*sqrt(2/N) + cos(3.0*π/3)*sqrt(2/N)*cos(3.0*π/3)*sqrt(2/N) + cos(5.0*π/3)*sqrt(2/N)*cos(5.0*π/3)*sqrt(2/N)

如果用上复数（欧拉公式），就可以用等比数列求和公式了！！

from math import cos, sqrt, pi
import numpy as npπ = pi; N = 3
m = None # DCT矩阵
mt = None # IDCT矩阵，m的转置
mc = None # "C"版本DCT矩阵
ms = None # string版DCT矩阵

ne = lambda a, b: not np.allclose(a, b, atol=1E-6)def generate_dct_matrices():global m; global mt; global mc; global msm = np.zeros((N, N))# [[0]*3]*3不对。[0]*3得[0,0,0]，再*3复制引用‌而非创建独立副本。修改任意一行会影响所有行；debug 1小时mc = [[0] * N for _ in range(N)]ms = [[0] * N for _ in range(N)]for k in range(N):for n in range(N):mc[k][n] = m[k,n] = cos((n + 0.5) * k * pi / N)ms[k][n] = f'cos({(n+0.5)*k}*π/{N})'m[0,:] *= sqrt(1/N)for n in range(N):mc[0][n] *= sqrt(1/N)ms[0][n] += f'*sqrt(1/N)'m[1:,:] *= sqrt(2/N)for k in range(1,N):for n in range(N):mc[k][n] *= sqrt(2/N)ms[k][n] += f'*sqrt(2/N)'mt = m.Tif ne(m, mc): raise ValueError()# x是一维数组，shape=(3,)，np将其视为‌列向量‌(3,1)
# 3x3 x 3x1 = 3x1 再摆平
dct  = lambda x: np.dot(m,  x)
idct = lambda y: np.dot(mt, y)
# So, dct和idct其实是一个函数

np.set_printoptions(suppress=True,  # 禁用科学计数法precision=3,    # 保留3位小数floatmode='fixed'  # 固定小数位数
)generate_dct_matrices()x = np.array([1.0, 2.0, 3.0]); print("输入:", x)
y = dct(x); print("输出:", y)
r = idct(y); print("重建:", r)
ee = np.eye(m.shape[0])
if ne(np.matmul(m, mt), ee): raise ValueError()
print()y2 = [0] * N
for i in range(N):for j in range(N): y2[i] += mc[i][j] * x[j]
if ne(y, y2): raise ValueError()y2 = [0] * N
for i in range(N):s = ''for j in range(N): s += f'{ms[i][j]}*{x[j]} + 's = s[:-3]y2[i] = eval(s)print(f'{[i]} {s} = {y2[i]:.3f}')
if ne(y, y2): raise ValueError()
print('')print('恭喜，你（基本）明白DCT了！可M为啥是正交矩阵？！')a = ms; b = np.array(ms).T; e = [[0] * N for _ in range(N)]
for i in range(N):for j in range(N):s = ''for k in range(N): s += f'{a[i][k]}*{b[k][j]} + 's = s[:-3]print(f'[{i},{j}] {s}')e[i][j] = eval(s)
if ne(e, ee): raise ValueError()
print()print('如果用上复数（欧拉公式），就可以用等比数列求和公式了！！')