用Python和NumPy动手实现8种DST变换：从公式到可视化基图像

用Python和NumPy动手实现8种DST变换从公式到可视化基图像在信号处理领域离散正弦变换DST是一组与离散余弦变换DCT齐名的重要工具。不同于DCT的对称延拓特性DST通过反对称延拓方式处理信号特别适合解决某些特定边界条件的问题。本文将带您用Python和NumPy实现全部8种DST变换并通过可视化技术揭示其数学本质。1. DST变换基础与准备工作1.1 理解DST的核心概念离散正弦变换家族包含8种变体DST-I至DST-VIII每种变体对应不同的边界条件和延拓方式。与DCT不同DST采用反对称延拓这使得它在处理特定类型信号时具有独特优势反对称延拓信号在边界点被镜像反射并取反纯正弦基函数变换核仅由正弦函数构成能量压缩特性适合处理具有不连续边界条件的信号在开始编码前我们需要搭建Python环境import numpy as np import matplotlib.pyplot as plt from scipy.fft import dct, idct # 用于对比验证1.2 创建测试信号为了验证我们的实现先构造一个典型的测试信号def generate_test_signal(N8): 生成包含多种频率成分的测试信号 n np.arange(N) signal 0.5 * np.sin(2*np.pi*n/N) signal 0.3 * np.cos(4*np.pi*n/N) signal 0.2 * np.random.randn(N) return signal2. DST-I到DST-IV的实现与解析2.1 DST-I最简单的反对称延拓DST-I的变换核定义如下$$ X[k] \sqrt{\frac{2}{N1}} \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(k1)(n1)}{N1}\right) $$Python实现def dst1(x): N len(x) k np.arange(1, N1)[:, None] n np.arange(1, N1) basis np.sin(np.pi * k * n / (N 1)) return np.sqrt(2/(N1)) * np.dot(basis, x)可视化基函数def plot_dst_basis(dst_func, N8): basis dst_func(np.eye(N)) plt.figure(figsize(10, 8)) for i in range(N): plt.subplot(N, 1, i1) plt.stem(basis[i], use_line_collectionTrue) plt.ylabel(fk{i}) plt.suptitle(DST Basis Functions) plt.show()2.2 DST-II视频编码的宠儿DST-II在H.265/HEVC等视频编码标准中广泛应用$$ X[k] \sqrt{\frac{2}{N}} \eta_k \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(k1)(2n1)}{2N}\right) $$其中$\eta_k$在$kN-1$时为$1/\sqrt{2}$否则为1。实现代码def dst2(x): N len(x) k np.arange(1, N1)[:, None] n np.arange(N) basis np.sin(np.pi * k * (2*n 1) / (2*N)) eta np.ones(N) eta[-1] 1/np.sqrt(2) return np.sqrt(2/N) * eta * np.dot(basis, x)2.3 DST-IIIDST-II的逆变换DST-III实际上是DST-II的逆变换$$ X[k] \sqrt{\frac{2}{N}} \eta_n \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(2k1)(n1)}{2N}\right) $$实现时注意$\eta_n$的处理def dst3(x): N len(x) k np.arange(N)[:, None] n np.arange(1, N1) basis np.sin(np.pi * (2*k 1) * n / (2*N)) eta np.ones(N) eta[-1] 1/np.sqrt(2) return np.sqrt(2/N) * eta * np.dot(basis, x)2.4 DST-IV对称性最强的变体DST-IV具有优美的对称特性$$ X[k] \sqrt{\frac{2}{N}} \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(2k1)(2n1)}{4N}\right) $$Python实现def dst4(x): N len(x) k np.arange(N)[:, None] n np.arange(N) basis np.sin(np.pi * (2*k 1) * (2*n 1) / (4*N)) return np.sqrt(2/N) * np.dot(basis, x)3. DST-V到DST-VIII的高级实现3.1 DST-V周期加倍变换DST-V需要特别注意其归一化因子$$ X[k] \frac{2}{\sqrt{2N1}} \sum_{n0}^{N-1} x[n] \sin\left(\frac{2\pi(k1)(n1)}{2N1}\right) $$实现代码def dst5(x): N len(x) k np.arange(1, N1)[:, None] n np.arange(1, N1) basis np.sin(2 * np.pi * k * n / (2*N 1)) return 2/np.sqrt(2*N 1) * np.dot(basis, x)3.2 DST-VI半采样偏移变体DST-VI引入了半采样偏移$$ X[k] \frac{2}{\sqrt{2N1}} \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(k1)(2n1)}{2N1}\right) $$Python实现def dst6(x): N len(x) k np.arange(1, N1)[:, None] n np.arange(N) basis np.sin(np.pi * k * (2*n 1) / (2*N 1)) return 2/np.sqrt(2*N 1) * np.dot(basis, x)3.3 DST-VIIH.266/VVC的新宠最新视频编码标准H.266/VVC采用了DST-VII$$ X[k] \frac{2}{\sqrt{2N1}} \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(2k1)(n1)}{2N1}\right) $$实现代码def dst7(x): N len(x) k np.arange(N)[:, None] n np.arange(1, N1) basis np.sin(np.pi * (2*k 1) * n / (2*N 1)) return 2/np.sqrt(2*N 1) * np.dot(basis, x)3.4 DST-VIII最复杂的变体DST-VIII的归一化最为复杂$$ X[k] \frac{2}{\sqrt{2N-1}} \eta_k \eta_n \sum_{n0}^{N-1} x[n] \sin\left(\frac{\pi(2k1)(2n1)}{4N-2}\right) $$Python实现def dst8(x): N len(x) k np.arange(N)[:, None] n np.arange(N) basis np.sin(np.pi * (2*k 1) * (2*n 1) / (4*N - 2)) eta_k np.ones(N) eta_k[-1] 1/np.sqrt(2) eta_n np.ones(N) eta_n[-1] 1/np.sqrt(2) return (2/np.sqrt(2*N - 1)) * eta_k * (eta_n * np.dot(basis, x))4. 可视化分析与实际应用4.1 基图像生成技术基图像能直观展示变换的特性。以下代码生成所有DST变体的基图像def plot_all_dst_bases(N8): dst_types [dst1, dst2, dst3, dst4, dst5, dst6, dst7, dst8] plt.figure(figsize(15, 20)) for i, dst_func in enumerate(dst_types): basis dst_func(np.eye(N)) # 计算基图像外积 full_img np.zeros((N*N, N*N)) for ki in range(N): for kj in range(N): full_img[ki*N:(ki1)*N, kj*N:(kj1)*N] np.outer(basis[ki], basis[kj]) plt.subplot(4, 2, i1) plt.imshow(full_img, cmapseismic, vmin-1, vmax1) plt.title(fDST-{i1} Basis Images) plt.colorbar() plt.tight_layout() plt.show()4.2 能量压缩特性对比不同DST变体对信号能量的压缩效率各异变换类型能量集中效率适用场景DST-I中等边界值为0的信号DST-II高视频编码DST-VII极高新一代视频编码DST-VIII高特殊边界条件4.3 实际信号处理示例def compare_dst_performance(): x generate_test_signal(32) dst_types [dst1, dst2, dst3, dst4, dst5, dst6, dst7, dst8] plt.figure(figsize(12, 8)) plt.subplot(2, 1, 1) plt.plot(x, o-, labelOriginal) plt.subplot(2, 1, 2) for i, dst_func in enumerate(dst_types): coeffs dst_func(x) energy np.cumsum(coeffs**2) / np.sum(coeffs**2) plt.plot(energy, labelfDST-{i1}) plt.legend() plt.title(Energy Compaction Comparison) plt.show()通过这个对比可以明显看出DST-VII和DST-II在能量压缩方面表现最为出色这也是它们被视频编码标准采用的原因。