MATLAB从零开始实现离散小波变化DWT

一、基础目标

在MATLAB中从零开始实现离散小波变换（DWT）是一项能让你深入理解小波分析核心思想的实践。下面将梳理其基本原理、具体的实现步骤、关键代码以及一些重要的注意事项。

二、离散小波变换基础

离散小波变换（DWT）的核心思想是通过一组低通和高通滤波器，将信号分解成不同频率的子带。这种多分辨率分析方法让我们能够同时观察信号的时域和频域特性。

核心概念与Mallat算法

DWT的数学本质是通过滤波器组对信号进行多尺度分解。Mallat算法是实现这一过程的高效方法，它包含两个核心步骤：

1. 分解（正向DWT）：信号通过低通滤波器（获得近似系数cA）和高通滤波器（获得细节系数cD），然后进行下采样（通常为2倍）。
2. 重构（逆向DWT）：对近似系数和细节系数进行上采样，然后通过重构滤波器组，最终恢复原始信号。

下表对比了DWT分解中的两种主要系数：

小波滤波器系数获取

在实现DWT前，我们需要先获得小波滤波器的系数。MATLAB的wfilters函数可以帮我们获取标准小波的滤波器系数。

function[Lo_D,Hi_D,Lo_R,Hi_R]=get_wavelet_filters(wavelet_name)% 获取小波滤波器系数% 输入：wavelet_name - 小波名称，如'db4'% 输出：Lo_D - 分解低通滤波器, Hi_D - 分解高通滤波器% Lo_R - 重构低通滤波器, Hi_R - 重构高通滤波器% 检查输入的小波名称是否有效valid_wavelets={'haar','db1','db2','db3','db4','db5',...'sym2','sym3','sym4','coif1','coif2'};if~any(strcmpi(wavelet_name,valid_wavelets))error('不支持的小波名称。请使用如db4、sym2等标准小波。');end% 获取滤波器系数[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wavelet_name);fprintf('小波%s的滤波器系数：\n',wavelet_name);fprintf('低通分解滤波器 (Lo_D): [%s]\n',num2str(Lo_D,'%.4f '));fprintf('高通分解滤波器 (Hi_D): [%s]\n',num2str(Hi_D,'%.4f '));fprintf('低通重构滤波器 (Lo_R): [%s]\n',num2str(Lo_R,'%.4f '));fprintf('高通重构滤波器 (Hi_R): [%s]\n',num2str(Hi_R,'%.4f '));end

三、从零实现DWT分解与重构

单层DWT分解实现

下面是单层DWT分解的具体实现，重点关注卷积和下采样过程：

function[cA,cD]=my_dwt(x,Lo_D,Hi_D,mode)% 单层一维离散小波变换分解% 输入：x - 输入信号, Lo_D - 分解低通滤波器, Hi_D - 分解高通滤波器% mode - 边界延拓模式 ('sym'对称延拓, 'zpd'零延拓, 'ppd'周期延拓)% 输出：cA - 近似系数, cD - 细节系数ifnargin<4mode='sym';% 默认对称延拓endN=length(x);L=length(Lo_D);% 边界延拓处理switchmodecase'sym'% 对称延拓 (MATLAB默认方式)x_ext=[fliplr(x(1:L-1)),x,fliplr(x(end-L+2:end))];case'zpd'% 零延拓x_ext=[zeros(1,L-1),x,zeros(1,L-1)];case'ppd'% 周期延拓x_ext=[x(end-L+2:end),x,x(1:L-1)];otherwiseerror('不支持的边界延拓模式');end% 卷积运算（不使用MATLAB的conv函数，直接实现）cA_conv=zeros(1,N+L-1);cD_conv=zeros(1,N+L-1);fori=1:N+L-1forj=1:Lif(i-j+1)>=1&&(i