小波应用与算法期末试题湖大研究生

更新时间：2023-11-02 02:45:01 阅读量：综合文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

小波理论及应用推荐度：
相关推荐

1. 什么是一维小波变换？相对于传统信号处理方法它有什么特点？为什么要对信号或图

像作多尺度分析？

小波变换是信号处理、图像压缩和模式识别等诸多领域中一个非常有效的数学分析工具，它是一种信号的时间—尺度（时间—频率）分析方法，它具有多分辨率分析的特点，而且与傅里叶变换不同，它具有时频两域都具有表征信号局部特征的能力。在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于探测正常信号中夹带的瞬变反常信号。小波应用于信号去噪已取得了极大的成功。但是在小波阈值去噪中信号边缘处出现的振荡现象一直困扰着小波去噪的应用。小波阈值去噪方法在边缘处产生振荡的原因是由于阈值去噪所采用的小波变换算法没有平移不变性，并且去噪结果依赖于小波基函数、尺度基函数与信号的空间结构的匹配程度，多尺度边缘检测方法能减轻阈值对边缘检测的负面影响，用强度阈值来提取重要边缘并剔出噪声。

2. 编程实现信号的多尺度分解与重构的快速算法（不使用函数dwt和idwt）。

一维小波分解的程序：

function [cA,cD] = mydwt(x,lpd,hpd,dim);

% 函数 [cA,cD]=MYDWT(X,LPD,HPD,DIM) 对输入序列x进行一维离散小波分解，输出分解序列[cA,cD]

% 输入参数：x——输入序列； % lpd——低通滤波器； % hpd——高通滤波器； % dim——小波分解级数。

% 输出参数：cA——平均部分的小波分解系数； % cD——细节部分的小波分解系数。 cA=x; % 初始化cA，cD cD=[];

for i=1:dim

cvl=conv(cA,lpd); % 低通滤波，为了提高运行速度，调用MATLAB提供的卷积函数conv()

dnl=downspl(cvl); % 通过下抽样求出平均部分的分解系数 cvh=conv(cA,hpd); % 高通滤波

dnh=downspl(cvh); % 通过下抽样求出本层分解后的细节部分系数 cA=dnl; % 下抽样后的平均部分系数进入下一层分解

cD=[cD,dnh]; % 将本层分解所得的细节部分系数存入序列cD end

function y=downspl(x);

% 函数 Y=DOWMSPL(X) 对输入序列进行下抽样，输出序列 Y。

% 下抽样是对输入序列取其偶数位，舍弃奇数位。例如 x=[x1,x2,x3,x4,x5]，则 y=[x2,x4]. N=length(x); % 读取输入序列长度

M=floor(N/2); % 输出序列的长度是输入序列长度的一半（带小数时取整数部分） i=1:M; y(i)=x(2*i);

而重构则是分解的逆过程，对低频系数、高频系数分别进行上抽样和低通、高通滤波处理。要注意重构时同一级的低频、高频系数的个数必须相等。

function y = myidwt(cA,cD,lpr,hpr);

% 函数 MYIDWT() 对输入的小波分解系数进行逆离散小波变换，重构出信号序列 y % 输入参数：cA —— 平均部分的小波分解系数； % cD —— 细节部分的小波分解系数；

% lpr、hpr —— 重构所用的低通、高通滤波器。

lca=length(cA); % 求出平均、细节部分分解系数的长度 lcd=length(cD);

while (lcd)>=(lca) % 每一层重构中，cA 和 cD 的长度要相等，故每层重构后， % 若lcd小于lca，则重构停止，这时的 cA 即为重构信号序列 y 。 upl=upspl(cA); % 对平均部分系数进行上抽样 cvl=conv(upl,lpr); % 低通卷积

cD_up=cD(lcd-lca+1:lcd); % 取出本层重构所需的细节部分系数，长度与本层平均部分系数的长度相等

uph=upspl(cD_up); % 对细节部分系数进行上抽样 cvh=conv(uph,hpr); % 高通卷积

cA=cvl+cvh; % 用本层重构的序列更新cA，以进行下一层重构 cD=cD(1:lcd-lca); % 舍弃本层重构用到的细节部分系数，更新cD

lca=length(cA); % 求出下一层重构所用的平均、细节部分系数的长度 lcd=length(cD);

end % lcd < lca，重构完成，结束循环

y=cA; % 输出的重构序列 y 等于重构完成后的平均部分系数序列 cA function y=upspl(x);

