小波变换

我给你大概标注了一下，但是你的程序有问题，

% 小波图像压缩 - RGB 图像

clear all

close all

% 读取图像

im = input('输入图像')%输入图像名称，要加分号

X=imread(im)

% 输入要分解的小波层数和小波

n=input('输入要分解的小波层数')%输入所要分解的层数

wname = input('输入小波名称')%输入小波名称，也要加分号

x = double(X)

NbColors = 255

map = gray(NbColors)

x = uint8(x)

%把RGB图像转换成灰度图

% x = double(X)

% xrgb = 0.2990*x(:,:,1) + 0.5870*x(:,:,2) + 0.1140*x(:,:,3)

% colors = 255

% x = wcodemat(xrgb,colors)

% map = pink(colors)

% x = uint8(x)

% 对图像x进行n维小波分解

[c,s] = wavedec2(x,n,wname)

% 使用默认参数选择各层不同的阈值芦辩

alpha = 1.5m = 2.7*prod(s(1,:))

[thr,nkeep] = wdcbm2(c,s,alpha,m)

% 使用上面的阈值和硬阈值处理进行图像压缩

[xd,cxd,sxd,perf0,perfl2] = wdencmp('lvd',c,s,wname,n,thr,'h')

disp('压缩效率')

disp(perf0)

% 重构（下面这个地方有问题，你这里是原始图像小波变换后进行重构，xd才是小波阀值压缩后重构的图像,cxd，sxd，是c，s经过阀值处理后得到的小波分解结构，也就是说xd=waverec2(cxd,sxd,wname)这个wdencmp函数不需要另外进行重构，你下面那些关于重构的都没用，而下面压缩后的缺清图像才是重构后的图像，）

R = waverec2(c,s,wname)

rc = uint8(R)

% 显示原始图像和压缩图陪扮缺像

subplot(221), image(x)

colormap(map)

title('原始图像')

subplot(222), image(xd)

colormap(map)

title('压缩后的图像')

% 显示结果

xlab1 = ['图像压缩后保留能量百分比',num2str(perfl2)]

xlab2 = ['小波阀值压缩后置零系数百分比 ',num2str(perf0), ' %']

xlabel([xlab1 xlab2])

subplot(223), image(rc)

colormap(map)

title('重构图像')

%计算图像大小

disp('原始图像')

imwrite(x,'original.tif')%将图像x保存为original.tif，下同

imfinfo('original.tif')%显示图片original.tif详细信息，下同

disp('压缩后的图像')

imwrite(xd,'compressed.tif')

imfinfo('compressed.tif')

disp('重构后的图像')

imwrite(rc,'decompressed.tif')

imfinfo('decompressed.tif')

生活中需要对一些图像进行处理，比如压缩，去噪，图像增强，图像锐化与钝化，图像融合，图像的分解等，以便对于图像的成分，边缘等细节信息有更加深刻的认识，小波分析由于其固有的时频特性，既可以对图像进行时域分析，也可以对图像进行频率分析，这使得小波分析在图像处理中得到了广泛的应用，本节对其中一些图像处理功能及函数进行讲解：

wavedec2函数用于对图像进行二维小波分解，其函数调用格式如下：

[c,l]=wavedec2(X,n,’wname’)

其中,X表示原始图像，n表示分解层数，wname表示小波函数，c表示各层系数，l表示各层系数对应的长度

ddencmp用于得到全局阀值，其调用格式如下：

[thr,sorh,keepapp]=ddencmp(‘cmp’,’wp’,X)

[thr,sorh,keepapp]=ddencmp(‘cmp’,’wv’,X)

其中cmp表示压缩，wp表示小波包，wv表示小波，X表示原始信号，thr表示阀值，sorh表示阀值类型，s表示软阀值，h表示硬阀值,keepapp=1表示保持近似系数不变

wdencmp用于对数据或图像进行阀值去噪或压缩，其调用格式如下：

[xcomp,c1,l1,perf0,perfl2]=wdencmp(‘gbl’,c,l,’wname’,n,thr,sorh,keepapp)

glb表示利用全局阀值，perf0表示恢复比，perfl2表示压缩比

示例：利用二维小波对图像进行压缩

编写对应的m文件如下：

clc

load woman

subplot(1,2,1)

imshow(X,map)

title('原始图像')

[c,l]=wavedec2(X,3,'sym4')

%%获取全指模局阀值%%

[thr,sorh,keepapp]=ddencmp('cmp','wp',X)

[xcmp,c1,l1,perf0,perfl2]=wdencmp('gbl',c,l,'sym4',3,thr,sorh,keepapp)

subplot(1,2,2)

imshow(xcmp,map)

title('压缩后图片')

程序运行结果如下图：

