用符号运算验证DCT的矩阵是正交矩阵

输入: [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 ms
  m = 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.T
  if ne(m, mc): raise ValueError()

# x是一维数组，shape=(3,)，np将其视为‌列向量‌(3,1)
# 3x3 x 3x1 = 3x1. 720p是1280x720
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('如果用上复数（欧拉公式），就可以用等比数列求和公式了！！')

cos²(θ/3) + cos²(θ) + cos²(5θ/3) = ? sympy不能化简。DeekSeek可以。

'''
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 pi as π, cos
import numpy as np
α = π / 3
N = 3
print(f'{(2/N)*(cos(α)**2 + cos(3*α)**2 + cos(5*α)**2):.3f}')

δ_plural = 0; βs = np.linspace(0, 2*π, 100);
for β in βs:
  a = np.cos(β) + np.cos(3*β) + np.cos(5*β)
  b = np.cos(3*β) * (2*np.cos(2*β) + 1)
  δ = abs(a-b); δ_plural += δ
  if δ > 1E-6: raise ValueError()
  print(f'{δ_plural / βs.size:.3f}')

# cos²γ = (1 + cos2γ) / 2 = 1/2 + cos2γ/2
β = 2 * α
print(f'{(cos(3*β) * (2*cos(2*β) + 1) + 3/2) * 2/N:.3f}')

噫，算出1可不易。

实部永远是实部，虚部永远是虚部，实对实来虚对虚。

我还是不认为AI有啥智能。它提到了等比数列求和，对我是个小小的打击（我还以为是我首创），无非是知道的人不少，而且发到网上了，AI看到了，我没看到而已。至于上学时么，可能老师没讲，更可能我没去上课或没仔细听，但我可以肯定我看到的课本上没讲。

AI还说：

① 设计正交矩阵时可利用三角函数的性质。矩阵的列向量（或行向量）必须两两正交且长度（模）为1。
② 最简单的正交矩阵是二维或三维的旋转矩阵（吗？）窃以为是单位矩阵。
③ 哈达玛矩阵（Hadamard Matrix）是由1和-1构成的正交矩阵。

在H.264/AVC等编码标准中，4×4和8×8 Hadamard变换被用于SATD计算流程：
1. 对残差信号矩阵执行二维变换，生成频域系数
2. 计算系数绝对值之和作为SATD值
3. 通过归一化处理消除阶数差异影响
作为DCT的替代方案，Hadamard变换在JPEG标准中实现能量集中效率达87%，较DCT低约9%，但运算速度提升显著。
微软BitNet v2框架用Hadamard变换优化LLM激活值量化流程。

〔为啥傅里叶变换在解密中也有用？〕〔矩阵动画〕