小波变换与时间序列分析

概述

小波变换（Wavelet Transform）是一种强大的时频分析方法，与傅里叶变换相比，它能够同时提供时间和频率信息，特别适合分析非平稳信号。在时间序列分析中，小波变换被广泛用于去噪、特征提取、多尺度分解等任务。¹

小波变换的核心思想是用一族小波基函数代替傅里叶变换中的正弦基函数，这些小波基函数具有有限的支撑（localized in time and frequency），使得小波变换能够同时捕获信号的局部特征和全局特征。

一、小波变换基础

1.1 从傅里叶变换到小波变换

傅里叶变换将信号分解为一组不同频率的正弦和余弦函数：

\hat{f} (ω) = \int_{- \infty}^{\infty} f (t) e^{- iω t} d t

问题：傅里叶变换只能告诉我们信号中包含哪些频率，无法告诉我们这些频率何时出现。

**短时傅里叶变换（STFT）**尝试通过加窗来获得时域信息：

STFT {f (t)} (u, ω) = \int_{- \infty}^{\infty} f (t) w (t - u) e^{- iω t} d t

问题：STFT 的窗口大小是固定的，无法同时获得好的时间分辨率和频率分辨率。

小波变换通过使用可伸缩的小波基函数解决了这个问题：

W_{f} (a, b) = \frac{1}{∣ a ∣} \int_{- \infty}^{\infty} f (t) ψ^{*} (\frac{t - b}{a}) d t

其中：

$ψ (t)$ 是母小波（Mother Wavelet）
$a$ 是尺度参数（Scale）
$b$ 是平移参数（Shift）

1.2 小波函数的条件

作为母小波， $ψ (t)$ 必须满足以下条件：

容许性条件：

C_{ψ} = \int_{0}^{\infty} \frac{∣ ψ ^ ( ω ) ∣ ^{2}}{ω} d ω < \infty

这意味着 $\hat{ψ} (0) = 0$ ，即小波函数必须具有零均值。

正则性条件：

小波函数应该足够光滑，以获得良好的频率局部化。

1.3 常用小波函数

小波	记号	特点	适用场景
Haar	haar	最早的小波，矩形波	边缘检测、基础分析
Daubechies	dbN	紧支撑、正则性好	通用信号处理
Symlet	symN	对称的 Daubechies	图像处理
Morlet	morl	余弦调制的高斯	连续小波分析
Mexican Hat	mexh	高斯的二阶导数	边缘检测
Meyer	meyr	频域定义	频谱分析

import pywt
import matplotlib.pyplot as plt
 
# 绘制不同小波函数
wavelets = ['haar', 'db4', 'sym4', 'morl', 'mexh']
fig, axes = plt.subplots(len(wavelets), 1, figsize=(10, 8))
 
for i, wavelet in enumerate(wavelets):
    wavelet_type, name = pywt.Wavelet(wavelet).wavename, pywt.Wavelet(wavelet).wavefun()
    psi, x = pywt.Wavelet(wavelet).wavefun()
    
    if len(psi) == 2:  # 小波函数
        axes[i].plot(x, psi[0], label=f'{wavelet}')
    else:
        axes[i].plot(psi, label=f'{wavelet}')
    axes[i].set_title(f'{wavelet} Wavelet')
    axes[i].grid(True)
 
plt.tight_layout()
plt.show()

二、连续小波变换（CWT）

2.1 定义

连续小波变换（CWT）对信号 $f (t)$ 进行连续尺度和位移的小波分解：

W_{f} (a, b) = \int_{- \infty}^{\infty} f (t) ψ_{a, b}^{*} (t) d t = \frac{1}{∣ a ∣} \int_{- \infty}^{\infty} f (t) ψ^{*} (\frac{t - b}{a}) d t

其中：

ψ_{a, b} (t) = \frac{1}{∣ a ∣} ψ (\frac{t - b}{a})

2.2 尺度和频率的关系

尺度和频率呈倒数关系：

大尺度 $a$ → 低频 → 粗分辨率 → 全局信息
小尺度 $a$ → 高频 → 细分辨率 → 局部信息

f_{a} = \frac{f _{c}}{a \cdot Δ t}