小波变换用与图像去噪，噪声会影响图像处理的输入，采集，处理的各个环节及输出结果等全过程，因此对于图像的噪声处理是一个不可忽略的重要的问题，去噪已经成为图像处理中不可或缺的一部分

示例：对图像进行二维小波去噪

编写对应的m文件如下：

load julia

%%产生噪声信号%%

init=3718025452

rand('seed',init)

xnoise=X+8*rand(size(X))

colormap(map)

subplot(1,3,1)

imshow(X,map)

title('原始信号')

subplot(1,3,2)

imshow(xnoise,map)

title('含有噪声的信号')

%%获取全局阀值%%

[thr,sorh,keepapp]=ddencmp('den','wp',xnoise)

[xden,c1,l1]=wdencmp('gbl',xnoise,'sym4'如消,3,thr,sorh,keepapp)

subplot(1,3,3)

imshow(xden,map)

title('去除噪声后信号')

程序运行结果如下图：

小波分析用于图像增强，图像增强是对图像进行一定处理，使图像比原图更加清晰，视觉渣逗知效果更好。

示例：利用小波分析对图像进行增强

编写对应的m文件如下：

clc

load facets

subplot(1,2,1)

imshow(X,map)

title('原始信号')

[c,l]=wavedec2(X,3,'sym4')

sizec=size(c)

fori=1:sizec(2)

if(c(i)>250)

c(i)=2*c(i)

else

c(i)=0.5*c(i)

end

end

y=waverec2(c,l,'sym4')

subplot(1,2,2)

imshow(y,map)

title('增强图像')

程序运行结果如下图：

图像钝化

图像的钝化可以在时域中，也可以在频域中，在时域中处理较为简单，只需要加一个平滑滤波器，使图像中每个点与其邻点做平滑处理即可，在此主要说明图像钝化在频域中的处理。图像钝化是为了突出低频信息，弱化高频信息。

示例：对图像进行频域钝化处理，

编写对应的m文件如下：

load chess

subplot(1,2,1)

imshow(X,map)

title('原始图像')

[c,l]=wavedec2(X,3,'db4')

sizec=size(c)

fori=1:sizec(2)

if(c(i)>280)

c(i)=c(i)*2

else

c(i)=c(i)*0.5

end

end

y=waverec2(c,l,'db4')

subplot(1,2,2)

imshow(y,map)

title('采用小波方法钝化图像')

程序运行结果如下图：

图像锐化，与图像钝化刚好相反，是为了突出高频信息，弱化低频信息，从快速变化的成分中分离出系统边界成分，以便进一步识别或者分割等 *** 作。

示例：对图像进行锐化处理

编写对应的m文件如下：

load chess

subplot(1,2,1)

imshow(X,map)

title('原始图像')

[c,l]=wavedec2(X,3,'db5')

sizec=size(c)

%%突出高频信息，弱化低频信息%%

fori=1:sizec(2)

if(abs(c(i))<280)

c(i)=c(i)*2

else

c(i)=c(i)*0.5

end

end

y=waverec2(c,l,'db5')

subplot(1,2,2)

imshow(y,map)

title('采用小波方法锐化图像')

程序运行结果如下图：

小波分析用于图像融合

图像融合是将同一图像的两个部分或者不同图像合成一张图，以便合成之后的图形比原来更容易理解。

示例：利用二维小波变换将两幅图像融合在一起

编写对应的m文件如下：

clear all

load bust

X1=X

map1=map

load woman

X2=X

map2=map

subplot(1,3,1)

imshow(X1,map1)

title('第一幅图像')

subplot(1,3,2)

imshow(X2,map2)

title('第二幅图像')

%%对第二幅图形低频部分和高频部分进行处理%%

fori=1:256

forj=1:256

if(X2(i,j)>120)

X2(i,j)=X2(i,j)*2

else

X2(i,j)=X2(i,j)*0.5

end

end

end

[c1,l1]=wavedec2(X1,3,'sym4')

[c2,l2]=wavedec2(X2,3,'sym4')

%%对图像进行融合%%

c=c1+c2

%%减少图像的亮度%%

c=c*0.5

y=waverec2(c,l1,'sym4')

subplot(1,3,3)

imshow(y,map2)

title('融合后图像')

程序运行结果如下图：

小波分析用于图像分解

对图像分解的目地在于可以更好的观察图像的细节，对图像做出更好的判断,swt2函数用于对图像进行分解，其调用格式如下：

[sa,sh,sv,sd]=swt2(X,N,’wname’)

其中sa，sh,sv,sd分别表示近似系数，水平系数，竖直系数，对角系数，x分解图像，N分解的层数，wname表示小波基名称

示例：对图像进行分解

编写对应的m文件如下：

clear all

load woman

[sa,sh,sv,sd]=swt2(X,3,'db3')

s=1

fori=1:3

subplot(3,4,s)

image(wcodemat(sa(:,:,i),192))