% 函数 Y = UPSPL(X) 对输入的一维序列x进行上抽样，即对序列x每个元素之间 % 插零，例如 x=[x1,x2,x3,x4],上抽样后为 y=[x1,0,x2,0,x3,0,x4]; N=length(x); % 读取输入序列长度

M=2*N-1; % 输出序列的长度是输入序列长度的2倍再减一

for i=1:M % 输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列元素 if mod(i,2)

y(i)=x((i+1)/2); else y(i)=0; end end

3. 编程实现一维信号的多尺度突变点检测。

% Mann-Kendall突变检测 % 数据序列y

% 结果序列UFk，UBk2

%-------------------------------------------- %读取excel中的数据，赋给矩阵y %获取y的样本数

%A为时间和径流数据列

A=xlswrite('数据.xls') x=A(:,1);%时间序列 y=A(:,2);%径流数据列 N=length(y); n=length(y);

% 正序列计算---------------------------------

% 定义累计量序列Sk，长度=y，初始值=0 Sk=zeros(size(y));

% 定义统计量UFk，长度=y，初始值=0 UFk=zeros(size(y)); % 定义Sk序列元素s s = 0;

% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0 % 此时UFk无意义，因此公式中，令UFk(1)=0 for i=2:n for j=1:i

if y(i)>y(j) s=s+1; else s=s+0; end; end; Sk(i)=s;

E=i*(i-1)/4; % Sk(i)的均值

Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差 UFk(i)=(Sk(i)-E)/sqrt(Var); end;

% ------------------------------正序列计算end % 逆序列计算--------------------------------- % 构造逆序列y2，长度=y，初始值=0 y2=zeros(size(y));

% 定义逆序累计量序列Sk2，长度=y，初始值=0 Sk2=zeros(size(y));

% 定义逆序统计量UBk，长度=y，初始值=0 UBk=zeros(size(y));

% s归0 s=0;

% 按时间序列逆转样本y

% 也可以使用y2=flipud(y);或者y2=flipdim(y,1); for i=1:n

y2(i)=y(n-i+1); end;

% i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var(1)均为0 % 此时UBk无意义，因此公式中，令UBk(1)=0

for i=2:n for j=1:i

if y2(i)>y2(j) s=s+1; else s=s+0; end; end; Sk2(i)=s;

E=i*(i-1)/4; % Sk2(i)的均值

Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差

%由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1， %则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反

% 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i)) % 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势 UBk(i)=0-(Sk2(i)-E)/sqrt(Var); end;

% ------------------------------逆序列计算end

% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量 % 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此 % 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用 UBk2=zeros(size(y));

% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1); for i=1:n

UBk2(i)=UBk(n-i+1); end;

% 做突变检测图时，使用UFk和UBk2 % 写入目标xls文件：f:\\test2.xls % 目标表单：Sheet1

% 目标区域：UFk从A1开始，UBk2从B1开始 xlswrite('f:\\test2.xls',UFk,'Sheet1','A1'); xlswrite('f:\\test2.xls',UBk2,'Sheet1','B1'); figure(3)%画图

plot(x,UFk,'r-','linewidth',1.5); hold on

plot(x,UBk2,'b-.','linewidth',1.5);

plot(x,1.96*ones(N,1),':','linewidth',1); axis([min(x),max(x),-5,5]);

legend('UF统计量','UB统计量','0.05显著水平');

xlabel('t (year)','FontName','TimesNewRoman','FontSize',12); ylabel('统计量','FontName','TimesNewRoman','Fontsize',12); %grid on hold on

plot(x,0*ones(N,1),'-.','linewidth',1); plot(x,1.96*ones(N,1),':','linewidth',1); plot(x,-1.96*ones(N,1),':','linewidth',1);

4. 编程实现图像的张量积分解与重构的快速算法（不使用函数dwt2和idwt2）。

首先，为了实现图像的N层分解，对一幅m行n列的黑白图像，我们要对其进行规范化处理，使其能被2的N次方整除。以下的modmat() 函数实现此功能： function y=modmat(x,dim)

% 函数 MODMAT() 对输入矩阵x进行规范化，使其行列数均能被 2^dim 整除 % 输入参数：x —— r*c 维矩阵； % dim —— 矩阵重构的维数

% 输出参数：y —— rt*ct 维矩阵，mod(rt,2^dim)=0，mod(ct,2^dim)=0

[row,col]=size(x); % 求出输入矩阵的行列数row,col

rt=row - mod(row,2^dim); % 将row,col分别减去本身模 2^dim 得到的数 ct=col - mod(col,2^dim); % 所得的差为rt、ct，均能被 2^dim 整除 y=x(1:rt,1:ct); % 输出矩阵 y 为输入矩阵 x 的 rt*ct 维子矩阵

