论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)
如果用上复数（欧拉公式），就可以用等比数列求和公式了！！

1. from math import cos, sqrt, pi
2. import numpy as np
4. π = pi; N = 3
5. m = None # DCT矩阵
6. mt = None # IDCT矩阵，m的转置
7. mc = None # "C"版本DCT矩阵
8. ms = None # string版DCT矩阵
10. ne = lambda a, b: not np.allclose(a, b, atol=1E-6)
12. def generate_dct_matrices():
13.   global m; global mt; global mc; global ms
14.   m = np.zeros((N, N))
15.   # [[0]*3]*3不对。[0]*3得[0,0,0]，再*3复制引用‌而非创建独立副本。修改任意一行会影响所有行；debug 1小时
16.   mc = [[0] * N for _ in range(N)]
17.   ms = [[0] * N for _ in range(N)]
18.   for k in range(N):
19.     for n in range(N):
20.       mc[k][n] = m[k,n] = cos((n + 0.5) * k * pi / N)
21.       ms[k][n] = f'cos({(n+0.5)*k}*π/{N})'
22.   m[0,:] *= sqrt(1/N)
23.   for n in range(N):
24.     mc[0][n] *= sqrt(1/N)
25.     ms[0][n] += f'*sqrt(1/N)'
26.   m[1:,:] *= sqrt(2/N)
27.   for k in range(1,N):
28.     for n in range(N):
29.       mc[k][n] *= sqrt(2/N)
30.       ms[k][n] += f'*sqrt(2/N)'
31.   mt = m.T
32.   if ne(m, mc): raise ValueError()
34. # x是一维数组，shape=(3,)，np将其视为‌列向量‌(3,1)
35. # 3x3 x 3x1 = 3x1 再摆平
36. dct  = lambda x: np.dot(m,  x)
37. idct = lambda y: np.dot(mt, y)
38. # So, dct和idct其实是一个函数
40. np.set_printoptions(
41.   suppress=True,  # 禁用科学计数法
42.   precision=3,    # 保留3位小数
43.   floatmode='fixed'  # 固定小数位数
44. )
46. generate_dct_matrices()
48. x = np.array([1.0, 2.0, 3.0]); print("输入:", x)
49. y = dct(x); print("输出:", y)
50. r = idct(y); print("重建:", r)
51. ee = np.eye(m.shape[0])
52. if ne(np.matmul(m, mt), ee): raise ValueError()
53. print()
55. y2 = [0] * N
56. for i in range(N):
57.   for j in range(N): y2[i] += mc[i][j] * x[j]
58. if ne(y, y2): raise ValueError()
60. y2 = [0] * N
61. for i in range(N):
62.   s = ''
63.   for j in range(N): s += f'{ms[i][j]}*{x[j]} + '
64.   s = s[:-3]
65.   y2[i] = eval(s)
66.   print(f'{[i]} {s} = {y2[i]:.3f}')
67. if ne(y, y2): raise ValueError()
68. print('')
70. print('恭喜，你（基本）明白DCT了！可M为啥是正交矩阵？！')
72. a = ms; b = np.array(ms).T; e = [[0] * N for _ in range(N)]
73. for i in range(N):
74.   for j in range(N):
75.     s = ''
76.     for k in range(N): s += f'{a[i][k]}*{b[k][j]} + '
77.     s = s[:-3]
78.     print(f'[{i},{j}] {s}')
79.     e[i][j] = eval(s)
80. if ne(e, ee): raise ValueError()
81. print()
83. print('如果用上复数（欧拉公式），就可以用等比数列求和公式了！！')

复制代码

 

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！