title(['第',num2str(i),'层近似系数'])

subplot(3,4,s+1)

image(wcodemat(sh(:,:,i),192))

title(['第',num2str(i),'层水平系数'])

subplot(3,4,s+2)

image(wcodemat(sv(:,:,i),192))

title(['第',num2str(i),'层竖直系数'])

subplot(3,4,s+3)

image(wcodemat(sd(:,:,i),192))

title(['第',num2str(i),'层对角系数'])

s=s+4

end

程序运行结果如下图：

设a为尺度，fs为采样频率，Fc为小波中心频率，则a对应的实际频率Fa为

Fa＝Fc×fs/a (1)

显然，为使小波尺度图的频率范围为(0,fs/2)，尺度范围应为(2*Fc,inf),其中inf表示为无穷大。在实际应用中，只需取尺度足够大即可。

没有说明和大家没有讲到是因为没人在小波变换中才去“确定”它们，说白了这两参数不是让你去确定的，是让你去设定的，你对小波的应用方法就没搞清，你要搞清楚的是它们如何影响处理结果的，按照你处理的目的设置不同的值，这牵扯到小波基的某些数学指标是如何影响处理结果的问题，要说的就多了。

    在cmorwavf函数的帮助文档中，列举了cmor1.5-1的函数波形，

中心频率（fc）可以这样看，从横轴0开始的波峰到横轴1的波峰，刚好是正弦波的一个完整周期，其经历的时间就应是频率值的倒数，那么中心频率刚好是1.

下面是cmor1.5-2的函数波形，从横轴0开始的波峰到横轴0.5的波峰也是一个完整周期，经历的时间为0.5，取倒数，中心频率刚好是2.这与你设定的fc一致，也就是说fc就是这么影响Morlet复小波的。当你要消噪或研究高频信息，对于同一个数据信号cmor1.5-2肯定比cmor1.5-1更能消除细小的噪声和得到更高频率的信息。

fb越大时域宽度越长，支撑长度越长，产生高幅值的小波系数也多，在检测信号奇异性的时候往往希望能有一定数量的波峰波谷（在小波中就意味着较长的支撑和较滑姿高的消失矩），当然也不是越多越好，这要看待分析信号的情况，所以这玩意是你先设定，做完CCWT后，看看效果是否满意，再来根据你要研究信号或处理的目的，更改fc和fb的值，不是开始就确定它们（再说在处理之前你如何确定，即使你确定了又有啥意义，在CCWT之前的一切确定是毫无意义的，只有出了结果反复修改你的设定才能最终用“确定”一词）。

总结：fc的大小影响的是小波的频率：fc越大，小波频率越大，因此当你要消噪或研究高频信息，fc增大的话更能消除细小的噪声和得到更高频率的信息。

fb影响的是小波的支撑长度，fb越大，时域宽度越长，支撑长度也越长，产生高幅值的小波系数也就越多，

－、绘制原理

1.需要用到的小波工具箱中的三个函数

COEFS = cwt(S,SCALES,'wname')

说明：该函数能实现 连续小波变换 ，其中S为输入信号，SCALES为尺度，wname为小波名

称信卖绝。

FREQ = centfrq('wname')

说明：该函数能求出以wname命名的母小波的中心频率。

F = scal2frq(A,'wname',DELTA)

说明：该函数能将尺度转换为实际频率，其中A为尺度，wname为小波名称，DELTA为采样

周期。

注：这三个函数还有其它格式，具体可参阅 matlab 的帮助文档。

2.尺度与频率之间的关系

设a为尺度，fs为采样频率，Fc为小波中心频率，则a对应的实际频率Fa为

Fa＝Fc×fs/a

(1)

显然，为使小波尺度图的频率范围为(0,fs/2)，尺度范围应为(2*Fc,inf),其中inf表示

为无穷大。在实际应用中，只需取尺度足够大即可。

3.尺度序列的确定

由式(1)可以看出，为使转换后的频率序列是一等差序列，尺度序列必须取为以下形式：

c/totalscal,...,c/(totalscal-1),c/4,c/2,c

(2)

其中，totalscal是对信号进行 小波变换 时所用尺度序列的长度(通常需要预先设定好)，

c为一常数。

下面讲讲c的求法。

根据式(1)容易看出，尺度c/totalscal所对应的实际频率应为fs/2，于是可得

c=2×Fc/totalscal

(3)

将式(3)代入式(2)便得到了所需的尺度序列。

4.时频图的绘制

确定了小波基和尺度后，就可以用cwt求小波系数coefs(系数是复数时要取模)，然后

用scal2frq将尺度序列转换为实际频率序列f，

最后结合时间序列t，用imagesc(t,f,abs(coefs))便能画出小波时频图。

注意：直接将尺度序列取为等差序列，例如1:1:64，将只能配培得到正确的尺度－时间－小

