小波分析-10利用滤波器组实现离散小波变换的matlab实现

现代信号处理

(ModernSignalProcessing)信号处理课程组信息与控制工程学院第七章小波变换主讲教师：张晓光利用滤波器组实现离散小波变换的MATLAB实现主讲教师：张晓光利用滤波器组实现离散小波变换的MATLAB实现一、一维信号的小波分解在MATLAB中，dwt函数可以用于计算一维信号的单级离散小波变换，即利用滤波器组实现离散小波变换的分解算法。该函数的基本语法为

[C,L]=dwt(X,'WaveletName',level);X是输入信号，通常是一个向量或矩阵。'WaveletName'是所使用的小波的名称，例如'db1'、'haar’等level是分解的层数，表示将信号分解为多少个子带。C是近似系数，表示信号的主要趋势和低频成分。L是细节系数，表示信号的高频成分。同时，当分解层数level=1时，即实现一级DWT，其框图如图所示；当分解层数level=3时，即实现三级DWT，其框图如图所示。例7-11试利用信号离散小波变换函数dwt对“doorbell.wav”一维语音信号进行小波分析。%读取音频文件'doorbell.wav’

[W,fs]=audioread('doorbell.wav');%绘制原始音频信号的时域波形W=W(:,1);figure;plot(W(:,1))%使用Coiflet4小波，对单声道音频信号进行一级离散小波分解dwtmode('sym‘)；[cA_L,cD_L]=dwt(W(:,1),'coif4');l=length(cA_L);%画图figure;subplot(211);plot(cA_L);title(‘近似信号’);%绘制低频近似信号subplot(212);plot(cD_L);title('细节信号');%绘制高频细节信号图7.6.8“doorbell”信号的一级DWT例7-12利用MATLAB和Daub10小波对chirp信号进行一级DWT分解，中间在样点250-270之间插入一段毛刺干扰。MATLAB提供了chirp函数，函数用法如下s=chirp(t,f0,t1,f1)其中，f0是在t=0时刻的起始频率，f1是在t=t1时刻的终止频率。MATLAB程序：fs=1000;T=1/fs;%设置采样频率fs=1000;计算采样周期t=0:T:1023*T%生成时间向量，总时长为1023*T，共1024个采样点s=chirp(t,0,1,100);%生成线性调频信号：起始频率0Hz，1秒时达到100Hzs(250:270)=1;%在第250至270采样点人为加入幅值为1的干扰（毛刺）subplot(311);plot(t/T,s);%绘制原始含噪信号的时域波形wname=['db10'];%定义小波基类型为Daubechies10阶[cA,cD]=dwt(s,wname);%对信号进行单层离散小波分解subplot(312);plot(cA);%绘制低频近似分量（包含主要趋势信息）subplot(313);plot(cD);%绘制高频细节分量（突出突变点特征）图7.6.9含毛刺正弦信号的一级DWT例7-13利用MATLAB提供的wavedec函数对例7-12题中的chirp信号进行三级离散小波分解，并绘出各级系数的波形。解:我们可以逐级利用dwt来进行计算，但函数wavede可以一次完成多级分解。函数wavedec如前所述将产生一个输出向量C，其包含了各级的输出系数。它也包含了一个标签向量L，以标示向量C中各个系数的长度。MATLAB程序：%设置小波基类型为Daubechies10阶；

wavef='db10';%设置小波分解层数为3层level=3;%定义采样参数fs=1000;T=1/fs;t=0:T:1023*T%生成线性调频信号（起始频率0Hz，1秒时达100Hz）；s=chirp(t,0,1,100);%在第250~270采样点人为叠加幅值为1的毛刺干扰s(250:270)=1;%绘制原始含噪信号的时域波形plot(t/T,s);%对信号进行多尺度小波分解[CL]=wavedec(s,level,wavef);%C:小波系数，L:分解层信息%创建新图形窗口并清空figure(1);clf;strA='近似系数level';%低频分量标题前缀strD='细节系数level';%高频分量标题前缀n=1024;end%遍历每一层分解结果fori=1:level%提取第i层近似系数（低频成分）

A=appcoef(C,L,wavef,i);subplot(3,2,2*i-1);plot(A);%绘制近似系数title([strAnum2str(i)]);%提取第i层细节系数（高频成分）D=detcoef(C,L,i);subplot(3,2,2*i);plot(D);%绘制细节系数title([strDnum2str(i)]);n=n/2;图7.6.10含毛刺正弦信号的三级DWT2.一维信号的小波重构MATLAB提供了idwt函数用于一维信号单级离散小波变换的信号重建

X=idwt(C,'wname')C：表示输入的小波系数，通常是一个包含近似系数和细节系数的向量或矩阵；'wname'为指定的小波函数名称，例如'db1'、'coif1‘等；返回值X：表示输出的重构信号。同时，MATLAB提供了waverec函数用于一维信号多级离散小波变换的信号重建

y=waverec(c,wname)C：表示一个包含小波系数的向量或矩阵；返回值y：表示输出的重构信号。重构算法也存在与分解算法类似的多级树形结构，当level=1时，即实现一级IDWT，其框图如图所示；当level=2时，即实现两级IDWT，其框图如图所示。例7-14考虑下面给定的信号对其0~4.095s的一段信号按频率1000Hz进行采样，再叠加一个频率为400Hz的正弦噪声信号。利用Daub4小波对此含有噪声的信号进行4级DWT，然后重建信号。MATLAB程序：%定义多项式系数向量xx=[1-12.453.61-95.05461.814-8.97];%设置采样参数fs=1000;%采样频率1000HzT=1/fs;%采样周期t=0:T:4095*T;%生成时间向量，总长度4096个采样点%初始化多项式ff=0;N=length(x);%多项式阶数+1（此处为6项，对应5次多项式）fori=1:N%累加各项构成多项式f(t)=x(1)*t^5+x(2)*t^4+...+x(6)*t^0f=f+x(i)*t.^(N-i);end%归一化处理并将值映射到[0,1]区间mx=max(f);%最大值mn=min(f);%最小值f=(f-mn)/(mx-mn);%归一化公式%添加400Hz高频噪声（幅度为0.3）fn=f+0.3*sin(2*pi*400*t);%设置小波参数并进行多级分解wavef='db10';%使用Daubechies10小波s=4;%分解层数[C,L]=wavedec(fn,s,wavef);%对含噪信号进行小波分解[C1,L1]=wavedec(f,s,wavef);%对无噪声信号进行小波分解%提取各层系数并绘图strA='近似系数level';%低频分量标题前缀strD='细节系数level';%高频分量标题前缀fori=1:s%提取第i层近似系数（低频成分）

A=appcoef(C,L,wavef,i);subplot(s,2,2*i-1);plot(A);%绘制近似系数

title([strAnum2str(i)]);

%提取第i层细节系数（高频成分）

D=detcoef(C,L,i);subplot(s,2,2*i);plot(D);%绘制细节系数

title([strDnum2str(i)]);end%信号重构fr=waverec(C,L,wavef);%重构含噪声信号fr1=waverec(C1,L1,wavef);%重构无噪声信号%重构结果绘图subplot(1,2,1);plot(t,fr);title('重建含噪声信号');subplot(1,2,2);plot(t,fr1);title('重建不含噪声信号');小波系数如图7.6.11所示，各级的近似系数在左边，细节系数在右边，第一级细节系数的数值最大因为这级对应信号中的噪声。含噪声信号和不含噪声信号的重建效果如图7.6.12所示。在上述代码中，D表示细节系数，一般d1表示第一级细节系数向量，d2表示第二级细节系数向量，以此类推。如果我们分别将d1、d2以及他们两个都设置为零后重建信号，给出设置第一级和第二级细节系数d1和d2为零后重建信号的MATLAB代码，其他类似，重建信号如图7.6.13所示%对含噪声信号进行4层小波分解，得到分解系数C和分解层信息L[C,L]=wavedec(fn,4,wavef);%提取第1至4层的细节系数（d1~d4）[d1,d2,d3,d4]=detcoef(C,L,1:4);%提取第4层的近似系数（a4）a4=appcoef(C,L,wavef,4);%重组系数：将近似系数a4放在首位，随后按d4到d1逆序排列细节系数C=[a4,d4,d3,d2,d1];%将第1层和第2层细节系数清零（去除高频噪声）d1=d1*0;d2=d1*0;%使用修改后的系数重构信号（抑制前两层高频成分）fr3=waverec(C,L,wavef);%绘制重构后的信号波形plot(t,fr3);title('d1和d2置零后重建的含噪声信号');3.二维信号的小波分解和重构MATLAB提供了dwt2和idwt2函数，用于计算二维单级离散小波变换的信号分解和重建。dwt2函数的基本语法为[LL,LH,HL,HH]=dwt2(X,'wname’)X：表示输入的二维信号，通常是一个矩阵，代表图像的灰度值；'wname'：表示指定的小波函数名称，例如'db1'、'coif1’等；LL：表示输出的近似系数，表示图像的低频成分；LH：表示输出的水平细节系数，表示图像在垂直方向的高频成分；HL：表示输出的垂直细节系数，表示图像在水平方向的高频成分；HH：表示输出的对角细节系数，表示图像在对角方向的高频成分。idwt2函数的基本语法为X=idwt2(C,S,'wname’)C：表示输入的小波系数，通常是一个包含近似系数和细节系数的向量或矩阵；S：表示输入的系数对应的尺寸；'wname'：表示指定的小波函数名称；返回值X：表示输出的重构信号或图像。同时，MATLAB还提供了wavedec2和waverec2函数用于计算二维多级离散小波变换的信号分解和重建。wavedec2函数的基本语法为[C,S]=wavedec2(X,N,'wname')X：表示输入的二维信号，通常是一个矩阵，代表图像的灰度值；N：表示分解的层数，即小波变换的级别数；'wname'：表示指定的小波函数的名称；C：表示输出的小波系数；S：表示输出的系数对应的尺寸。waverec2函数的基本语法为X=waverec2(C,S,'wname')C：表示输入的小波系数；S：表示输入的系数对应的尺寸；'wname'：表示指定的小波函数的名称；返回值X：表示输出的重构信号或图像。例7-15用dwt2和idwt2函数，对二维图像“玛丽莲∙爱因斯坦”进行单级离散小波分解与重构。MATLAB程序如下:%定义图像路径imgPath='C:/tupianmatlab/ma.jpg';%读取原始图像并转换为灰度图（如果是彩色图）x=imread(imgPath);ifsize(x,3)==3Y=rgb2gray(x);%RGB转灰度，保留亮度信息elseY=x;%已是灰度图则直接使用endY=im2double(Y);%转换为双精度浮点型（[0,1]范围），便于后续计算%%小波分解准备dwtmode('zpd');%设置离散小波变换模式为"零填充"%%单层二维离散小波分解[cA,cH,cV,cD]=dwt2(Y,'coif3');%输出4个子带：%cA-低频近似分量（图像轮廓）%cH-水平细节分量（水平方向边缘）%cV-垂直细节分量（垂直方向边缘）%cD-对角细节分量（斜向边缘和纹理）%%显示分解后的子带图像figure('Name','分解后的子带图像');subplot(2,2,1);imshow(mat2gray(cA),[]);title('低频近似(cA)');subplot(2,2,2);imshow(mat2gray(abs(cH)),[]);title('水平细节(cH)');subplot(2,2,3);imshow(mat2gray(abs(cV)),[]);title('