用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(N=8): """生成包含多种频率成分的测试信号""" 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 signal

2. DST-I到DST-IV的实现与解析

2.1 DST-I：最简单的反对称延拓

DST-I的变换核定义如下：

$$ X[k] = \sqrt{\frac{2}{N+1}} \sum_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(k+1)(n+1)}{N+1}\right) $$

Python实现：

def dst1(x): N = len(x) k = np.arange(1, N+1)[:, None] n = np.arange(1, N+1) basis = np.sin(np.pi * k * n / (N + 1)) return np.sqrt(2/(N+1)) * np.dot(basis, x)

可视化基函数：

def plot_dst_basis(dst_func, N=8): basis = dst_func(np.eye(N)) plt.figure(figsize=(10, 8)) for i in range(N): plt.subplot(N, 1, i+1) plt.stem(basis[i], use_line_collection=True) plt.ylabel(f'k={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_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(k+1)(2n+1)}{2N}\right) $$

其中$\eta_k$在$k=N-1$时为$1/\sqrt{2}$，否则为1。

实现代码：

def dst2(x): N = len(x) k = np.arange(1, N+1)[:, 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-III：DST-II的逆变换

DST-III实际上是DST-II的逆变换：

$$ X[k] = \sqrt{\frac{2}{N}} \eta_n \sum_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(2k+1)(n+1)}{2N}\right) $$

实现时注意$\eta_n$的处理：

def dst3(x): N = len(x) k = np.arange(N)[:, None] n = np.arange(1, N+1) 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_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(2k+1)(2n+1)}{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{2N+1}} \sum_{n=0}^{N-1} x[n] \sin\left(\frac{2\pi(k+1)(n+1)}{2N+1}\right) $$

实现代码：

def dst5(x): N = len(x) k = np.arange(1, N+1)[:, None] n = np.arange(1, N+1) 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{2N+1}} \sum_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(k+1)(2n+1)}{2N+1}\right) $$

Python实现：

def dst6(x): N = len(x) k = np.arange(1, N+1)[:, 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-VII：H.266/VVC的新宠

最新视频编码标准H.266/VVC采用了DST-VII：

$$ X[k] = \frac{2}{\sqrt{2N+1}} \sum_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(2k+1)(n+1)}{2N+1}\right) $$

实现代码：

def dst7(x): N = len(x) k = np.arange(N)[:, None] n = np.arange(1, N+1) 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_{n=0}^{N-1} x[n] \sin\left(\frac{\pi(2k+1)(2n+1)}{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(N=8): 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:(ki+1)*N, kj*N:(kj+1)*N] = np.outer(basis[ki], basis[kj]) plt.subplot(4, 2, i+1) plt.imshow(full_img, cmap='seismic', vmin=-1, vmax=1) plt.title(f'DST-{i+1} Basis Images') plt.colorbar() plt.tight_layout() plt.show()

4.2 能量压缩特性对比

不同DST变体对信号能量的压缩效率各异：

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-', label='Original') 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, label=f'DST-{i+1}') plt.legend() plt.title('Energy Compaction Comparison') plt.show()

通过这个对比可以明显看出，DST-VII和DST-II在能量压缩方面表现最为出色，这也是它们被视频编码标准采用的原因。