波系数图，而无法将其转换为频率－时间－小波系数图。这是因为此时的频率间隔不为

常数。

此时，可通过查表的方法将尺度转化为频率或直接修改尺度轴标注。同理，利用本帖所

介绍的方法只能得到频率－时间－小波系数图，不能得到正确的尺度－时间－小波系数

图。

二、应用例子

下面给出一实际例子来说明小波时频图的绘制。所取仿真信号是由频率分别为100Hz和2

00Hz的两个正弦分量所合成的信号。

clear

clc

fs=1024%采样频率

f1=100

f2=200

t=0:1/fs:1

s=sin(2*pi*f1*t)+sin(2*pi*f2*t)

%两个不同频率正弦信号合成的仿真信号

%%%%%%%%%%%%%%%%%小波时频图绘制%%%%%%%%%%%%%%%%%%

wavename='cmor3-3'

totalscal=256

%尺度序列的长度，即scal的长度

wcf=centfrq(wavename)

%小波的中心频率

cparam=2*wcf*totalscal

%为得到合适的尺度所求出的参数

a=totalscal:-1:1

scal=cparam./a

%得到各个尺度，以使转换得到频率序列为等差序列

coefs=cwt(s,scal,wavename)

%得到小波系数

f=scal2frq(scal,wavename,1/fs)

%将尺度转换为频率

imagesc(t,f,abs(coefs))

%绘制色谱图

colorbar

xlabel('时间 t/s')

ylabel('频率 f/Hz')

title('小波时频图')

程序运行结果如下：

说明：(1)应用时只须改变wavename和totalscal两个参数即可。

(2)在这个例子中，最好选用复的morlet小波，其它小波的分析效果不好，而且morlet小

波的带宽参数和中心频率取得越大，时频图上反映的时频聚集性越好。

首先生成一个信号：

设置fc=0.1:0.5:8,fb=10:10:80,观察生成的时频图，查看参数对小波变换的影响。

1. 采样频率和信号点数之间的关系的影响

当采样频率fs与信号点数相同时，即倍数相同时，发现：

时间上为1s,  频率上显示的时真实值，这个时候时频图比较完美。此时参数为：fs=2^16,N=2^16,fc=1.5,fb=3,totalscal=2^7=128

增大fs,当fs为信号点数的两倍时，此时时间上为0.5s,频率坐标范围增大1倍，频率显示的仍是真实值。

减小fs,当fs为信号点数的一半时，此时时间上为2s,频率坐标范围减小1倍，频率显示的仍是真实值。仔细观察，此时的时频图上显示的有杂波。

   如果接着减小，当减小到原信号的1/4时，频率轴坐标为采样频率的一半，因此频率轴范围跟着缩小相同的倍数，这时由于采样频率小于信号的最大频率，此时频率轴显示的并不是频率的真实值。而且时频图出现了错误。所以 采样频率一定要大于信号的最大频率的2倍以上。

下面观察采样频率与信号点数之间的关系，尝试增大采样频率来找到与信号点数之间的限制。

     当信号点数与信号的最大频率近似相等时，设置采样频率为信号点数的2倍时，此时效果并不好，一方面有杂波，另一方面，采样频率与信号最高频率的两倍相离过小。

   当信号点数与信号的最大频率近似相等时，设置采样频率为信号点数的4倍时，此时效果比2倍时要好一点，但存在高频信号分辨率低，有少量杂波存在的问题

当采样频率增加过大时，会出现下列情况，效果反而不好

首先以下列参数设置为基准

fc的值增大从1.5增大到2，发现出现了严重的杂波：

将fc的值增大从1.5增大到3：

fc增大到5

fc增大到7：

当fc从1.5减少到1时，低频处有杂波出现而且频率分辨率明显降低。

当从1.5减少到0.5时：

     下面四张图分别是fb取20，40，80，120的值时的时频图，从中可以看出，fb的取值增大可以增大频率分辨率，但是不像fc的值那样敏感，当增大的范围过大时，也会在不同频率上出现杂波，但是相比fc变动引起的杂波来说很小。

     fb值减小时，频率分辨率会降低，有杂波出现，但是和增大fb时一样，杂波成分分布广但是较小。

     fb值减小时，频率分辨率会降低，有杂波出现，但是和增大fb时一样，杂波成分分布广但是较小。下图fb取值为1.

基准：fs=2^15,N=2^15,fc=1.5,fb=3,totalscal=2^7=128  

增大尺度参数：从128增加到256

从128增加到512

从128增加到1024：

可以看出，当增大尺度值时，低频杂波分量出现，且数值较大，直接盖过信号频率分量

下面是减小尺度值的情况：

从128减小到64：

从128减小到32：

减小到16：

可以减尺度值，分辨率逐渐变差，但是无低频杂波分量出现。

---