然后，将规范化后的图像的数据格式由适合显示图像的uint8格式转换为适合数值处理的double格式，再调用二维小波分解函数进行图像分解，最后为了清晰地显示分解图像的塔式结构，在图像的相应区域绘制若干分界线。具体程序如下： function y=mywavedec2(x,dim)

% 函数 MYWAVEDEC2() 对输入矩阵 x 进行 dim 层分解，得到相应的分解系数矩阵 y % 输入参数：x —— 输入矩阵； % dim —— 分解层数。

% 输出参数：y —— 分解系数矩阵

x=modmat(x,dim); % 首先规范化输入矩阵，使其行列数均能被 2^dim 整除 subplot(121);imshow(x);title('原始图像'); % 画出规范化后的源图像 [m,n]=size(x); % 求出规范化矩阵x的行列数

xd=double(x); % 将矩阵x的数据格式转换为适合数值处理的double格式 for i=1:dim

xd=modmat(xd,1);

[dLL,dHL,dLH,dHH]=mydwt2(xd); % 矩阵小波分解

tmp=[dLL,dHL;dLH,dHH]; % 将分解系数存入缓存矩阵

xd=dLL; % 将缓存矩阵左上角部分的子矩阵作为下一层分解的源矩阵 [row,col]=size(tmp); % 求出缓存矩阵的行列数

y(1:row,1:col)=tmp; % 将缓存矩阵存入输出矩阵的相应行列 end

yd=uint8(y); % 将输出矩阵的数据格式转换为适合显示图像的uint8格式 for i=1:dim % 对矩阵 yd 进行分界线处理，画出分解图像的分界线 m=m-mod(m,2); n=n-mod(n,2); yd(m/2,1:n)=255; yd(1:m,n/2)=255; m=m/2;n=n/2;

end

subplot(122);imshow(yd);title([ num2str(dim) ' 维小波分解图像']); 以下是程序的运行结果。

图2

上述的图像分解程序，其输出数据是double格式的，以便作为重构程序的输入。接下来我们讨论图像的重构。我编写的重构程序中，为了比较分解图像和重构图像，首先绘出经过小波分解的图像，然后再进行重构。在绘制分解图像和重构图像的过程中，要注意数据格式的转换。 function y=mywaverec2(x,dim)

% 函数 MYWAVEREC2() 对输入的分解系数矩阵x进行 dim 层重构，得到重构矩阵 y % 输入参数：x ——分解系数矩阵； % dim ——重构层数； % 输出参数：y ——重构矩阵。 % 绘制分解图像

xd=uint8(x); % 将输入矩阵的数据格式转换为适合显示图像的uint8格式 [m,n]=size(x); % 求出输入矩阵的行列数

for i=1:dim % 对转换矩阵xd进行分界线处理 m=m-mod(m,2); n=n-mod(n,2); xd(m/2,1:n)=255; xd(1:m,n/2)=255; m=m/2;n=n/2; end figure;

subplot(121);imshow(xd);title([ num2str(dim) ' 层小波分解图像']); % 画出带有分界线的分解图像 % 重构图像

xr=double(x); % 将输入矩阵的数据格式转换回适合数值处理的double格式 [row,col]=size(xr); % 求出转换矩阵xr的行列数

for i=dim:-1:1 % 重构次序是从内层往外层进行，所以先抽取矩阵 xr 的最内层分解矩阵进行重构

tmp=xr(1:floor(row/2^(i-1)),1:floor(col/2^(i-1))); % 重构的内层矩阵的行列数均为矩阵xr的2^(i-1)

[rt1,ct1]=size(tmp); % 读取待重构矩阵 tmp 的行列数 rt=rt1-mod(rt1,2);ct=ct1-mod(ct1,2);

rLL=tmp(1:rt/2,1:ct/2); % 将待重构矩阵 tmp 分解为四个部分 rHL=tmp(1:rt/2,ct/2+1:ct); rLH=tmp(rt/2+1:rt,1:ct/2); rHH=tmp(rt/2+1:rt,ct/2+1:ct);

tmp(1:rt,1:ct)=myidwt2(rLL,rHL,rLH,rHH); % 将重构结果返回到矩阵 tmp

xr(1:rt1,1:ct1)=tmp; % 把矩阵 tmp 的数据返回到矩阵 xr 的相应区域，准备下一个外层的重构 end

y=xr; % 重构结束后得到的矩阵xr即为输出矩阵 y