其中 $f_{c}$ 是小波的中心频率， $Δ t$ 是采样间隔。

2.3 Python 实现

import numpy as np
import pywt
import matplotlib.pyplot as plt
 
def continuous_wavelet_transform(signal, wavelet='morl', scales=None, 
                                 sampling_period=1.0):
    """
    连续小波变换
    
    Args:
        signal: 输入信号
        wavelet: 小波类型
        scales: 尺度和频率的对应关系
        sampling_period: 采样周期
    Returns:
        coeffs: 小波系数
        frequencies: 对应的频率
    """
    # 如果未指定 scales，根据信号自动生成
    if scales is None:
        # 典型范围：1 到 信号长度的一半
        scales = np.arange(1, len(signal) // 2)
    
    # 连续小波变换
    coefficients, frequencies = pywt.cwt(signal, scales, wavelet, 
                                         sampling_period=1.0)
    
    return coefficients, frequencies
 
 
def plot_scalogram(signal, wavelet='morl', cmap='jet'):
    """绘制小波尺度图（scalogram）"""
    coeffs, freqs = continuous_wavelet_transform(signal, wavelet)
    
    plt.figure(figsize=(12, 6))
    plt.imshow(np.abs(coeffs), extent=[0, len(signal), freqs[-1], freqs[0]],
               cmap=cmap, aspect='auto', interpolation='bilinear')
    plt.yscale('log')  # 对数尺度更好地显示频率
    plt.colorbar(label='Magnitude')
    plt.ylabel('Frequency')
    plt.xlabel('Time')
    plt.title(f'CWT Scalogram ({wavelet} wavelet)')
    plt.show()
 
 
# 示例：分析多频率信号
t = np.linspace(0, 1, 1000)
# 信号包含三个频率成分
signal = (np.sin(2 * np.pi * 5 * t) +   # 5 Hz
           0.5 * np.sin(2 * np.pi * 20 * t) +  # 20 Hz
           0.3 * np.sin(2 * np.pi * 50 * t))   # 50 Hz
 
plot_scalogram(signal, wavelet='morl')

三、离散小波变换（DWT）

3.1 定义

离散小波变换（DWT）使用离散的尺度和位移：

ψ_{j, k} (t) = 2^{j /2} ψ (2^{j} t - k)

DWT 的计算通过滤波器组实现，使用塔式算法（Pyramid Algorithm）。

3.2 塔式算法（Mallat 算法）

┌─────────────┐
│   信号 x    │
└──────┬──────┘
       │
       ├──────────→ conv(g) ──→ 近似系数 cA
       │
       └──────────→ conv(h) ──→ 细节系数 cD
                        │
                        └─→ 下采样(2) ──→ 下一层

其中：

$g$ 是低通滤波器（尺度函数）
$h$ 是高通滤波器（小波函数）

3.3 分解与重构

def wavelet_decompose(signal, wavelet='db4', level=3):
    """
    小波分解
    
    Args:
        signal: 输入信号
        wavelet: 小波类型
        level: 分解层数
    Returns:
        coeffs: 小波系数 [(cA_n, cD_n, ..., cD_1)]
    """
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    return coeffs
 
 
def wavelet_reconstruct(coeffs, wavelet='db4'):
    """小波重构"""
    return pywt.waverec(coeffs, wavelet)
 
 
# 示例分解
signal = np.random.randn(1024) + 0.5 * np.sin(2 * np.pi * 5 * np.linspace(0, 1, 1024))
coeffs = wavelet_decompose(signal, wavelet='db4', level=3)
 
# 分解结果
cA3, cD3, cD2, cD1 = coeffs
print(f"cA3 (低频近似): length = {len(cA3)}")  # 约 128
print(f"cD3 (第3层细节): length = {len(cD3)}")  # 约 128
print(f"cD2 (第2层细节): length = {len(cD2)}")  # 约 256
print(f"cD1 (第1层细节): length = {len(cD1)}")  # 约 512

3.4 多分辨率分析

小波分解提供了信号的多分辨率表示：

层    频率范围          系数
─────────────────────────────────────
cA3   [0, f_s/16)     低频近似
cD3   [f_s/16, f_s/8) 低频细节
cD2   [f_s/8, f_s/4)  中频细节
cD1   [f_s/4, f_s/2)  高频细节

四、时间序列去噪

4.1 小波阈值去噪

小波阈值去噪是最常用的去噪方法，假设噪声主要存在于高频系数中：

对信号进行小波分解
对细节系数应用阈值处理
重构去噪信号

def wavelet_denoise(signal, wavelet='db4', level=3, threshold_mode='soft'):
    """
    小波阈值去噪
    
    Args:
        signal: 含噪信号
        wavelet: 小波类型
        level: 分解层数
        threshold_mode: 'soft' (软阈值) 或 'hard' (硬阈值)
    Returns:
        denoised: 去噪后的信号
        threshold: 使用的阈值
    """
    # 小波分解
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    
    # 计算阈值（使用通用阈值估计）
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745  # MAD 估计
    threshold = sigma * np.sqrt(2 * np.log(len(signal)))
    
    # 对细节系数应用阈值
    denoised_coeffs = [coeffs[0]]  # 保留近似系数
    for detail_coeff in coeffs[1:]:
        # 软阈值: sign(x) * max(|x| - threshold, 0)
        if threshold_mode == 'soft':
            thresholded = pywt.threshold(detail_coeff, threshold, mode='soft')
        else:
            thresholded = pywt.threshold(detail_coeff, threshold, mode='hard')
        denoised_coeffs.append(thresholded)
    
    # 重构
    denoised = pywt.waverec(denoised_coeffs, wavelet)
    
    # 截断到原始长度
    return denoised[:len(signal)], threshold
 
 
def wavelet_denoise_visu(signal, wavelet='db4', level=3):
    """
    VisuShrink 阈值去噪
    使用更保守的阈值估计
    """
    # 小波分解
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    
    # 分层阈值（每层有不同阈值）
    thresholds = []
    for i, coeff in enumerate(coeffs):
        if i == 0:
            thresholds.append(0)  # 近似系数不处理
        else:
            sigma = np.median(np.abs(coeff)) / 0.6745
            t = sigma * np.sqrt(2 * np.log(len(coeff)))
            thresholds.append(t)
            coeffs[i] = pywt.threshold(coeff, t, mode='soft')
    
    denoised = pywt.waverec(coeffs, wavelet)
    return denoised[:len(signal)], thresholds

4.2 贝叶斯小波去噪

from scipy.stats import norm
 
def bayesian_wavelet_denoise(signal, wavelet='db4', level=3):
    """
    贝叶斯小波去噪
    
    假设：信号系数服从 Laplace 分布，噪声系数服从 Gaussian 分布
    """
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    
    # 估计噪声方差
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745
    
    # 分层处理
    denoised_coeffs = []
    for i, coeff in enumerate(coeffs):
        if i == 0:
            denoised_coeffs.append(coeff)
        else:
            # 贝叶斯估计
            # posterior mean of Laplace distribution under Gaussian noise
            noise_var = sigma ** 2
            signal_var = np.var(coeff)
            
            # posterior = signal_likelihood * signal_prior / evidence
            # 这里使用简化的 posterior mean
            shrinkage = noise_var / (noise_var + min(signal_var, np.var(coeff)))
            denoised_coeff = coeff * (1 - shrinkage)
            denoised_coeffs.append(denoised_coeff)
    
    denoised = pywt.waverec(denoised_coeffs, wavelet)
    return denoised[:len(signal)]

五、特征提取与多尺度分析

5.1 多尺度能量特征

def extract_multiscale_energy(signal, wavelet='db4', max_level=5):
    """
    提取多尺度能量特征
    
    Returns:
        energies: 各尺度的能量
        ratios: 相邻尺度的能量比
    """
    coeffs = pywt.wavedec(signal, wavelet, level=max_level)
    
    energies = [np.sum(c**2) for c in coeffs]
    ratios = [energies[i] / (energies[i+1] + 1e-10) for i in range(len(energies)-1)]
    
    return energies, ratios
 
 
def extract_multiscale_statistics(signal, wavelet='db4', level=5):
    """
    提取多尺度统计特征
    """
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    
    features = {}
    
    # 近似系数统计
    features['approx_mean'] = np.mean(coeffs[0])
    features['approx_std'] = np.std(coeffs[0])
    features['approx_energy'] = np.sum(coeffs[0]**2)
    
    # 细节系数统计
    for i, c in enumerate(coeffs[1:], 1):
        features[f'detail_{i}_mean'] = np.mean(c)
        features[f'detail_{i}_std'] = np.std(c)
        features[f'detail_{i}_energy'] = np.sum(c**2)
        features[f'detail_{i}_entropy'] = -np.sum(c**2 * np.log(c**2 + 1e-10))
    
    return features

5.2 边缘检测与突变点识别

def detect_change_points(signal, wavelet='db4', threshold_factor=3.0):
    """
    使用小波变换检测信号中的突变点
    """
    # 使用一阶导数小波（Haar）检测突变
    coeffs = pywt.wavedec(signal, wavelet, level=3)
    
    # 最后几层细节系数的模极大值
    detail_coeffs = np.abs(coeffs[-1])
    
    # 计算阈值
    threshold = threshold_factor * np.std(detail_coeffs)
    
    # 找极大值点
    change_points = np.where(detail_coeffs > threshold)[0]
    
    return change_points
 
 
def wavelet_edge_detection(signal, wavelet='haar'):
    """
    使用小波进行边缘检测
    
    等价于信号的一阶导数
    """
    # 一层分解
    cA, cD = pywt.dwt(signal, wavelet)
    
    # cD 是信号的近似导数
    # 边缘位置在 cD 的极值点
    
    edges = np.where(np.abs(np.diff(cD)) > np.std(cD))[0]
    
    return edges

六、小波神经网络

6.1 可学习小波层

class LearnableWaveletLayer(nn.Module):
    """
    可学习的小波变换层
    
    替代固定的小波基函数，
    端到端学习最适合任务的小波
    """
    def __init__(self, in_channels, out_channels, kernel_size=31):
        super().__init__()
        
        self.kernel_size = kernel_size
        
        # 可学习的低通滤波器 (尺度函数)
        self.low_filter = nn.Parameter(
            torch.randn(in_channels, out_channels, kernel_size) * 0.01
        )
        
        # 可学习的高通滤波器 (小波函数)
        self.high_filter = nn.Parameter(
            torch.randn(in_channels, out_channels, kernel_size) * 0.01
        )
        
        # 初始化为近似 Haar 小波
        self._initialize_haar()
    
    def _initialize_haar(self):
        """初始化为 Haar 小波"""
        with torch.no_grad():
            # Haar 低通: [0.5, 0.5]
            self.low_filter.copy_(torch.tensor([[[0.5, 0.5]]] * self.low_filter.shape[0] 
                                               * self.low_filter.shape[1]).reshape_as(self.low_filter))
            # Haar 高通: [0.5, -0.5]
            self.high_filter.copy_(torch.tensor([[[0.5, -0.5]]] * self.high_filter.shape[0] 
                                                * self.high_filter.shape[1]).reshape_as(self.high_filter))
    
    def forward(self, x):
        """
        Args:
            x: [batch, channels, length]
        Returns:
            low: 低频成分 [batch, channels, length//2]
            high: 高频成分 [batch, channels, length//2]
        """
        # 1D 卷积实现小波变换
        low = F.conv1d(x, self.low_filter, padding=self.kernel_size//2, stride=2)
        high = F.conv1d(x, self.high_filter, padding=self.kernel_size//2, stride=2)
        
        return low, high
 
 
class InverseLearnableWaveletLayer(nn.Module):
    """逆小波变换层"""
    def __init__(self, in_channels, out_channels, kernel_size=31):
        super().__init__()
        
        self.kernel_size = kernel_size
        
        # 可学习的重构滤波器
        self.low_filter = nn.Parameter(torch.randn(out_channels, in_channels, kernel_size) * 0.01)
        self.high_filter = nn.Parameter(torch.randn(out_channels, in_channels, kernel_size) * 0.01)
    
    def forward(self, low, high):
        """
        Args:
            low: 低频成分 [batch, channels, length//2]
            high: 高频成分 [batch, channels, length//2]
        Returns:
            x: 重构信号 [batch, channels, length]
        """
        batch, channels, half_len = low.shape
        
        # 上采样并卷积
        low_up = F.interpolate(low, scale_factor=2, mode='linear', align_corners=True)
        high_up = F.interpolate(high, scale_factor=2, mode='linear', align_corners=True)
        
        low_conv = F.conv1d(low_up, self.low_filter, padding=self.kernel_size//2)
        high_conv = F.conv1d(high_up, self.high_filter, padding=self.kernel_size//2)
        
        return low_conv + high_conv

6.2 WTConv 层

WTConv 是 2024 年 ECCV 论文提出的卷积层，使用小波变换实现多频率响应。²

class WTConv(nn.Module):
    """
    Wavelet Transform Convolution
    
    核心思想：用小波变换替代池化操作，
    同时捕获多个频率的信息
    """
    def __init__(self, in_channels, out_channels, kernel_size=3, 
                 wavelet='db1', levels=2):
        super().__init__()
        
        self.wavelet = wavelet
        self.levels = levels
        
        # 小波分解
        self.decompose = pywt.Wavelet(wavelet)
        
        # 卷积层
        self.conv = nn.ModuleList()
        for _ in range(levels + 1):  # 多频率分支
            self.conv.append(
                nn.Sequential(
                    nn.Conv1d(in_channels, out_channels, kernel_size, 
                             padding=kernel_size//2),
                    nn.BatchNorm1d(out_channels),
                    nn.GELU()
                )
            )
    
    def forward(self, x):
        """
        Args:
            x: [batch, channels, length]
        Returns:
            out: [batch, out_channels, length]
        """
        # 分解为多个频率
        coeffs = pywt.wavedec(x.cpu().numpy(), self.wavelet, level=self.levels)
        coeffs_cpu = coeffs
        
        outputs = []
        for i, c in enumerate(coeffs_cpu):
            # 转换回 tensor
            c_tensor = torch.tensor(c, dtype=x.dtype, device=x.device)
            # 插值到原始长度
            if c_tensor.shape[-1] != x.shape[-1]:
                c_tensor = F.interpolate(
                    c_tensor.unsqueeze(0), 
                    size=x.shape[-1], 
                    mode='linear', 
                    align_corners=True
                ).squeeze(0)
            # 应用卷积
            out = self.conv[i](c_tensor.unsqueeze(0).transpose(1, 2)).transpose(1, 2)
            outputs.append(out)
        
        # 融合多频率输出
        out = torch.stack(outputs, dim=-1).mean(dim=-1)
        
        return out

6.3 完整的小波增强预测模型

class WaveletTimeSeriesModel(nn.Module):
    """
    小波增强的时间序列预测模型
    
    结合小波变换的多尺度特征提取能力和深度学习的预测能力
    """
    def __init__(self, input_dim, hidden_dim, output_dim, 
                 wavelet='db4', n_levels=3):
        super().__init__()
        
        # 小波分解层
        self.wavelet = wavelet
        
        # 编码器
        self.encoder = nn.ModuleList()
        for _ in range(n_levels + 1):  # 多尺度分支
            self.encoder.append(
                nn.Sequential(
                    nn.Linear(input_dim, hidden_dim),
                    nn.GELU(),
                    nn.Linear(hidden_dim, hidden_dim)
                )
            )
        
        # 融合层
        self.fusion = nn.MultiheadAttention(
            hidden_dim, num_heads=4, batch_first=True
        )
        
        # 预测头
        self.decoder = nn.Sequential(
            nn.Linear(hidden_dim, hidden_dim),
            nn.GELU(),
            nn.Linear(hidden_dim, output_dim)
        )
    
    def forward(self, x):
        """
        x: [batch, seq_len, input_dim]
        """
        # 转换为 [batch, input_dim, seq_len] 用于小波变换
        x_wav = x.transpose(1, 2)
        
        # 小波分解
        coeffs = pywt.wavedec(
            x_wav.cpu().numpy(), 
            self.wavelet, 
            level=len(self.encoder) - 1
        )
        
        # 各尺度编码
        scale_outputs = []
        for i, c in enumerate(coeffs):
            c_tensor = torch.tensor(c, dtype=x.dtype, device=x.device)
            # 插值到统一长度
            c_tensor = F.interpolate(
                c_tensor.unsqueeze(0), 
                size=x.shape[1], 
                mode='linear', 
                align_corners=True
            ).squeeze(0).transpose(0, 1)  # [input_dim, seq_len]
            
            encoded = self.encoder[i](c_tensor.transpose(0, 1))  # [seq_len, hidden_dim]
            scale_outputs.append(encoded)
        
        # 注意力融合
        fused, _ = self.fusion(scale_outputs[0], scale_outputs, scale_outputs)
        
        # 预测
        output = self.decoder(fused)  # [seq_len, output_dim]
        
        return output

七、实践指南

7.1 小波选择建议

信号类型	推荐小波	原因
突变信号	Haar	矩形波，适合检测边缘
光滑信号	Db4, Db6	正则性好，频率局部化强
对称要求	Sym4, Sym6	近对称，减少相位失真
振荡信号	Morlet	天然振荡形式
图像/2D	Db8, Coif4	常用的标准选择

7.2 分解层数选择

def optimal_wavelet_level(signal_length, min_coef_length=20):
    """
    自动确定最优分解层数
    
    约束：最细尺度的系数长度至少为 min_coef_length
    """
    import math
    level = int(math.log2(signal_length / min_coef_length))
    return max(1, level)

7.3 去噪参数调优

def tune_denoise_parameters(signal_noisy, signal_clean, 
                           wavelets=['db4', 'db6', 'sym4'],
                           levels=[2, 3, 4, 5],
                           modes=['soft', 'hard']):
    """
    通过网格搜索找到最优去噪参数
    """
    best_mse = float('inf')
    best_params = None
    
    for wavelet in wavelets:
        for level in levels:
            for mode in modes:
                try:
                    denoised, _ = wavelet_denoise(
                        signal_noisy, wavelet, level, mode
                    )
                    mse = np.mean((denoised - signal_clean)**2)
                    
                    if mse < best_mse:
                        best_mse = mse
                        best_params = {'wavelet': wavelet, 
                                      'level': level, 
                                      'mode': mode}
                except:
                    continue
    
    return best_params, best_mse

Metaphor

探索

小波变换与时间序列分析

概述

一、小波变换基础

1.1 从傅里叶变换到小波变换

1.2 小波函数的条件

1.3 常用小波函数

二、连续小波变换（CWT）

2.1 定义

2.2 尺度和频率的关系

2.3 Python 实现

三、离散小波变换（DWT）

3.1 定义

3.2 塔式算法（Mallat 算法）

3.3 分解与重构

3.4 多分辨率分析

四、时间序列去噪

4.1 小波阈值去噪

4.2 贝叶斯小波去噪

五、特征提取与多尺度分析

5.1 多尺度能量特征

5.2 边缘检测与突变点识别

六、小波神经网络

6.1 可学习小波层

6.2 WTConv 层

6.3 完整的小波增强预测模型

七、实践指南

7.1 小波选择建议

7.2 分解层数选择

7.3 去噪参数调优

八、参考

相关阅读

关系图谱

目录

反向链接

Metaphor

探索

小波变换与时间序列分析

概述

一、小波变换基础

1.1 从傅里叶变换到小波变换

1.2 小波函数的条件

1.3 常用小波函数

二、连续小波变换（CWT）

2.1 定义

2.2 尺度和频率的关系

2.3 Python 实现

三、离散小波变换（DWT）

3.1 定义

3.2 塔式算法（Mallat 算法）

3.3 分解与重构

3.4 多分辨率分析

四、时间序列去噪

4.1 小波阈值去噪

4.2 贝叶斯小波去噪

五、特征提取与多尺度分析

5.1 多尺度能量特征

5.2 边缘检测与突变点识别

六、小波神经网络

6.1 可学习小波层

6.2 WTConv 层

6.3 完整的小波增强预测模型

七、实践指南

7.1 小波选择建议

7.2 分解层数选择

7.3 去噪参数调优

八、参考

相关阅读

Footnotes

关系图谱

目录

反向链接