第5章图像正交变换

信号不仅可以在空间域表示，也可以在频率域表示，后者将有利于许多问题的分析讨论，因此，常需要对信号进行正交变换，变换到频率域进行处理。对图像进行正交变换在图像增强、图像复原、图像特征提取、图像编码等处理中都经常采用。常用的
图像正交变换有多种，本章主要介绍图像的离散傅里叶变换、离散余弦变换、KL变换、Radon变换和小波变换。
5.1离散傅里叶变换
离散傅里叶变换(Discrete Fourier Transform，DFT)是直接处理离散时间信号的傅里叶变换，在数字信号处理中应用广泛。本节学习离散傅里叶变换的定义、性质、实现及其在图像处理中的应用。
5.1.1离散傅里叶变换的定义

对于有限长数字序列f(x)，x=0，1，…，N-1，一维离散傅里叶变换及离散傅里叶反变换
（IDFT）定义为： 



F(u)=∑N-1x=0f(x)e-j2πuxN
u=0，1，2，…，N-1

f(x)=1N∑N-1u=0F(u)ej2πuxN
x=0，1，2，…，N-1(51)



f(x)和F(u)为离散傅里叶变换对，表示为
F［f(x)］=F(u)或f(x)F(u)。
设W=e-j2πN，则一维的离散傅里叶变换和
离散傅里叶反变换可以表示为：



F(u)=∑N-1x=0f(x)Wuxu=0，1，2，…，N-1

f(x)=1N∑N-1u=0F(u)W-uxx=0，1，2，…，N-1(52)


数字图像为二维数据，二维离散傅里叶变换是由一维离散傅里叶变换推广而来的。二维离散傅里叶变换和离散傅里叶反变换定义为： 




F(u，v)=∑M-1x=0∑N-1y=0f(x，y)e-j2π
xuM+yvN
x，u=0，1，2，…，M-1

f(x，y)=1MN∑M-1u=0∑N-1v=0F(u，v)ej2πxuM+yvN

y，v=0，1，2，…，N-1(53)



其中，f(x，y)是二维离散信号，F(u，v)为f(x，y)的频谱，u、v为频域采样值； f(x，y)和F(u，v)为二维离散傅里叶变换对，
记为F［f(x，y)］=F(u，v)或f(x，y)F(u，v)。
F(u，v)一般为复数，表示为： 



F(u，v)=R(u，v)+jI(u，v)=|F(u，v)|ejφ(u，v)(54)


其中，|F(u，v)|称为f(x，y)的傅里叶谱，如式（55）所示； (u，v)称为f(x，y)的相位谱
，如式（56）所示； f(x，y)的功率谱定义为傅里叶谱的平方，如式(57)所示。



|F(u，v)|=R2(u，v)+I2(u，v)(55)

(u，v)=arctanI(u，v)R(u，v)(56)

E(u，v)=|F(u，v)|2=R2(u，v)+I2(u，v)(57)






5.1.2离散傅里叶变换的实现
将二维离散傅里叶变换变换式做如下变换： 


F(u，v)=∑M-1x=0∑N-1y=0f(x，y)e-j2πxuMe-j2πyvN

=∑M-1x=0∑N-1y=0f(x，y)e-j2πyvNe-j2πxuM

=∑M-1x=0Fy［f(x，y)］e-j2πxuM

=FxFy［f(x，y)］
(58)

式(58)称为二维离散傅里叶变换的可分性，表明二维离散傅里叶变换可用一维离散傅里叶变换来实现，即先对f(x，y)的每一列进行一维离散傅里叶变换，得到Fy［f(x，y)］，再对该中间结果的每一行进行一维离散傅里叶变换，得到F(u，v)，运算过程中每个一维离散傅里叶变换可以采用快速傅里叶变换(Fast Fourier Transformation，FFT)实现快速运算。相反的顺序(先行后列)也可以。
MATLAB中提供了实现离散傅里叶变换的相关函数，列举如下。

(1) fft函数： 实现一维离散傅里叶变换。
① fft(X)： 对向量X进行离散傅里叶变换。若X为二维矩阵，则对它的每一列进行变换； 若为N维矩阵，则对其第一个长度非1的维进行变换。
② fft(X，N)： N点离散傅里叶变换，若X长度小于N，则补0； 若X长度大于N，则截断。
③ fft(X，［］，DIM)或者fft(X，N，DIM)： 在第DIM维进行离散傅里叶变换。
(2) ifft函数： 实现一维离散傅里叶反变换。
① ifft(X)： X的离散傅里叶反变换。
② ifft(X，N)： 通过补0或截断实现X的N点离散傅里叶反变换。
③ ifft(X，［］，DIM)或者ifft(X，N，DIM)： X在DIM维的离散傅里叶反变换。
(3) fft2函数： 实现二维离散傅里叶变换。
① fft2(X)： 对矩阵X实现二维离散傅里叶变换。
② fft2(X，MROWS，NCOLS)： 通过补0或截断对X进行MROWS×NCOLS的二维离散傅里叶变换。
(4) ifft2函数： 实现二维离散傅里叶反变换。
① ifft2(F)： 对矩阵F实现二维离散傅里叶反变换。
② ifft2(F，MROWS，NCOLS)： 通过补0或截断对F进行MROWS×NCOLS的二维离散傅里叶反变换。
（5） fftshift(X)： 将零频率分量移动到频谱中心。
（6） ifftshift(X)： 将零频率分量移回原位，取消fftshift的效果。
【例51】对灰度图像进行离散傅里叶变换并显示频谱图。
程序如下： 

clear,clc,close all;

grayI=imread('cameraman.tif');%打开灰度图像

DFT=fft2(grayI);%计算离散傅里叶变换

ADFT=abs(DFT);%计算傅里叶谱

top=max(ADFT(:));

bottom=min(ADFT(:));

ADFT1=(ADFT-bottom)/(top-bottom)*100;%把傅里叶谱系数规格化到［0 100］，便于观察

ADFT2=fftshift(ADFT1);%将规格化频谱图移位，低频移至频谱图中心

subplot(131),imshow(Image),title('原图');

subplot(132),imshow(ADFT1),title('原频谱图');

subplot(133),imshow(ADFT2),title('移位频谱图');


程序运行结果如图51所示。图51(a)是原图； 将傅里叶谱规格化到［0 100］的频谱图如图51(b)所示，四角部分对应低频成分，中央部分对应高频成分； 采用fftshift函数将频谱图进行移位，如图51(c)所示，频谱图中间为低频部分，越靠外边频率越高。图像中的能量主要集中在低频区，高频能量很少或为0。


图51灰度图像傅里叶频谱图


【例52】对彩色图像进行离散傅里叶变换并显示频谱图。

彩色图像有3个色彩通道，其数据为3个二维矩阵，因此，需要进行3个二维离散傅里叶变换，将3个频谱图合成彩色频谱图。fft2函数对于M×N×3的矩阵，即按照上述过程进行离散傅里叶变换，变换后的频谱矩阵同样为M×N×3的矩阵，为彩色频谱图。
程序如下： 

clear,clc,close all;

Image=imread('desert.jpg');%彩色图像

DFT=fftshift(fft2(Image)); %离散傅里叶变换并搬移频谱

DFT=abs(DFT); %取傅里叶谱

ADFT=(DFT-min(DFT(:)))/(max(DFT(:))-min(DFT(:)))*100;

subplot(121),imshow(Image),title('原图');

subplot(122),imshow(ADFT),title('彩色频谱图');

程序运行结果如图52所示。从MATLAB工作区可以看到Image为231×352×3的uint8型矩阵，离散傅里叶变换同样是231×352×3的矩阵。


图52彩色图像傅里叶频谱图


【例53】对图像进行离散傅里叶变换、
显示频谱图和重建。
程序如下： 

clear,clc,close all;

Image=imread('peppers.jpg');

DFT=fftshift(fft2(Image));

ADFT=abs(DFT);

top=max(ADFT(:));

bottom=min(ADFT(:));

ADFT=(ADFT-bottom)/(top-bottom)*100;%计算傅里叶谱

recI=ifft2(ifftshift(DFT));%离散傅里叶反变换

recI=abs(recI);%取复数的模

recI=uint8(recI);%变double型数据为整型数据，以便显示图像

subplot(131),imshow(Image),title('原图');

subplot(132),imshow(ADFT),title('彩色频谱图');

subplot(133),imshow(recI),title('DFT重建图');

程序运行结果如图53所示。


图53离散傅里叶变换变换并重建


5.1.3离散傅里叶变换的性质
离散傅里叶变换有许多重要性质，这些性质给离散傅里叶变换的运算和实际应用提供了极大的便利。这里主要介绍几个和二维离散傅里叶变换在图像处理中的应用密切相关的性质。
1. 线性和周期性

若F［f(x，y)］=F(u，v)，0≤x，u≤M，0≤y，v≤N，则



F［a1f1(x，y)+a2f2(x，y)］=a1F［f1(x，y)］+a2F［f2(x，y)］(59)


F(u，v)=F(u+M，v)=F(u，v+N)=F(u+M，v+N)


f(x，y)=f(x+M，y)=f(x，y+N)=f(x+M，y+N)
(510)


式(510)表明，尽管F(u，v)对无穷多个u和v的值重复出现，但只需根据在任一个周期里的值就可从F(u，v)得到f(x，y)； 同样只需一个周期里的变换就可将F(u，v)在频域里完全确定。
2. 几何变换性
1) 共轭对称性
若F［f(x，y)］=F(u，v)，F*(-u，-v)是f(-x，-y)的离散傅里叶变换的共轭函数，则


F(u，v)=F*(-u，-v)(511)


2) 平移性
若F［f(x，y)］=F(u，v)，则


f(x-x0，y-y0)F(u，v)e-j2πx0uM+y0vN

f(x，y)ej2πxu0M+yv0NF(u-u0，v-v0)
(512)

式(512)表示对图像平移不影响其傅里叶变换的幅值，只改变相位谱。当u0=M/2，v0=N/2时，ej2π(u0x/M+v0y/N)=ejπ(x+y)=(-1)x+y，则f(x，y)(-1)x+yF(u-M/2，v-N/2)，频域的坐标原点从起始点(0，0)移至中心点，只要将f(x，y)乘以(-1)x+y因子，再进行傅里叶变换即可实现，即
例51中的频谱搬移。
3) 旋转性
把f(x，y)和F(u，v)表示为极坐标形式，若f(γ，θ)F(k，φ)，则


f(γ，θ+θ0)F(k，φ+θ0)(513)


若将空间域函数旋转角度θ0，那么在变换域此函数的离散傅里叶变换也旋转同样的角度； 反之，若将变换域函数旋转某一角度，则空间域函数也旋转同样的角度。
4) 比例变换特性
若F［f(x，y)］=F(u，v)，则


f(ax，by)1|ab|Fua，vb(514)


对图像f(x，y)在空间尺度的缩放导致其傅里叶变换F(u，v)在频域尺度的相反缩放。

【例54】对一幅图像进行几何变换，再进行离散傅里叶变换运算，验证以上性质。
程序如下： 

clear,clc,close all;

Image=rgb2gray(imread('block.bmp'));

［h,w,c］=size(Image);

scale=imresize(Image,0.5,'bilinear');%缩小变换

rotate=imrotate(Image,30,'bilinear','crop');%旋转变换

tform=affine2d(［1 0 0;0 1 0;20 20 1］);

R=imref2d(［h,w］,［1 w］,［1 h］);

trans=imwarp(Image,tform,'FillValue',0,'OutputView',R);%平移变换

Originaldft=abs(fftshift(fft2(Image)));%原图离散傅里叶变换

Scaledft=abs(fftshift(fft2(scale)));%缩小图离散傅里叶变换

Rotatedft=abs(fftshift(fft2(rotate)));%旋转图离散傅里叶变换

Transdft=abs(fftshift(fft2(trans)));%平移图离散傅里叶变换

figure,imshow(Image),title('原图');

figure,imshow(scale),title('缩小变换');

figure,imshow(rotate),title('旋转变换');

figure,imshow(trans),title('平移变换'); 

figure,imshow(Originaldft,［］),title('原图DFT');

figure,imshow(Scaledft,［］),title('比例变换DFT');

figure,imshow(Rotatedft,［］),title('旋转变换DFT');

figure,imshow(Transdft,［］),title('平移变换DFT');


程序运行结果如图54所示。可以看出，缩小变换后图像的频谱图尺度展宽，旋转后图像的频谱图随着旋转，平移后图像的频谱图没有变化。


图54图像几何变换及其傅里叶频谱图


3. Parseval定理

若F［f(x，y)］=F(u，v)，则


∑M-1x=0∑N-1y=0|f(x，y)|2=∑M-1u=0∑N-1v=0|F(u，v)|2(515)



Parseval定理也称为能量守恒定理，这个性质说明信号经傅里叶变换后不损失能量，只是改变了信号的表现形式，是变换编码的基本条件。
【例55】对图像进行离散傅里叶变换并计算时域和频域能量。
程序如下： 

clear,clc,close all;

grayI=double(imread('cameraman.tif'));

［N,M,C］=size(grayI);

SI=grayI.^2;%计算|fx,y|2

energyT=sum(SI(:));%时域能量

DFT=abs(fft2(grayI));%离散傅里叶变换

SDFT=DFT.^2;%计算|Fu,v|2

energyF=sum(SDFT(:));%频域能量

energyF=energyF/(N*M);%频域能量除以采样点数,消除采样增益

diff=abs(energyF-energySI);%时域、频域能量差

运行程序，在MATLAB工作区可以看到diff值约为0，即时域和频域能量一致。
4. 卷积定理

若F［f(x，y)］=F(u，v)，F［g(x，y)］=G(u，v)，则


f(x，y)*g(x，y)F(u，v)·G(u，v)

f(x，y)·g(x，y)F(u，v)*G(u，v)
(516)


在以上几个性质中，可分性使得二维离散傅里叶变换可以通过一维快速傅里叶变换快速实现； 共轭对称性、平移性、旋转性、比例变换特性使得二维离散傅里叶变换具有一定的几何变换不变性，可以作为一种图像特征； Parseval定理是变换编码的基本条件； 卷积定理可以降低某些复杂图像处理算法的计算量，这几个性质在图像处理中应用较多。
5.1.4离散傅里叶变换在图像处理中的应用
离散傅里叶变换在图像处理中的应用主要包括用于描述图像信息、滤波、压缩及便于运算几个方面。

1. 傅里叶描绘子
从原始图像中产生的数值、符号或者图形称为图像特征，反映了原图像的重要信息和主要特性，以便于让计算机有效地识别目标。这些表征图像特征的一系列符号称为描绘子。

描绘子应具有几何变换不变性，即在图像内容不变，仅产生几何变换(平移、旋转、缩放等)的情况下，描绘子不变，以保证识别结果的稳定性。离散傅里叶变换在图像特征提取方面应用较多，本小节主要介绍离散傅里叶变换直接作为特征的应用
——傅里叶描绘子。

一个闭合区域，区域边界上的点(x，y)用复数表示为
x+jy。沿边界跟踪一周，得到一个复数序列z(n)=x(n)+jy(n)，n=0，1，…，N-1，z(n)为周期信号，其离散傅里叶变换系数用Z(k)表示，Z(k)称为傅里叶描绘子。

根据离散傅里叶变换特性，Z(k)幅值具有旋转和平移不变性，相位信息具有缩放不变性，在一定程度上满足了描绘子的几何变换不变性，可以作为一种图像特征。

2. 离散傅里叶变换在图像滤波中的应用

经过离散傅里叶变换后，在傅里叶频谱中，中间部分为低频部分，越靠外边频率越高。因此，可以在离散傅里叶变换后，设计相应的滤波器，实现高通滤波、低通滤波等的处理。
【例56】对图像进行离散傅里叶变换及频域滤波。
程序如下： 

clear,clc,close all;

Image=imread('desert.jpg');

grayIn=rgb2gray(Image);

［h,w］=size(grayIn);

DFTI=fftshift(fft2(grayIn));%离散傅里叶变换及频谱搬移

cf=30;%截止频率

HDFTI=DFTI;

HDFTI(h/2-cf:h/2+cf,w/2-cf:w/2+cf)=0;%低频置为0

grayOut1=uint8(abs(ifft2(ifftshift(HDFTI))));%离散傅里叶反变换

LDFTI=zeros(h,w);

LDFTI(h/2-cf:h/2+cf,w/2-cf:w/2+cf)=DFTI(h/2-cf:h/2+cf,w/2-cf:w/2+cf);  %高频置为0

grayOut2=uint8(abs(ifft2(ifftshift(LDFTI))));%离散傅里叶反变换

subplot(131),imshow(Image),title('原图');

subplot(132),imshow(uint8(grayOut1)),title('高通滤波');

subplot(133),imshow(uint8(grayOut2)),title('低通滤波');

程序运行结果如图55所示。程序中进行的实际是理想高通和低通滤波器，为提高滤波性能，也有其他类型的滤波器，详见第6章相关内容。


图55图像频域滤波示例


3. 离散傅里叶变换在图像压缩中的应用

由Parseval定理知，信号经傅里叶变换前后能量不发生损失，只是改变了信号的表现形式，离散傅里叶变换系数表现的是各个频率点上的幅值。高频反映细节、低频反映景物概貌，往往认为可将高频系数置为0，降低数据量； 同时由于人眼的惰性，合理地设置高频系数为0，可使图像质量一定范围内的降低不会被人眼察觉到。因此，离散傅里叶变换可以方便地进行压缩编码。
4. 离散傅里叶变换卷积性质的应用
抽象来看，图像处理算法可以认为是图像信息经过了滤波器的滤波(如平滑滤波、锐化滤波等)，空间域滤波通常需要进行卷积运算。如果滤波器的结构比较复杂，可以利用离散傅里叶变换的卷积性质，把空间域卷积变为变换域的相乘，以简化运算，如式(517)所示。



fg=g*f

Fg(u，v)=G(u，v)·F(u，v)

fg=IDFT(Fg)

(517)



式中，f为原图像，g为滤波器，利用g对f滤波，是用g和f卷积得到fg。这个过程可以改变为先对f、g进行离散傅里叶变换，把G和F相乘得Fg，再进行傅里叶反变换，以降低卷积计算量。

注意： 由于离散傅里叶变换和离散傅里叶反变换都是周期函数，因此在计算卷积时，需要让这两个离散函数具有同样的周期，否则将产生错误。利用快速傅里叶变换计算卷积，为防止频谱混叠误差，需对离散的二维函数补零，即周期延拓，两个函数同时周期延拓，使其具有相同的周期。

5.2离散余弦变换

在傅里叶级数展开式中，如果被展开的函数是实偶函数，那么其傅里叶级数中只包含余弦项，再将其离散化可导出余弦变换，因此称为DCT(Discrete Cosine Transform，离散余弦变换)。
5.2.1离散余弦变换的定义
1. 一维DCT


对于有限长数字序列f(x)，x=0，1，…，N-1，其一维DCT
和
IDCT（离散余弦反变换）定义为：



F(u)=C(u)2N∑N-1x=0f(x)cos(2x+1)uπ2N

f(x)=2N∑N-1u=0C(u)F(u)cos(2x+1)uπ2N
(518)



其中，x，u=0，1，…，N-1，C(u)=
1/2u=0


1u=1，2，…，N-1



将一维DCT变换表示成矩阵运算形式，即


F=Af(519)

A=2N1212…12

cos12Nπcos32Nπ…cos(2N-1)2Nπ

︙︙︙

cosN-12Nπcos3(N-1)2Nπ…cos(2N-1)(N-1)2Nπ(520)


其中，F为变换系数矩阵，A为正交变换矩阵，f为时域数据矩阵。

一维DCT逆变换的矩阵形式表示为： 


f=ATF(521)


2. 二维DCT

数字图像为二维数据，把一维DCT推广到二维，二维DCT变换和反变换定义为： 



F(u，v)=2MNC(u)C(v)∑M-1x=0∑N-1y=0f(x，y)cosπ(2x+1)u2Mcos
π(2y+1)v2N

f(x，y)=2MN∑M-1u=0∑N-1v=0C(u)C(v)F(u，v)cosπ(2x+1)u2Mcos
π(2y+1)v2N

(522)



式中，x，u=0，1，2，…，M-1


y，v=0，1，2，…，N-1，C(u)，C(v)=1/2u，v=0


1，其他
二维DCT的矩阵形式表示为： 


F=AfAT(523)


二维IDCT的矩阵形式表示为： 



f=ATFA(524)


其中，F为DCT变换系数矩阵，f为空域数据矩阵，A为正交变换矩阵。
5.2.2离散余弦变换的实现
【例57】用矩阵运算对图像进行DCT变换。
程序如下： 

clear,clc,close all;

Image=double(imread('cameraman.tif'));

［N,M,C］=size(Image);%图像尺寸，本例中图像为256×256，N和M相等

A=zeros(N,N);

A(1,1:N)=1/sqrt(2);%系数矩阵A的第1行

for j=2:N

for i=1:N

A(j,i)=cos((2*i-1)*pi*(j-1)/(2*N));

end

end

A=sqrt(2/N)*A;%生成系数矩阵A

DCT=A*Image*A';%二维DCT矩阵运算

DCT=abs(DCT);

top=max(DCT(:));

bottom=min(DCT(:));

DCT=(DCT-bottom)/(top-bottom)*100;%DCT系数规格化

imshow(DCT),title('DCT频谱图');

程序运行结果如图56所示。从所得的结果可以看出，离散余弦变换具有使信息集中的特点，
对图像进行DCT变换后，在变换域中，矩阵左上角低频的幅值大，而右下角高频的幅值小。


图56图像DCT变换



MATLAB中提供了实现DCT变换的相关函数，列举如下。
(1) dct函数： 实现DCT变换。
① Y=dct(X)： 对向量X进行DCT变换，返回Y； 若X是二维矩阵，则对其每一列进行DCT变换； 若X是N维矩阵，则对其第1个长度非1的维进行变换。

② Y=dct(X，N)： 通过补0或截断对向量X进行N点DCT变换。
③ Y=dct(X，［］，DIM)或Y=dct(X，N，DIM)： 对X的第DIM维进行DCT变换。

④ Y=dct(…，'Type'，K)： 指定计算方式执行DCT变换，K可取1、2、3、4，分别代表dctI、dctII、dctIII、dctIV变换，默认值为2。
(2) idct函数： 实现IDCT变换。
① X=idct(Y)： IDCT变换。
② X=idct(Y，N)： 通过补0或截断对向量Y进行N点IDCT变换。
③ X=idct(Y，［］，DIM)或X=idct(Y，N，DIM)： 对Y的第DIM维进行IDCT变换。
④  X=idct(…，'Type'，K)： 指定计算方式执行IDCT变换，K可取1、2、3、4，分别代表idctI、idctII、idctIII、idctIV变换，默认值为2。
(3) dct2函数： 实现二维DCT变换。
① B=dct2(A)： 对矩阵A进行二维DCT变换，返回矩阵B。
② B=dct2(A，［M N］)或B=dct2(A，M，N)： 通过补0或截断对矩阵A进行M×N的DCT变换。
（4） idct2函数： 实现二维IDCT变换。
① B=idct2(A)： 对矩阵A进行二维IDCT变换。
② B=idct2(A，［M N］)或 B=idct2(A，M，N)： 通过补0或截断对矩阵A进行M×N的IDCT变换。
（5） D=dctmtx(N)： 返回N×N的DCT变换矩阵D(即式(520)中的变换矩阵A)。

（6） B=blockproc(A，［M N］，FUN)： 对图像A中的每一M×N块执行FUN定义的操作。函数将图像A中的每一数据块打包为"block struct"传递给用户定义的函数FUN，返回矩阵(向量、标量)Y，Y=FUN(BLOCK_STRUCT)。"block struct"是MATLAB中定义的结构体，包括该块的信息，相关字段如表51所示。


表51block struct结构体字段



参数含义

border二维向量［V H］，指定数据块的水平、垂直边界
blockSize二维向量［rows cols］，指定块大小，若'border'被指定，块大小不包括边界像素
dataM×N或者M×N×P的块数据矩阵
imageSize二维向量［rows cols］，指定输入图像大小
location二维向量［row col］，指定输入图像中数据块的第1个像素(第1行第1列)，不包括边界像素

【例58】对灰度图像进行DCT变换并显示频谱图。
程序如下： 

clear,clc,close all;

grayI=imread('cameraman.tif');

DCT=dct2(grayI); %进行离散余弦变换

ADCT=abs(DCT);%求模

top=max(ADCT(:));

bottom=min(ADCT(:));

result=(ADCT-bottom)/(top-bottom)*100; %把模规格化到［0 100］

subplot(121),imshow(grayI),title('原灰度图');

subplot(122),imshow(result),title('灰度图DCT频谱图');

程序运行结果如图57所示。从3幅灰度图的频谱可以看出，能量主要集中在左上角低频分量处。


图57灰度图像的DCT频谱图



【例59】对彩色图像进行DCT变换、显示频谱图并重建原图。
程序如下： 

clear,clc,close all;

RGBI=imread('peppers.jpg');
%打开彩色图像

［N,M,C］=size(RGBI);

DCTS=zeros(N,M,C);%DCT频谱图初始化

RGBOut=zeros(N,M,C);%重建图初始化

for i=1:C%分色彩通道处理

channel=RGBI(:,:,i);

DCT=dct2(channel);%DCT变换

ADCT=abs(DCT);

top=max(ADCT(:));bottom=min(ADCT(:));

DCTS(:,:,i)=(ADCT-bottom)/(top-bottom)*100;

RGBOut(:,:,i)=abs(idct2(DCT));%IDCT变换

end

RGBOut=uint8(RGBOut);

subplot(131),imshow(RGBI),title('原彩色图');

subplot(132),imshow(DCTS),title('彩色图DCT频谱图');

subplot(133),imshow(RGBOut),title('DCT重建图');

程序运行结果如图58所示。


图58彩色图像DCT变换并重建


5.2.3离散余弦变换在图像处理中的应用

离散余弦变换在图像处理中主要用于对图像(包括静止图像和运动图像)进行有损数据压缩。例如，静止图像编码标准JPEG、运动图像编码标准MPEG中都使用了离散余弦变换，这是由于离散余弦变换具有很强的“能量集中”特性，大多数的能量都集中在离散余弦变换后的低频部分，压缩编码效果较好。
具体的做法一般是先把图像分成8×8的块，对每一个方块进行二维DCT变换，变换后的能量主要集中在低频区。对DCT系数进行量化，给高频系数大间隔量化，低频部分小间隔量化，舍弃绝大部分取值很小或为0的高频数据，降低数据量，同时保证重构图像不会发生显著失真。
【例510】对图像进行DCT变换，将高频置为0，并进行反变换重建。
程序如下： 

clear,clc,close all;

grayIn=imread('cameraman.tif');

［h,w］=size(grayIn);

DCTI=dct2(grayIn);

cf=90;%截止频率

FDCTI=zeros(h,w);

FDCTI(1:cf,1:cf)=DCTI(1:cf,1:cf);%将高频系数置0

grayOut=uint8(abs(idct2(FDCTI)));

subplot(121),imshow(grayIn),title('原图');

subplot(122),imshow(grayOut),title('压缩重建');

程序运行结果如图59所示。因为丢失了部分高频系数，所以重建图比原图模糊。程序中将大于截止频率的高频系数置0，并没有考虑系数的大小，截止频率的指定也缺乏依据。


图59DCT压缩重建示例一



【例511】将彩色图像转换到YCbCr空间，分别对亮度和色度数据进行DCT变换，采用合理的方式将高频系数变为0，并进行反变换重建。
本例中通过将DCT系数除以合适的数据以便将高频系数变为0，称为量化； 各DCT系数除以的数据称为量化步长，将在第11章具体学习。
程序如下： 

clear,clc,close all;

ImageIn=imread('peppers.jpg');
%读取彩色图像

YCbCrIn=double(rgb2ycbcr(ImageIn));%转换到YCbCr空间

YQT=［1611101624405161;

1212141926586055;

1413162440576956;

1417222951878062;

182237566810910377;

243555648110411392;

49647887103121120101;

7292959811210010399］;%亮度量化表，即块内各数据对应的量化步长

CQT=［1718244799999999;

1821266699999999;

2426569999999999;

4766999999999999;

9999999999999999;

9999999999999999;

9999999999999999;

9999999999999999］;%色度量化表

blocksize=8;%定义块大小

A=dctmtx(blocksize);%计算变换矩阵

FUN1=@(block_struct) A*block_struct.data*A';%定义块DCT变换函数

FUN2=@(block_struct) A'*block_struct.data*A;%定义块IDCT变换函数

YCbCrOut=zeros(size(YCbCrIn));%对输出的YCbCr图像初始化

for i=1:3

channel=YCbCrIn(:,:,i);%获取YCbCr各通道

DCT=blockproc(channel,［blocksize,blocksize］,FUN1);%块DCT变换

if i==1

QT=YQT;

else

QT=CQT;

End%选择亮度或色度量化表

FUN3=@(block_struct) round(block_struct.data./QT);%定义量化函数，点除运算

QDCT=blockproc(DCT,［blocksize,blocksize］,FUN3);%块DCT系数量化

FUN4=@(block_struct) block_struct.data.*QT;%定义反量化函数，点乘运算

IQDCT=blockproc(QDCT,［blocksize,blocksize］,FUN4); %块量化后的数据反量化

IDCT=blockproc(IQDCT,［blocksize,blocksize］,FUN2);%块系数IDCT变换

YCbCrOut(:,:,i)=IDCT;%存储YCbCr图像数据

end

YCbCrOut=uint8(YCbCrOut);

ImageOut=ycbcr2rgb(YCbCrOut);%变换回RGB空间

figure,imshow(ImageIn),title('原图');

figure,imshow(ImageOut),title('重建图像');

程序运行结果如图510所示。程序中实现了色彩空间变换、DCT变换、量化、反量化、IDCT变换、色彩空间逆变换等操作，数据损失发生在量化一步。从图510可以看出，重建后的图像质量有所下降，但不明显，比直接指定截止频率更合理。


图510DCT压缩重建示例二

5.3KL变换
KL变换是建立在统计特性基础上的一种变换，又称为霍特林(Hotelling)变换或主成分分析(Principal Component Analysis，PCA)，其突出优点是相关性好，是均方误差(Mean Square Error，MSE)意义下的最佳变换，在数据压缩技术中占有重要地位。
5.3.1KL变换原理
1. KL变换的定义


用确定的正交归一向量系uj，j=1，2，…，∞展开向量X： 


X=∑∝j=1ajuj(525)


用有限的m项来估计向量X，即


X^=∑mj=1ajuj(526)


表示成矩阵形式为： 


X^=UA(527)


为了找到合适的变换矩阵U，计算用有限项展开代替无限项展开引起的均方误差，使均方误差最小的U最优。均方误差如下所示： 


ε2=E［(X-X^)T(X-X^)］=E∑∞j=m+1ajuj·∑∞j=m+1ajuj


利用已知条件求解均方误差。u为正交归一向量系，
uTiuj=1i=j


0i≠j； 且aj=uTjX，所以


ε2=E∑∞j=m+1a2j
=E∑∞j=m+1uTjXXTuj
=∑∞j=m+1uTjE［XXT］uj


令ψ=E［XXT］，则


ε2=∑∞j=m+1uTjψuj


利用拉格朗日乘数法求均方误差取极值时的u，拉格朗日函数为： 


h(uj)=∑∞j=m+1
uTjψuj-∑∞
j=m+1λ［uTjuj-1］


对uj求导数，得


(ψ-λjI)uj=0，j=m+1，…，∞


其中，λj是X的自相关矩阵ψ的特征值，uj是对应的特征向量。

则有： 


ε2=∑∞j=m+1uTjψuj=∑∞j=m+1uTjλj
uj=∑∞j=m+1λj(528)


得到以下结论： 
以X的自相关矩阵ψ的m个最大特征值对应的特征向量来逼近X时，其截断均方误差具有极小性质。这m个特征向量所组成的正交坐标系U称作X所在的n维空间的m维KL变换坐标系。X在KL坐标系上的展开系数向量
A称作X的KL变换，满足： 



A=UTX


X=UA(529)


其中，U=(u1u2…um)。

2. KL变换的性质

因ψuj=λjuj，则ψU=UDλ，Dλ为对角矩阵，其互相关成分都应为0，即



Dλ=λ10…0

0λ2…0

︙︙︙

00…λn(530)


因U为正交矩阵，所以ψ=UDλUT。
因X=UA，则ψ=E［XXT］=
E［UAATUT］=UE［AAT］UT，所以，



E［AAT］=Dλ(531)


由式(531)可知，变换后的向量A的自相关矩阵ψA是对角矩阵，且对角元素是X的自相关矩阵ψ的特征值。显然，通过KL变换，消除了原有向量X的各分量之间的相关性，即变换后的数据A的各分量之间的信息是相互独立的。
3. 信息量分析
通过前文的分析可知，数据X的KL坐标系的产生矩阵采用的是自相关矩阵ψ=E［XXT］，由于总体均值向量μ常常没有什么意义，
因此也常常把数据的协方差矩阵作为KL坐标系的产生矩阵，如式（532）所示。



Σ=E［(X-μ)(X-μ)T］(532)


已知a1=uT1X，
计算a1的方差如下： 


var(a1)=E［a21］-E［a1］2=E［uT1XXTu1］-E［uT1X］E［XTu1］=uT1Σu1


u1为Σ的特征向量，λ1为对应的特征值，则var(a1)=uT1Σu1=λ1uT1u1=λ1，即a1的方差为Σ最大的特征值，a1称作第一主成分。

采用大特征值对应的特征向量组成变换矩阵，能对应地保留原向量中方差最大的成分，KL变换起到了减小相关性、突出差异性的效果，称之为主成分分析。
计算主成分的方差之和为： 


∑nj=1var(aj)=∑nj=1λj=tr(Σ)=∑nj=1var(Xj)


上式说明，n个互不相关的主成分的方差之和等于原数据的总方差，即n个互不相关的主成分包含了原数据中的全部信息，各主成分的贡献率依次递减，第一主成分贡献率最大，数据的大部分信息集中在较少的几个主成分上。
主成分ai的贡献率为： 


λi∑nj=1λj，i=1，2，…，n(533)


前m个主成分的累积贡献率反映前m个主成分综合原始变量信息的能力，定义如下： 


∑mi=1λi∑nj=1λj(534)


【例512】
设向量集为
ω1： (000)T，(101)T，(100)T，(110)T

ω2： (001)T，(011)T，(010)T，(111)T，采用其自相关矩阵作为产生矩阵对其进行KL变换，实现向二维降维，即求其前两个主成分。
程序如下： 

clear,clc,close all;

X=［0 0 0;1 0 1;1 0 0;1 1 0;0 0 1;0 1 1;0 1 0;1 1 1］';

［n,N］=size(X);%n为维数，N为样本数

V=X*X'/N;%求自相关矩阵

［coeff,D］=eigs(V);%求特征值和特征向量

［D_sort,index］=sort(diag(D),'descend');%特征值降序排列

D=D(index,index);%按序调整特征值对角矩阵

coeff=coeff(:,index);%按序调整特征向量矩阵

score=coeff'*X;%K-L变换

figure; plot(score(1,:),score(2,:),'ko'),title('K-L变换');

xlabel('第一主成分得分');ylabel('第二主成分得分');

reconstructed=score'*coeff';%K-L逆变换

程序运行结果如图511所示。


图511KL变换第一、二主成分


5.3.2图像的KL变换及其实现

图像的KL变换是指将二维的图像转化为一维的向量，通常采用奇异值分解进行KL变换。
1. 原理
图像KL变换的原理是
将二维图像采用行堆叠或列堆叠转换为一维处理。设一幅大小为M×N的图像f(x，y)在某个传输通道上传输了L次，由于受到各种因素的随机干扰，接收的图像是一个图像集合，即


{f1(x，y)，f2(x，y)，…，fL(x，y)}


采用列堆叠将每一个M×N的图像表示为MN维的向量，即


fi=(fi(0，0)fi(0，1)…fi(M-1，N-1))T


图像向量f=(f1f2…fL)的协方差矩阵和相应变换核矩阵为


Σf=E［(f-μf)(f-μf)T］≈1L∑Li=1fi
fTi-μfμTf(535)


其中，f-μf为原始图像f减去平均值向量μf，称为中心化图像向量； Σf是MN×MN维的矩阵。
设λi和ui为Σf的特征值和特征向量，且降序排列，即


λ1>λ2>λ3>λ4>…>λM×N


KL变换矩阵U为


U=(u1，u2，…，uM×N)=u11u21…uMN1

u12u22…uMN2

︙︙︙

u1MNu2MN…uMNMN


则二维KL变换表示为：


F=UT(f-μf)(536)


离散KL变换向量
F是中心化向量f-μf与变换核矩阵U相乘所得的结果。
2. 奇异值分解
如前文所述，f向量的协方差矩阵Σf是MN×MN维的矩阵，由于图像的维数MN一般很高，
因此直接求解Σf的特征值和特征向量不现实。本小节简单介绍奇异值分解(Singular Value Decomposition，SVD)的方法，其详细的数学理论可以参看矩阵论的相关资料。
1) 原理
奇异值分解将一个大矩阵分解为几个小矩阵乘积，如式（537）所示。


B=PDQT(537)


其中，B为m×n的矩阵； 
P为m×m的方阵，其列向量正交，称为左奇异向量； 
D为m×n的矩阵，仅对角线元素不为0，对角线上的元素称为奇异值； Q
T为n×n的方阵，其列向量正交，称为右奇异向量。
式(537)中小矩阵的求解可以采用下列方法。
设R = BTB，得到一个方阵，且RT = (BTB)T = BTB=R，即R为n阶厄米特矩阵，可以证明R的特征值均为非负值。对矩阵R求特征值，如式（538）所示。


(BTB)qi=λiqi(538)


右奇异矩阵Q由qi组成。
由式(539)可得左奇异矩阵P。



σi=λi


pi=Bqi/σi(539)


其中，左奇异矩阵P由pi组成； σ即矩阵B的奇异值，在矩阵D中从大到小排列，且减小很快，可以用前r个大的奇异值来近似描述矩阵，
如式（540）所示。


Bm×n≈Pm×rDr×rQTn×r(540)


需注意，Q为n×r的矩阵，QT为r×n的矩阵。
2) 图像KL变换实现
将中心化图像向量f-μf进行奇异值分解，即B=f-μf，用前r个大的奇异值来近似描述，如式（541）所示。


BMN×L≈PMN×rDr×rQTL×r(541)


将式(541)两边同时右乘QL×r，即BMN×LQL×r≈PMN×rDr×rQTL×rQL×r。

由于Q为正交矩阵，所以QTQ为单位阵，所以： 



BMN×LQL×r≈PMN×rDr×r=
B～MN×r(542)


由式(538)求出矩阵Q，进而求出B～MN×r，实现列压缩。
将式(541)两边同时左乘PTMN×r，即PTMN×rBMN×L≈PTMN×rPMN×rDr×rQTL×r。
由于P为正交矩阵，所以PTP为单位阵，所以： 



PTMN×rBMN×L≈Dr×rQTL×r=B～r×L(543)


由式(538)和式(539)求出矩阵Q、D、P，进而求出B～r×L，实现行压缩。
【例513】打开人脸图像，采用SVD方法对其进行KL变换，并显示变换结果。
程序如下： 

clear,clc,close all;

fmt={'*.jpg','JPEG image(*.jpg)';'*.*','All Files(*.*)'};

［FileName,FilePath］=uigetfile(fmt,'导入数据','face*.jpg','MultiSelect','on');

%利用uigetfile函数交互式选取训练样本图像

if ~isequal(［FileName,FilePath］,［0,0］)

FileFullName=strcat(FilePath,FileName);

else 

return; %若选择了图像，则获取图像文件的完整路径，否则退出程序，不再运行

end 

N=length(FileFullName); 

for k=1:N

Image=im2double(rgb2gray(imread(FileFullName{k})));

X(:,k) = Image(:);%把图像放在矩阵X的第k列

InImage(:,:,:,k)=Image;

end

figure,montage(InImage,'size',［1 NaN］),title('原图');

［h,w,c］=size(Image);

averagex = mean(X')';%计算图像的平均向量μ

X=X-averagex;%求中心化图像向量

R = X'*X;%奇异值分解中的矩阵
R = BTB

［Q,D］ = eig(R);%求矩阵R的特征值和特征向量

［D_sort,index］ = sort(diag(D),'descend');%特征值从大到小排序

D=D(index,index);

Q = Q(:,index);%按从大到小顺序重排特征值矩阵D和特征向量矩阵Q

P = X*Q*(abs(D))^-0.5;%求左奇异矩阵P

total = 0.0;

count = sum(D_sort);

for r =1:N

total = total + D_sort(r);

if total/count > 0.95%取占全部奇异值之和95%的前r个奇异值

break;

end

end

KLCoefR = P'*X;

figure; plot(KLCoefR(1,:),KLCoefR(2,:),'ko'),title('K-L变换行压缩');

xlabel('第一主成分得分');ylabel('第二主成分得分');

Y= P(:,1:2)*KLCoefR(1:2,:)+averagex;%基于前2个奇异值重建人脸图像

outImage=zeros(h,w,1,N);

for j=1:N

outImage(:,:,1,j)=reshape(Y(:,j),h,w);

end

figure,montage(outImage,'size',［2 NaN］),title('基于左奇异矩阵前两个奇异值重建图像');

Z= P(:,1:r)*KLCoefR(1:r,:)+averagex;%基于前r个奇异值重建人脸图像

for j=1:N

outImage(:,:,1,j)=reshape(Z(:,j),h,w);

end

figure,montage(outImage,'size',［1 NaN］,'DisplayRange',［］);

title('基于左奇异矩阵前r个奇异值重建的人脸图像');

KLCoefC = X*Q; %使用右奇异矩阵进行K-L变换

for j =1:N

outImage(:,:,1,j)=reshape(KLCoefC(:,j),h,w);

end

figure,montage(outImage,'size',［1 NaN］,'DisplayRange',［］);

title('基于右奇异矩阵的K-L变换');

程序运行结果如图512所示。


图512人脸图像KL变换


3. KL变换函数

MATLAB中对于KL变换提供了相应的函数，列举如下。
(1) pcacov函数： 根据相关系数矩阵或协方差矩阵进行主成分分析。
① COEFF=pcacov(V)： 根据n×n的矩阵V进行主成分分析，返回主成分系数COEFF。COEFF为n×n的矩阵，每一列为一个主成分的系数向量，按主成分方差递减顺序排列，即原理中的变换矩阵U。
② ［COEFF，LATENT］=pcacov(V)： LATENT是由主成分的方差构成的向量，即由V的n个特征值(从大到小)构成的向量。
③ ［COEFF，LATENT，EXPLAINED］=pcacov(V)： EXPLAINED是由n个主成分的百分比贡献率构成的向量。
（2） pca函数： 对原数据进行主成分分析。
① COEFF=pca(X)： 返回N×n的数据矩阵X的主成分系数COEFF，N为样本数目，n为样本维数，X的每一行对应一个样本； COEFF含义同pcacov函数中的参数。默认情况下，pca函数将数据变为中心化数据，使用SVD算法。

② ［COEFF，SCORE］=pca(X)： 增加返回主成分得分SCORE，是原理中矩阵A的转置。需要注意原理讲解中一般都是用列向量表示一个样本，N个n维样本则是一个n×N的矩阵； 但pca函数中输入的X为N×n的数据矩阵，是原理中数据矩阵的转置，SCORE的计算采用X×COEFF，实际是原理中矩阵A的转置； 默认情况下，是先将X变为中心化数据再乘以COEFF计算出SCORE，重建中心化数据则采用SCORE×COEFF'。

③ ［COEFF，SCORE，LATENT］=pca(X)： 增加返回主成分方差LATENT，即X的协方差矩阵的n个特征值构成的向量。

④ ［COEFF，SCORE，LATENT，TSQUARED］=pca(X)： 增加返回X中每个样本的Hotelling T平方统计量。TSQUARED
为有N个元素的向量，第i个元素是第i个样本与数据集的中心之间的距离。



T2i=∑nj=1SCORE2ijλji=1，2，…，N


其中，SCOREij(i=1，2，…，N;j=1，2，…，n)为第i个样本的第j个主成分得分，λj为样本协方差矩阵的n个降序排列的特征值。

⑤ ［COEFF，SCORE，LATENT，TSQUARED，EXPLAINED］=pca(X)： 增加返回各主成分的百分比贡献率组成的向量EXPLAINED。
⑥ ［COEFF，SCORE，LATENT，TSQUARED，EXPLAINED，MU］=pca(X)： 当输入参数'Centered'设为true时，返回X中样本均值MU； 否则MU为全0。
⑦ ［…］=pca(…，'PARAM1'，val1，'PARAM2'，val2，…)： 指定各参数实现主成分分析，用以控制计算、处理特殊数据类型。各参数含义如表52所示。


表52pca函数参数



参数含义

Algorithmpca函数实现主成分分析的算法，可选'svd'(奇异值分解)、'eig'(协方差矩阵特征值分解)、'als'(交替最小二乘法)，默认值为'svd'
Centered逻辑值，默认为true，pca函数将X变为中心化数据； 为false，不对数据进行中心化处理
Economy逻辑值，默认为true，函数仅返回结果的前一部分； 为false，则返回全部结果
NumComponents正整数，指明需要的主成分数目K，返回COEFF和SCORE的前K列
续表



参数含义

RowsX中含有NaN情况时的处理方法，可取'complete'(默认值，计算前去除值为NaN的样本，运算后在SCORE中的对应位置插入)； 'pairwise'(切换Algorithm为'eig'，计算协方差矩阵时用非NaN值取代NaN)； 'all'(确认X中没有数据缺失，采用所有数据进行运算，遇NaN终止运算)。采用'als'算法时本参数无效
Weights样本权向量，元素均为正值
VariableWeights分量权向量，元素均为正值
配合'als'算法的参数略

(3) pcares函数： 重建数据，并求样本观测值矩阵中每个样本的每一个分量所对应的残差。
① RESIDUALS=pcares(X，NDIM)： 返回利用前NDIM个主成分重建数据时的残差。NDIM≤n； RESIDUALS为N×n的矩阵，其元素为X中相应元素所对应的残差。
② ［RESIDUALS，RECONSTRUCTED］=pcares(X，NDIM)： 增加返回用前NDIM个主成分的得分重建的数据RECONSTRUCTED，是X的一个近似。pcares只接受原始样本观测数据作为它的输入，不会自动对数据做标准化变换。

【例514】利用函数对不同字体数字的图像进行主成分分析并重建。
本例中对白色背景下不同字体的黑色数字的图像进行处理，设计分为以下3部分。
(1) 生成样本。读取图像文件，反色，将数字目标变为白色； 通过确定数字的外接矩形截取数字所在区域； 将数字区域归一化为16×16的子图像； 将16×16的数据变换为1×256的一维数据，生成样本。

clc;clear;close all;

fmt={'*.jpg','JPEG (*.jpg)';'*.bmp','BMP (*.bmp)';'*.*','All Files(*.*)'};

［FileName,FilePath］=uigetfile(fmt,'导入外部数据','*.jpg','MultiSelect','on');

if ~isequal(［FileName,FilePath］,［0,0］);

FileFullName=strcat(FilePath,FileName);

else

return;

end%选择数据并判断，同例5-13

N=length(FileFullName);%选择的图像数目

Image=zeros(50);%存储原始图像数据，大小为50×50

n=16;n1=20;%规格化图像为n×n，显示用图像为n1，外围增加了2行2列

BWI=zeros(n);%存储二值化图像数据，大小为16×16

training=zeros(1,n*n);%存储训练样本1×256

showImage=zeros(n1,n1*N);%显示用图像，将各个样本小图像拼到一幅大图显示

for j=1:N%对每一幅图像进行预处理

Image=imread(FileFullName{j});%读取每一幅图像数据

Image=255-Image;%反色，将数字目标变为白色

Image=im2bw(Image,0.4);%二值化

［y,x］=find(Image==1);

BWI=Image(min(y):max(y),min(x):max(x));%截取数字所在区域

BWI=imresize(BWI,［n,n］);%子区域归一化为16×16

showImage((n1-n)/2+1:n+(n1-n)/2,(j-1)*n1+1:(j-1)*n1+n)=BWI;
%当前图像拼到大图中

training(j,:)=double(BWI(:)');%转化为一维的样本

end

(2) 进行主成分分析。采用pca函数进行主成分分析，并通过绘图观察样本的第一、二主成分。

［coeff,score,latent,tsquare］=pca(training); %主成分分析

figure,plot(score(:,1),score(:,2),'ko');%绘制第一、二主成分得分图

title('主成分分析'),xlabel('第一主成分得分'),ylabel('第二主成分得分'); 

(3) 利用主成分重建图像。采用pcares函数分别基于第一、第一、二、第一、二、三主成分重建图像，并通过绘图观察重建图像与原始图像的区别。

［residuals1,reconstructed1］=pcares(training,1);%利用第一主成分重建数据

［residuals2,reconstructed2］=pcares(training,2);  %利用第一、二主成分重建数据

［residuals3,reconstructed3］=pcares(training,3);  %利用第一、二、三主成分重建数据

showResult1=zeros(n1,n1*N);

showResult2=zeros(n1,n1*N);

showResult3=zeros(n1,n1*N);

for j=1:N

tempI=reconstructed1(j,:);%重建数据一

tempI=reshape(tempI,n,n);%转换为二维图像

showResult1((n1-n)/2+1:n+(n1-n)/2,(j-1)*n1+1:(j-1)*n1+n)=tempI;

tempI=reconstructed2(j,:);

tempI=reshape(tempI,n,n);

showResult2((n1-n)/2+1:n+(n1-n)/2,(j-1)*n1+1:(j-1)*n1+n)=tempI;

tempI=reconstructed3(j,:);

tempI=reshape(tempI,n,n);

showResult3((n1-n)/2+1:n+(n1-n)/2,(j-1)*n1+1:(j-1)*n1+n)=tempI;

end

figure,

subplot(221),imshow(showImage),title('原始二值数字图像');

subplot(222),imshow(showResult1),title('第一主成分重建的数字图像');

subplot(223),imshow(showResult2),title('第一、二主成分重建的数字图像');

subplot(224),imshow(showResult3),title('第一、二、三主成分重建的数字图像');

程序运行结果如图513所示，随着重建采用的主成分数目增多，重建图像清晰度逐渐提高。


图513图像主成分分析及重建


5.3.3KL变换在图像处理中的应用

KL变换在图像处理中主要用于图像数据压缩和特征提取。

如前文所述，KL变换矩阵由特征向量组成，特征向量按特征值递减顺序排列。由于能量集中在特征值较大的系数中，因此丢掉特征值小的特征向量构成变换矩阵。KL变换的结果F是原图像f的低维投影，减少了数据量，在保留的主成分的贡献率不低于一定程度的情况下，不影响重建图像的质量。

KL变换常作为一种特征提取方法，从一组特征中计算出一组按重要性从大到小排列的新特征，是原有特征的线性组合，并且相互之间是不相关的，实现了数据的降维。例如，在人脸识别中，可以用KL变换对人脸图像的原始空间进行转换，即构造人脸图像数据集的协方差矩阵，求出协方差矩阵的特征向量，再依据特征值的大小对这些特征向量进行排序，每一个特征向量的维数与原始图像一致，可以看作是一个图像，被称作特征脸。每一个人脸图像都可以确切地表示为一组特征脸的线性组合。
5.4Radon变换
图像的Radon变换也是一种重要的图像处理研究方法，将图像函数f(x，y)沿其所在平面内的不同直线做线积分，即进行投影变换，可以获取图像在该方向上的突出特性，在去噪、重建、检测、复原中多有应用。本节介绍Radon变换的原理、实现及其应用。
5.4.1Radon变换的原理

如图514(a)所示，直线L的方程可以表示为
ρ=xcosθ+ysinθ，其中，ρ代表坐标原点到直线L的距离，θ∈［0，π］是直线法线与x轴的夹角
。要将函数f(x，y)沿直线L做线积分，即进行Radon变换，变换式表示为： 


R(ρ，θ)=∫Lf(x，y)ds(544)



采用Delta函数求解该线积分。Delta函数是一个广义函数，在非零点取值为0，而在整个定义域的积分为1，
用最简单的Delta函数表示，如式(545)所示。


δ(t)=
0t≠0


1t=0(545)



对于直线L，直线上的点(x，y)满足δ(t)=1，非直线上的点满足δ(t)=0，即



δ(xcosθ+ysinθ-ρ)=
0xcosθ+ysinθ-ρ≠0


1xcosθ+ysinθ-ρ=0(546)


Radon变换表达式可更改为： 


R(ρ，θ)=∫∞-∞∫∞-∞f(x，y)δ(xcosθ+ysinθ-ρ)dxdy(547)


其中，R(ρ，θ)是f(x，y)的Radon变换，表示为R［f(x，y)］=R(ρ，θ)。

给定一组(ρ，θ)，即可得出一个沿Lρ，θ的积分值。如图514(b)所示，n条与直线L平行的线
具有相同的θ角，但ρ不同，对每一条线都做f(x，y)的线积分，有n条投影线，即对一幅图像在某一特定角度下的Radon变换会产生n个线积分值，构成一个n维的向量，称为f(x，y)在角度θ下的投影。


图514Radon变换坐标系图


Radon变换可以看成是xy空间向ρθ空间的投影，ρθ空间上的每一点对应xy空间的一条直线。图像中高灰度值的线段会在ρθ空间形成亮点，而低灰度值的线段在ρθ空间形成暗点。因而，对图像中线段的检测可转化为在变换空间对亮点、暗点的检测。

二维Radon变换的反变换如式(548)所示。


f(x，y)=12π2∫π0dθ∫∞-∞R/ρxcosθ+ysinθ-ρdρ(548)


5.4.2Radon变换的实现
在给定θ方向的情况下，数字图像f(x，y)沿直线L的线积分，可以通过坐标系旋转后按列累加来实现。如图515所示，θ是
直线L法线与x轴的夹角，坐标系xoy顺时针旋转θ角变为x′oy′，x′轴与L垂直，y′轴与L平行，将f(x′，y′)沿y′方向求和，实现图像在θ方向上的Radon变换。


图515坐标系旋转示意图


由图515可知，坐标系旋转前后点的对应关系如式（549）所示。



x′=xcosθ+ysinθ


y′=ycosθ-xsinθ(549)


实际是图像的逆时针旋转。
因此，图像的Radon变换可以按下列步骤实现。
(1) 计算图像对角线的长度，即是ρ的取值范围； 
(2) 设定旋转方向θ∈［0，π］； 
(3) 将图像中的点按式(549)变为新坐标系中的点(旋转变换)； 
(4) 将f(x′，y′)沿y′方向求和，即新图像的垂直投影。
【例515】对一幅图像进行指定方向上的Radon变换，并显示变换结果。
程序如下： 

clc;clear;close all;

Image=rgb2gray(imread('block.bmp'));

［N,M］=size(Image);

len=floor(sqrt(N^2+M^2)+1);%对角线长度，ρ的取值范围

x=-len/2:len/2-1;%投影后向量的横坐标范围，用于绘图

R=zeros(1,len);%初始化投影向量

free=floor((len-M+1)/2);%将投影值置于向量中间，求出前后空置的元素数目

R(free+1:free+M)=sum(Image,1);%原图垂直投影，即0°方向上的Radon变换 

%以下为45°方向上的Radon变换

theta1=45;

Image1=imrotate(Image,theta1);

［N1,M1］=size(Image1);

free1=floor((len-M1+1)/2); 

R1=zeros(1,len);

R1(free1+1:free1+M1)=sum(Image1,1);

%以下为70°方向上的Radon变换

theta2=70;

Image2=imrotate(Image,theta2);

［N2,M2］=size(Image2);

free2=floor((len-M2+1)/2);

R2=zeros(1,len);

R2(free2+1:free2+M2)=sum(Image2,1);

subplot(231),imshow(Image),title('原图');

subplot(232),imshow(Image1),title('旋转45°');

subplot(233),imshow(Image2),title('旋转70°');

subplot(234),plot(x,R),title('0°方向上的Radon变换');

subplot(235),plot(x,R1),title('45°方向上的Radon变换');

subplot(236),plot(x,R2),title('70°方向上的Radon变换');

程序运行结果如图516所示。

MATLAB提供了如下实现Radon变换和反变换的函数。



图516Radon变换的实现


(1) R=radon(I，THETA)： 对灰度图像I进行Radon变换。THETA为投影夹角，若为一标量，则R为一列向量； 若THETA为一向量，则R为一矩阵，其每一列对应THETA中一个角度的Radon变换值，默认情况下THETA取值为0:179。
(2) ［R，Xp］=radon(…)： 向量Xp是R每一行对应的径向坐标，以图像的中心像素为原点。
(3) I=iradon(R，THETA)： 利用二维矩阵R中的投影数据重建图像I。THETA指定投影的角度，若为向量，其中的元素为等间距角度； 若为一个标量D_theta，投影进行的角度为m×D_theta(m=0，1，2，…，size(R，2)-1)； 若输入为空矩阵，D_theta默认为180/size(R，2)。
(4) I=iradon(R，THETA，INTERPOLATION，FILTER，FREQUENCY_SCALING，OUTPUT_SIZE)： 指定参数实现Radon反变换，参数如表53所示。


表53iradon函数参数



参数含义

INTERPOLATION插值运算方法，可取'nearest'、'linear'、'spline'、'pchip'，默认为'linear'
FILTER指定在频率域滤波的滤波器，可取'RamLak'、'SheppLogan'、'Cosine'、'Hamming'、'Hann'、'none'，默认为'RamLak'
FREQUENCY_SCALING(0，1］之间的数，用于修正滤波器，默认为1，若小于1，滤波器被修正为适于频域范围［0，FREQUENCY_SCALING］，大于FREQUENCY_SCALING的频率被置为0
OUTPUT_SIZE指定重建图像的尺寸，若未指定，重建图像尺寸为2×floor(size(R，1)/(2×sqrt(2)))

【例516】采用radon函数对图像进行Radon变换，并显示变换结果。
程序如下： 

clear,clc,close all;

Image=rgb2gray(imread('block.bmp'));

［R1,X1］=radon(Image,0);

［R2,X2］=radon(Image,45);

［R3,X3］=radon(Image,70);

subplot(221),imshow(Image),title('原图');

subplot(222),plot(X1,R1),title('0°方向上的Radon变换曲线');

subplot(223),plot(X2,R2),title('45°方向上的Radon变换曲线');

subplot(224),plot(X3,R3),title('70°方向上的Radon变换曲线');

程序运行结果如图517所示。



图517指定方向上的Radon变换


【例517】对一幅图像进行Radon正变换和逆变换，并显示变换结果。
程序如下： 

clear,clc,close all;

Image=rgb2gray(imread('block.bmp'));

theta1=0:0.1:180;%间隔0.1°进行Radon变换

［R1,X1］=radon(Image,theta1);

result1=iradon(R1,theta1);%采用默认设置进行Radon逆变换

theta2=0:1:180; %间隔1°进行Radon变换

［R2,X2］=radon(Image,theta2);

result2=iradon(R2,theta2);

theta3=0:10:180;%间隔10°进行Radon变换

［R3,X3］=radon(Image,theta3);

result3=iradon(R3,theta3);

figure,imshow(Image),title('原图');

figure,colormap(gray);

subplot(231),imagesc(theta1,X1,R1),title('间隔0.1°投影的Radon变换曲线集合');

subplot(232),imagesc(theta2,X2,R2),title('间隔1°投影的Radon变换曲线集合');

subplot(233),imagesc(theta3,X3,R3),title('间隔10°投影的Radon变换曲线集合');

subplot(234),image(result1),title('间隔0.1°投影的重建图像');

subplot(235),image(result2),title('间隔1°投影的重建图像');

subplot(236),image(result3),title('间隔10°投影的重建图像');

程序运行结果如图518所示。R1、R2、R3为二维矩阵，表示多条变换后的曲线，如图518(a)~(c)所示，横轴表示180°，纵轴表示每条曲线的高度，将R中的元素数值按大小转化为不同颜色，在坐标轴对应位置处以该颜色染色； 投影间隔角度越小，投影后的变换曲线越多，染色越细腻。图518(d)~(f)为重建图，与原图有差别，是由于Radon变换过程中数据损失造成的，投影间隔越小，损失越小，重建图像和原图越接近。


图518Radon变换和逆变换


5.4.3Radon变换的应用
Radon变换可用来检测图像中的线段。将原来的
xy平面内的点映射到ρθ平面上，原xy平面一条线段上所有的点都将投影到ρθ平面上的同一点。记录ρθ平面上
点的累积程度，累积程度足够的点所对应的ρθ值即是xy平面上线段的参数。
Radon变换与第9章要讲的Hough变换检测线段的原理一样，可用于需要进行线检测的相关应用中，如线轨迹检测、滤波、倾斜校正等。

采用Radon变换计算出原图中各方向上的投影值，可以作为方向特征用于目标检测和识别，如应用于掌纹、静脉识别。

Radon变换改变图像的表现形式，为相关处理提供便利，如图像复原中在Radon域用高阶统计量估计点扩散函数，提高了算法的运算速度。

【例518】对一幅图像进行Radon正变换，检测其中的线段。
程序如下： 

clear,clc,close all;

Image=rgb2gray(imread('rail.jpg'));

theta=0:10:180;

［R1,X1］=radon(Image,theta);%Radon变换

［N,M］=size(R1);%矩阵R1的尺寸，M为theta中元素的个数

R2=reshape(R1,1,N*M);%变换R1的形式，为找M个最大值做准备

sortR=sort(R2,'descend');%降序排列

R2=R1;

R2(R1<sortR(M))=0;%只保留最大的M个值，其余置为0 

R3=R1;

R3(R1<max(R1(:)))=0;%只保留最大的一个值，其余置为0

result1=iradon(R2,theta);

result2=iradon(R3,theta);%两种情况下的Radon变换

subplot(221),imshow(Image),title('原图');

colormap(gray);

subplot(222),imagesc(theta,X1,R1),title('间隔10°投影的Radon变换曲线');

subplot(223),image(result1),title('根据theta个数保留最大值对应线段');

subplot(224),image(result2),title('仅保留R最大值对应线段');

程序运行结果如图519、图520所示。


图519Radon变换检测线段示例一




图520Radon变换检测线段示例二

5.5小波变换
作为重要的数学工具，小波变换被应用到数字图像处理的多个方面，如图像平滑、边缘检测、图像分析及压缩编码等。本节介绍小波变换的基本原理、特性及其在图像处理中的应用。
5.5.1小波
波(Wave)被定义为时间或空间的一个振荡函数，如一条正弦曲线。小波(Wavelet)是“小的波”，具有在时间上集中能量的能力，是分析瞬变的、非平稳的或时变现象的工具。如图521所示，正弦曲线在-∞≤t≤∞上等振幅振荡，具有无限能量，而小波具有围绕一点集结的有限能量。


图521波和小波


1. 定义

定义1： 设函数ψ(t)满足下列条件： 


∫Rψ(t)dt=0(550)


对其进行平移和伸缩产生如下函数族ψa，b(t)。


ψa，b(t)=1aψt-baa，b∈R，a≠0(551)


ψ(t)称为基小波或母小波，a称为伸缩因子(尺度因子)，b为平移因子，ψa，b(t)称为ψ(t)生成的连续小波。由傅里叶变换性质可得： 


Ψa，b(ω)=aΨ(aω)e-jωb(552)


定义2： 若函数ψ(t)的傅里叶变换Ψ(ω)满足： 


Cψ=∫R|Ψ(ω)|2|ω|dω<∞(553)


则称ψ(t)为允许小波，式(553)称为允许性条件。其中Ψ(ω)=∫Rψ(t)e-jωtdt。
因Ψ(ω)|ω=0=∫Rψ(t)dt=0，故允许小波一定是基小波。
2. 特点
1) 紧支撑性

小波函数ψ(t)满足式(550)，即均值为0，ψ(t)应具有振荡性，即在图形上具有“波”的形状。ψ(t)满足
∫R|ψ(t)|dt<∞、∫R|ψ(t)|2dt<∞，因此，ψ(t)仅在小范围内波动，且能量有限，即小波函数ψ(t)的定义域是紧支撑的，超出一定范围时，波动幅度迅速衰减，具有速降性。
2) 变化性
小波函数ψa，b(t)及它的频谱Ψa，b(ω)随尺度因子a的变化而变化。由式(551)可知，随着a的减小，ψa，b(t)的支撑区随之变窄，其幅值变大，如图522所示。

由傅里叶变换的尺度变换性质： 若F［f(t)］=F(ω)，则F［f(αt)］=F(ω/α)/|α|，可知，Ψa，b(ω)随着a的减小而向高频端展宽。


图522Marr小波参数a，b取不同值的波形


3) 消失矩
若小波ψ(t)满足式(554)，则称该小波具有K阶消失矩。


∫Rtkψ(t)dt=0k=0，1，…，K-1(554)


这时，Ψ(ω)在ω=0处K次可微，即Ψk(0)=0，k=1，2，…，K，小波ψ(t)随着K的增加，波形振荡越来越强烈。

3. 一维小波实例
1) Haar小波
Haar小波是最简单的小波，其表达式为： 


ψH(t)=10≤t<1/2

-11/2≤t<1

0其他

ΨH(ω)=1-2e-iω2+e-iωωi(555)


Haar小波具有紧支性(长度为1)和对称性，消失矩为1，即
∫RψH(t)dt=0。Haar小波不是连续函数，应用有限，但结构简单，一般作为原理示意或说明。

2) Morlet小波
Morlet小波是用高斯函数构造的一种小波，其时域、频域可表示为： 



ψ(t)=π-1/4(e-iω0t-e-ω20/2)e-t2/2

Ψ(ω)=π-1/4［e-(ω-ω0)2/2-e-ω20/2e-ω2/2］
(556)


由上式可以看出，Morlet小波满足允许条件，即
Ψ(0)=0。
当ω0≥5时，e-ω20/2≈0，所以，式(556)的第2项可以忽略，Morlet小波可以近似表示为： 



ψ(t)=π-1/4e-iω0te-t2/2


Ψ(ω)=π-1/4e-(ω-ω0)2/2
(557)



Morlet小波在时域、频域都具有较好的局部性，是很常用的小波。
3) Mexican草帽小波

Mexican草帽小波成比例于高斯函数二阶导数，也称为Marr小波，其表达式为： 



ψ(t)=23π-1/4(1-t2)e-t2/2(558)



Mexican草帽小波具有对称性，支撑区间是无限的，有效支撑区间为［-5 5］，在视觉信息处理方面获得了很多应用。

MATLAB提供了相应函数生成一维小波，其调用格式如下。
(1) ［PSI，X］=morlet(LB，UB，N)： 返回Morlet小波在区间［LB UB］内N个规则网格点的值PSI。小波有效紧支集为［-4 4］，函数返回的PSI是按照cos(5t)e-t2/2计算的，即式(557)中ω0=5，只取了实部，且未考虑系数π-1/4。

(2) ［PSI，X］=mexihat(LB，UB，N)： 返回Mexican草帽小波在区间［LB UB］内N个规则网格点的值PSI。小波有效紧支集为［-5 5］。
【例519】绘制Haar小波、Morlet小波及Marr小波的函数图形及其频谱图。
程序如下： 

clear,clc,close all;

t=0:0.01:1;

N=length(t); 

Hpsi(1:N/2)=1;Hpsi(N/2:N)=-1; %Haar小波

HDFT=abs(fftshift(fft(Hpsi)))*2/N;

x=linspace(0,1,N);

subplot(121),plot(x,Hpsi,'k',［0,1.5］,［0,0］,':k',［1,1］,［-1,0］,':k'),…


ylim(［-2 2］),title('Haar小波');


subplot(122),plot(x,HDFT,'k'),title('Haar小波频谱');

t=-4:0.01:4;

omega0=5;

Realpsi=pi^(-0.25)*(cos(omega0*t)-exp(-omega0^2/2)).*exp(-t.^2/2); %Morlet小波实部

Impsi=pi^(-0.25)*(-sin(omega0*t)-exp(-omega0^2/2)).*exp(-t.^2/2);  %Morlet小波虚部

figure;

subplot(121),plot(t,Realpsi,'k',t,Impsi,'k:'),title('Morlet小波'); 

omega=0:0.01:10;

Fm=pi^(-0.25)*(exp(-(omega-omega0).^2/2)-exp(-omega0^2/2-omega.^2/2));

subplot(122),plot(omega,Fm,'k'),title('Morlet小波频谱');

［Mpsi,t］=mexihat(-5,5,256); %Mexican小波

figure;subplot(121),plot(t,Mpsi,'k'),title('Mexican小波');

MDFT=abs(fftshift(fft(Mpsi)));

x=linspace(-5,5,length(MDFT));

subplot(122),plot(x,MDFT,'k'),title('Mexican小波频谱');

程序运行结果如图523所示。


图5233种常见小波及其频谱


5.5.2一维小波变换

将ψa，b(t)中的参数a，b离散化，取a=aj0，j=0，±1，±2，…，b=kaj0，j，k∈Z，则离散化后的小波函数为： 



ψj，k(t)=a-j/20ψ(a-j0t-k)，j，k∈Z(559)


用ψj，k(t)将函数f(t)展开，即


f(t)=∑j∑kαj，kψj，k(t)(560)


展开系数αj，k的集合称为f(t)的离散小波变换(Discrete Wavelet Transform，DWT)，如式（561）所示。


αi，k=〈f(t)，ψi，k(t)〉=∫Rf(t)ψ*j，k(t)dt(561)


1. 正交小波变换与多分辨分析

定义1： 设平方可积函数空间L2(R)中的函数ψ(t)是一个允许小波，若其二进伸缩平移系(a0=2)
构成L2(R)的标准正交基，则称ψ(t)为正交小波，称ψj，k(t)是正交小波函数，称相应的离散小波变换αj，k=〈f(t)，ψj，k(t)〉为正交小波变换，如式（562）所示。


ψj，k(t)=2-j/2ψ(2-jt-k)，j，k∈Z(562)


例如，Haar小波的二进伸缩平移系如下： 


ψj，k(t)=2-j/2ψ(2-jt-k)=2-j/22jk≤t<(2k+1)2j-1

-2-j/2(2k+1)2j-1≤t≤(k+1)2j

0其他(563)



可以验证构成L2(R)的一个标准正交基。

多分辨分析(MultiResolution Analysis，MRA)也称多尺度分析，是构造正交小波基的一般方法。

定义2： 若L2(R)中一个子空间序列
{Vj}j∈Z及一个函数φ(t)满足以下条件，则称其为一个正交多分辨分析。
(1) VjVj-1，j∈Z； 
(2) f(t)∈Vjf(2t)∈Vj-1； 
(3) ∩j∈ZVj=0，
∪j∈ZVj=L2(R
)； 
(4) φ(t)∈V0，且{φ(t-k)}k∈Z是V0的标准正交基，称φ(t)是此多分辨分析的尺度函数。
2. 尺度函数与小波函数
由多分辨分析性质(2)、(4)可知，对于任何φ(t)∈V0，有φ(2-jt)∈Vj； 
{φ(t-k)}k∈Z是V0的标准正交基，函数系{2-j/2φ(2-jt-k)}k∈Z则构成了Vj的一组标准正交基； 即{φ(t-k)}k∈Z张成L2(R)的子空间V0，{2-j/2φ(2-jt-k)}k∈Z张成了Vj。

φ(t)∈V0，而V0V-1，因此，φ(t)∈V-1。因函数系{21/2φ(2t-k)}k∈Z构成了V-1的一组标准正交基，因此，φ(t)可以借助于{21/2φ(2t-k)}k∈Z的加权和表示。


φ(t)=∑k∈Zhk2φ(2t-k)(564)


其中，系数{hk}k∈Z称为尺度函数(尺度滤波器)系数，满足



hk=12∫Rφ(t)φ*(2t-k)dt

H(ω)=12∑khke-ikω(565)


式(564)称为双尺度方程，其频域形式为： 


Φ(2ω)=H(ω)Φ(ω)(566)


定义函数ψj，k(t)，张成尺度函数在不同尺度下张成的空间之间的差空间为{Wj}j∈Z，称Wj为尺度为j的小波空间，Vj为尺度为j的尺度空间。

由于VjVj-1，即Vj-1=Vj+Wj，且Wj⊥Vj，j∈Z，显然，当m，n∈Z，m≠n时，有Wm⊥Wn，
如式（567）所示。


Vj-1=Vj+Wj=Vj+1+Wj+1+Wj=…

=Vj+s+Wj+s+Wj+s-1+…+Wj+1+Wj

(567)


令s→+∞，则Vj-1=+∞m=jWm
令j→-∞，则L2(R)=+∞m=-∞Wm

因W0V-1，张成W0的小波函数ψ(t)可以由V-1的标准正交基{21/2φ(2t-k)}k∈Z表示。


ψ(t)=∑k∈Zgk2φ(2t-k)(568)


式(568)也称为双尺度方程，其频域表示为： 


Ψ(2ω)=G(ω)Φ(ω)(569)


{gk}k∈Z满足： 



gk=12∫Rψ(t)φ*(2t-k)dt

G(ω)=12∑kgke-ikω
(570)


式(564)和式(568)这两个双尺度方程是多分辨分析赋予尺度函数φ(t)和小波函数ψ(t)的最基本性质。

综上所述，多分辨分析的基本思想其实就是： 为有效地寻找空间L2(R)的基底，从L2(R)的某个子空间出发，在这个子空间中建立基底，然后利用简单的变换把该基底扩充到L2(R)中去。

3. 常见小波函数的尺度函数计算
MATLAB提供了wavefun函数，用于计算小波函数的近似值及对应的尺度函数； 提供了waveinfo函数查看小波的相关信息，其部分调用如下。
(1) ［PHI，PSI，XVAL］=wavefun('wname'，ITER)： 对正交小波、Meyer小波适用，返回尺度函数PHI和小波函数PSI的取值； ITER为迭代次数； XVAL为取值点，支撑区间内有2ITER个点。
(2) ［PSI，XVAL］=wavefun('wname'，ITER)： 对没有尺度函数的小波适用，如Morlet小波、Mexican草帽小波、高斯导数小波或复小波。输出PSI为实数或复数向量。
(3) …=wavefun(…，'plot')： 计算并绘制函数。
(4)  waveinfo('wname')： 提供wname指定的小波的信息，包括haar、db、sym、coif、bior、rbio、meyr、dmey、gaus、mexh、morl、cgau、cmor、shan、fbsp、fk等。
【例520】使用wavefun函数绘制Haar小波、Morlet小波、Meyer小波、Daubechies(dbN)小波的尺度函数及小波函数图形。
程序如下： 

clear,clc,close all;

figure,［phi1,psi1,xval1］=wavefun('haar',8,'plot');

figure,［psi2,xval2］=wavefun('morl',8,'plot');

figure,［phi3,psi3,xval3］=wavefun('meyr',8,'plot');

figure,［phi4,psi4,xval4］=wavefun('db4',8,'plot');

程序运行结果如图524所示。


图524常见小波的尺度函数和小波函数


4. 函数的正交小波分解

由以上讨论可知，给定一个多分辨分析({Vk}k∈Z，φ(t))，可确定一个小波函数ψ(t)和其伸缩系{ψj，k(t)=2-j/2ψ(2-jt-k)}j，k∈Z，并张成小波空间{Wj}j∈Z，因Wi⊥Wj(i≠j)，且L2(R)=j∈ZWj，所以{ψj，k(t)=2-j/2ψ(2-jt-k)}j，k∈Z构成L2(R)的标准正交基。因此，对任何f(t)∈L2(R)，有： 


f(t)=∑j，kdj，kψj，k(t)
(571)


其中，dj，k=〈f(t)，ψj，k(t)〉j，k∈Z是f(t)的离散小波变换，且是正交小波变换，式(571)是f(t)的重构公式，也称为f(t)的正交小波分解。

由多分辨分析性质(2)，即f(t)∈Vjf(2t)∈Vj-1可知： Vj的频率范围是Vj-1的一半，且是Vj-1中的低频表现部分，而Vj-1=Vj+Wj，所以，Wj的频率表现在Vj与Vj-1之间的部分，而且Wj的频带互不重叠，因此，通常认为Vj表现了Vj-1的“概貌”，Wj表现了Vj-1的不同频带中的“细节”。记dj，kψj，k(t)=wj(t)，则wj(t)∈Wj，式(571)可写成： 


f(t)=∑jwj(t)(572)


式(572)说明任何一个函数f(t)∈L2(R)可以分解成不同频带的细节之和。

实际情况中，函数f(t)∈L2(R)仅有有限的细节。由式(567)得V0=Vs+Ws+Ws-1+…+W1。设一个函数f(t)∈V0，则f(t)可分解为f(t)=fs(t)+ws(t)+ws-1(t)+…+w1(t)，即


f(t)=∑k∈Zcs，kφs，k(t)+∑sj=1∑k∈Zdj，kψj，k(t)(573)


其中，cs，k=〈f(t)，φs，k(t)〉，k∈Z，dj，k=〈f(t)，ψj，k(t)〉，k∈Z。称式中第1项为f(t)的不同尺度s(s≥1)下的近似式，是f(t)中频率不超过2-s的成分； 称第2项中的∑k∈Zdj，kψj，k(t)为f(t)的不同尺度j下的细节，是f(t)中频率2-j到2-j+1之间的细节成分。

根据以上分析可知，当尺度函数φ(t)、小波函数ψ(t)确定后，通过计算{cs，k}k∈Z、{dj，k}k∈Z，即可得到函数f(t)∈L2(R)的近似和细节。



cj+1，k=〈f(t)，φj+1，k(t)〉=∫Rf(t)φ*j+1，k(t)dt=∫Rf(t)2-(j+1)/2φ*(2-(j+1)t-k)dt

=∫Rf(t)2-(j+1)/2∑n∈Zh*n2φ*［2(2-(j+1)t-k)-n］dt

=∫Rf(t)2-j/2∑n∈Zh*nφ*［2-jt-(2k+n)］dt

=∑m∈Zh*m-2k∫Rf(t)2-j/2φ*(2-jt-m)dt

=∑m∈Zh*m-2k∫Rf(t)φ*j，m(t)dt

=∑m∈Zh*m-2k〈f(t)，φj，m(t)〉

=∑m∈Zh*m-2kcj，m



同理，得到dj+1，k=〈f(t)，ψj+1，k(t)〉=∑m∈Zg*m-2kcj，m
因此，得到正交小波分解的Mallat快速算法，如式（574）所示。



cj+1，k=∑n∈Zh*n-2kcj，n

dj+1，k=∑n∈Zg*n-2kcj，nk∈Z(574)


因此，只要知道双尺度方程中的传递系数{hk}k∈Z(gk=(-1)kh*1-k)，就可计算出一系列正交小波分解系数，过程如图525所示。


图525正交小波分解算法示意图


初始值c0，k=〈f(t)，φ0，k(t)〉=∫Rf(t)φ*(t-k)dt，对于离散序列，f(t)→fn=f(nΔt)，因此，



c0，k≈∑nfnφ(n-k)(575)


根据信号处理理论，利用序列{hk}k∈Z对一个离散信号{xn}n∈Z∈l2(Z)进行滤波，则


yk=hk*xk=∑n∈Zhk-nxn(576)


比较式(574)和式(576)发现，式(576)卷积式中k对所有的n值做卷积运算，而式(574)卷积式中是2k对所有的n值做卷积运算，缺少了奇数(2k+1)的部分，即卷积运算或滤波处理之后所得的序列抽去了k的奇数部分，只剩下偶数部分，这一过程称为再抽样，抽样率为2。所以，分辨率j的近似分量cj，k分解为分辨率为j+1的近似分量cj+1，k和细节分量dj+1，k的分解方法可以用图526所示的滤波过程来表示。


图526近似分量cj，k分解为cj+1，k和dj+1，k(2↓代表再抽样，抽样率为2)


5. 函数的正交小波重构

所谓重构，即已知近似序列{cj+1，k}k∈Z和细节序列{dj+1，k}k∈Z，求出序列{cj，k}k∈Z。
由正交小波分解式可知： 



fj(t)=∑k∈Z〈f(t)，φj，k(t)〉φj，k(t)=∑k∈Zcj，kφj，k(t)(577)


由于Vj=Vj+1+Wj+1，所以，fj(t)=fj+1(t)+wj+1(t)，而


fj+1(t)=∑k∈Z〈f(t)，φj+1，k(t)〉φj+1，k(t)=∑k∈Zcj+1，kφj+1，k(t)

wj+1(t)=∑k∈Z〈f(t)，ψj+1，k(t)〉ψj+1，k(t)=∑k∈Zdj+1，kψj+1，k(t)

fj+1(t)+wj+1(t)=∑k∈Zcj+1，kφj+1，k(t)+∑k∈Zdj+1，kψj+1，k(t)

=∑k∈Zcj+1，k2-(j+1)/2φ(2-(j+1)t-k)+∑k∈Zdj+1，k2-(j+1)/2ψ(2-(j+1)t-k)

=∑k∈Zcj+1，k2-(j+1)/2∑n∈Zhn2φ［2(2-(j+1)t-k)-n］+

∑k∈Zdj+1，k2-(j+1)/2∑n∈Zgn2φ［2(2-(j+1)t-k)-n］

=∑k∈Zcj+1，k2-j/2∑n∈Zhnφ［2-jt-(2k+n)］+

∑k∈Zdj+1，k2-j/2∑n∈Zgnφ［2-jt-(2k+n)］

=∑k∈Zcj+1，k2-j/2∑m∈Zhm-2kφ［2-jt-m］+

∑k∈Zdj+1，k2-j/2∑m∈Zgm-2kφ［2-jt-m］

=∑k∈Zcj+1，k∑m∈Zhm-2kφj，m(t)+∑k∈Zdj+1，k∑m∈Zgm-2kφj，m(t)

=∑m∈Z(∑k∈Zcj+1，khm-2k+∑k∈Zdj+1，kgm-2k)φj，m(t)(578)


由式(577)和式(578)可得Mallat小波重构算法。


cj，k=∑k∈Zcj+1，khn-2k+∑k∈Zdj+1，kgn-2k(579)


其重构过程如图527所示。


图527小波重构算法示意图



比较式(579)和式(576)发现，式(576)卷积式中k对所有的n值做卷积运算，而式(579)卷积式中是n对k的偶数序列2k做卷积运算，从而造成cj+1，k、dj+1，k的取值个数比hn-2k、gn-2k的取值个数多出一倍，
可将(2k+1)对应的cj+1，k、dj+1，k当作0值来处理，即在两个数值之间插入一个0，这一过程称为插值抽样，抽样率为2。所以，分辨率j+1的近似分量cj+1，k和细节分量dj+1，k重构分辨率j级近似分量cj，k的重构方法可以用图528所示的滤波过程来表示。


图528cj+1，k和dj+1，k重构cj，k(2↑代表插值抽样，抽样率为2)


6. 正交小波变换的实现

MATLAB提供了一系列的函数实现小波分解、重构等功能。
(1) L=wmaxlev(S，'wname')： 计算尺寸为S的信号使用wname指定的小波分解时的最大分解级数。
(2) wfilters函数： 计算小波对应的滤波器。
① ［LO_D，HI_D，LO_R，HI_R］=wfilters('wname')： 计算wname指定的正交或双正交小波对应的分解及重构滤波器。
② ［F1，F2］=wfilters('wname'，'type')： 根据type的不同，返回不同的滤波器。type为d，返回分解滤波器LO_D和HI_D； type为r，返回重构滤波器LO_R和HI_R； type为l，返回低通滤波器LO_D 和 LO_R； type为h，返回高通滤波器HI_D和HI_R。
（3） dwt函数： 一级一维离散小波变换。
① ［CA，CD］=dwt(X，'wname')： 对向量X进行wname指定的小波分解，计算近似系数向量CA和细节系数向量CD。
② ［CA，CD］=dwt(X，Lo_D，Hi_D)： 使用给定的分解滤波器实现向量X的小波分解。
(4) idwt函数： 一级一维小波逆变换。
① X=idwt(CA，CD，'wname')： 基于近似系数向量CA和细节系数向量CD，使用wname指定的小波实现近似系数向量X的一级重构。
② X=idwt(CA，CD，Lo_R，Hi_R)： 给定重构滤波器实现一级近似系数向量X的计算。

（5） wavedec函数： 实现多级一维小波分解。
① ［C，L］=wavedec(X，N，'wname')： 使用wname指定的小波实现信号X在N级上的小波分解。N为正整数，可使用wmaxlev函数计算，确保小波系数不受边界效应的影响； 若不考虑边界效应，可设置N≤fix(log2(length(X)))。输出向量C为小波的分解系数，L为各级系数的长度。
② ［C，L］=wavedec(X，N，Lo_D，Hi_D)： 使用指定的低通分解滤波器Lo_D和高通分解滤波器Hi_D实现信号X在N级上的小波分解。
(6) detcoef函数： 提取一维小波分解的细节系数。
① D=detcoef(C，L，N)： 从小波分解结构［C，L］中提取第N级细节系数。N为整数，N≥1且N≤length(L)-2。如果N为整数向量，则提取各元素所示级别的细节系数。
② D=detcoef(C，L)： 从小波分解结构［C，L］中提取最高级细节系数。
(7) appcoef函数： 提取一维小波分解的近似系数。
① A=appcoef(C，L，'wname'，N)： 从小波分解结构［C，L］中提取第N级近似系数。N为整数，N≥1且N≤length(L)-2。

② A=appcoef(C，L，'wname')： 提取第length(L)-2级近似系数。
（8） wrcoef函数： 由一维小波系数进行单支重构。
① X=wrcoef('type'，C，L，'wname'，N)： 基于小波分解结构［C，L］在N级计算重构系数向量。type取a，重构近似向量； 取b，重构细节向量； N为整数，且N≤length(L)-2。
② X=wrcoef('type'，C，L，Lo_R，Hi_R，N)： 基于重构低通滤波器Lo_R和重构高通滤波器Hi_R实现重构系数向量计算。
（9） waverec函数： 多级一维小波重构。
① X=waverec(C，L，'wname')： 基于多级小波分解结构［C，L］ 重构信号X。
② X=waverec(C，L，Lo_R，Hi_R)： 给定重构滤波器实现重构系数向量计算。
(10) upcoef函数： 一维小波系数直接重构。
① Y=upcoef(O，X，'wname'，N)： 计算向量X向上N步的重构系数。N为正整数。若O取'a'，近似系数被重构； 若O取'd'，细节系数被重构。
② Y=upcoef(O，X，'wname'，N，L)： 重构的同时取出结果中长度为L的中间部分。
③ Y=upcoef(O，X，Lo_R，Hi_R，N)或Y=upcoef(O，X，Lo_R，Hi_R，N，L)： 使用重构滤波器实现直接重构。
【例521】装载sumsin信号，使用dwt和idwt函数实现一级一维小波分解与重构。sumsin是3个正弦波的叠加，即sumsin(t)=sin(3t)+sin(0.3t)+sin(0.03t)。
程序如下： 

clear,clc,close all;

load sumsin;%装载sumsin信号

signal=sumsin;

［CA,CD］=dwt(signal,'db4');%基于db4小波实现一级一维小波分解

RecS=idwt(CA,CD,'db4');%重构

subplot(221),plot(signal),title('原始信号');

subplot(222),plot(CA),title('近似系数');

subplot(223),plot(CD),title('细节系数');

subplot(224),plot(RecS),title('重构信号');

程序运行结果如图529所示。


图529一级一维小波分解与重构


【例522】装载sumsin信号，实现多级一维小波分解与重构。
程序如下： 

clear,clc,close all;

load sumsin;

signal=sumsin;

［C,L］=wavedec(signal,2,'db4');%基于db4小波实现2级一维小波分解

［CD1,CD2］=detcoef(C,L,［1 2］);%提取细节系数

CA2=appcoef(C,L,'db4',2);%提取近似系数

［Lo_R,Hi_R］=wfilters('db4','r');%计算重构滤波器系数

Recs1=waverec(C,L,Lo_R,Hi_R);%采用重构滤波器实现信号重构

Recs2= wrcoef('a',C,L,'db4',2);%基于近似系数重构信号

subplot(231),plot(signal),title('原始信号');

subplot(234),plot(CA2),title('二级近似系数');

subplot(232),plot(CD2),title('二级细节系数');

subplot(235),plot(CD1),title('一级细节系数');

subplot(233),plot(Recs2),title('近似系数重构');

subplot(236),plot(Recs1),title('重构信号');

程序运行结果如图530所示。


图530多尺度一维小波分解与重构



小波变换还有很多别的很有用的理论和特点，如小波包、多带小波、多小波等，因篇幅关系，不再深入分析。
5.5.3二维小波变换

图像为二维信号，用二元函数f(x，y)∈L2(R2)表示，可以对其进行二维小波变换和多分辨分析。
1. 二维多分辨分析

设({V1j}j∈Z，φ1(t))，({V2j}j∈Z，φ2(t))是L2(R)的两个多分辨分析，ψ1(t)，ψ2(t)分别是相应的正交小波函数，则V～j是V1j与V2j的张量积空间，如式
（580）所示。


V～j=V1jV2j=f1(x)f2(y)|f1(x)∈V1j，f2(y)∈V2j(580)



{φ1j，l(x)}l∈Z，{φ2j，m(y)}m∈Z是V1j与V2j的标准正交基，则{φ1j，l(x)φ2j，m(y)}l，m∈Z是V～j的标准正交基。
设W1j是V1j在V1j-1中的正交补，W2j是V2j在V2j-1中的正交补，则


V～j-1=V1j-1V2j-1=(V1jW1j)(V2jW2j)

=(V1jV2j)(V1jW2j)(W1jV2j)(W1jW2j)

=V～jW～1jW～2jW～3j(581)



其中，W～1j，W～2j，W～3j称为二维小波空间。它们的标准正交基依次为{φ1j，l(x)ψ2j，m(y)}l，m∈Z，{ψ1j，l(x)φ2j，m(y)}l，m∈Z和{ψ1j，l(x)ψ2j，m(y)}l，m∈Z。
记


ψ1(x，y)=φ1(x)ψ2(y)

ψ2(x，y)=ψ1(x)φ2(y)

ψ3(x，y)=ψ1(x)ψ2(y)

φ(x，y)=φ1(x)φ2(y)(582)


则φ(x，y)，ψ1(x，y)，ψ2(x，y)，ψ3(x，y)的伸缩平移系分别构成V～j，W～1j，W～2j，W～3j的标准正交基。
由式(581)可知： 


L2(R2)=∑+∞j=-∞W～j(583)


其中，W～j=W～1jW～2jW～3j。对于任何f(x，y)∈L2(R2)，有： 



f(x，y)=∑+∞j=-∞wj(x，y)(584)


其中，wj(x，y)∈W～j。
因此，二维小波变换的重构公式为： 


f(x，y)=∑+∞j=-∞∑l，mαjl，mφ1j，l(x)ψ2j，m(y)+βjl，mψ1j，l(x)φ2j，m(y)

+γjl，mψ1j，l(x)ψ2j，m(y)(585)


其中，αjl，m、βjl，m、γjl，m是f(x，y)的二维离散小波变换，
αjl，m=R2f(x，y)
φ1*j，l(x)ψ2*j，m(y)
dxdy，βjl，m=R2f(x，y)
ψ1*j，l(x)φ2*j，m(y)dxdy，γjl，m=R2f(x，y)ψ1*j，l(x)ψ2*j，m(y)dxdy。
实际问题中，二元函数f(x，y)只有有限分辨率，设f(x，y)∈V～0，因此，


f(x，y)=fs(x，y)+ws(x，y)+ws-1(x，y)+…+w1(x，y)(586)

f(x，y)=∑l，m［λsl，mφ1s，l(x)φ2s，m(y)］+

∑sj=1∑l，m［αjl，mφ1j，l(x)ψ2j，m(y)+βjl，mψ1j，l(x)φ2j，m(y)+γjl，mψ1j，l(x)ψ2j，m(y)］

(587)


其中，λsl，m=R2f(x，y)φ1*s，l(x)φ2*s，m(y)dxdy。
式(587)中第一项fs(x，y)是f(x，y)在尺度s下的近似； 后三项称为f(x，y)在不同尺度j下的细节。

2. 二维正交小波分解

由于φ1(x)，φ2(y)，ψ1(x)，ψ2(y)满足双尺度方程： 



φ1(x)=2∑lh1lφ1(2x-l)

ψ1(x)=2∑lg1lφ1(2x-l)

φ2(y)=2∑mh2mφ2(2y-m)

ψ2(y)=2∑mg2mφ2(2y-m)
(588)


将式(588)代入式(582)，得



φ(x，y)=φ1(x)φ2(y)=2∑l，mh1lh2mφ1(2x-l)φ2(2y-m)=2∑l，mh1lh2mφ(2x-l，2y-m)

ψ1(x，y)=φ1(x)ψ2(y)=2∑l，mh1lg2mφ1(2x-l)φ2(2y-m)=2∑l，mh1lg2mφ(2x-l，2y-m)

ψ2(x，y)=ψ1(x)φ2(y)=2∑l，mg1lh2mφ1(2x-l)φ2(2y-m)=2∑l，mg1lh2mφ(2x-l，2y-m)

ψ3(x，y)=ψ1(x)ψ2(y)=2∑l，mg1lg2mφ1(2x-l)φ2(2y-m)=2∑l，mg1lg2mφ(2x-l，2y-m)
(589)



则


λj+1l，m=R2f(x，y)φ1*j+1，l(x)φ2*j+1，m(y)dxdy

=R2f(x，y)2-(j+1)/2φ1*［2-(j+1)x-l］2-(j+1)/2φ2*［2-(j+1)y-m］dxdy

=R2f(x，y)2-j/2∑nh1*nφ1*［2(2-(j+1)x-l)-n］

2-j/2∑kh2*kφ2*［2(2-(j+1)y-m)-k］dxdy

=R2f(x，y)2-j/2∑ph1*p-2lφ1*(2-jx-p)2-j/2∑qh2*q-2mφ2*(2-jy-q)dxdy

=∑p，qh1*p-2lh2*q-2mR2f(x，y)φ1*j，p(x)φ2*j，q(y)dxdy

=∑p，qh1*p-2lh2*q-2mλjp，q
(590)


同理，得



αj+1l，m=∑p，qh1*p-2lg2*q-2mλjp，q

βj+1l，m=∑p，qg1*p-2lh2*q-2mλjp，q

γj+1l，m=∑p，qg1*p-2lg2*q-2mλjp，q(591)


由式(590)和式(591)可以看出，分辨率j的近似分量λjp，q分解为分辨率为j+1的近似分量λj+1l，m和细节分量αj+1l，m，βj+1l，m，γj+1l，m的分解方法可以用
如图531所示的滤波过程来表示： 首先对水平方向进行滤波，再对垂直方向进行滤波，得到4个不同的频带； 若对近似分量λj+1l，m继续进行这样的滤波过程，即可得到如图532所示的塔形分解。


图531二维小波变换近似分量λjp，q分解为λj+1l，m和αj+1l，m、βj+1l，m、γj+1l，m


若对一幅二维图像进行3层分解，可得图532。其中，L代表低频分量，H代表高频分量； LH代表垂直方向上的高频信息； HL频带存放的是图像水平方向的高频信息； HH频带存放图像在对角线方向的高频信息。



图532二维图像3层小波分解示意图


3. 二维正交小波重构
因为V～j=V～j+1W～1j+1W～2j+1W～3j+1，所以，


fj(x，y)=fj+1(x，y)+w1j+1(x，y)+w2j+1(x，y)+w3j+1(x，y)


而


fj+1(x，y)=∑l，m∈Z〈f(x，y)，φj+1，l，m(x，y)〉φj+1，l，m(x，y)=∑l，m∈Zλj+1l，mφ1j+1，l(x)φ2j+1，m(y)

=∑l，m∈Zλj+1l，m［2-(j+1)/2φ1(2-(j+1)x-l)］［2-(j+1)/2φ2(2-(j+1)y-m)］

=∑l，m∈Zλj+1l，m
2-j/2∑n∈Zh1nφ1［2(2-(j+1)x-l)-n］


2-j/2∑n∈Zh2nφ2［2(2-(j+1)y-m)-n］

=∑l，m∈Zλj+1l，m2-j/2∑p∈Zh1p-2lφ1{2-jx-p}2-j/2∑q∈Zh2q-2mφ2(2-jy-q)

=∑l，m∈Zλj+1l，m∑p，q∈Zh1p-2lh2q-2mφ1j，p(x)φ2j，q(y)(592)


同理，得


w1j+1(x，y)=∑l，m∈Z〈f(x，y)，ψ1j+1，l，m(x，y)〉ψ1j+1，l，m(x，y)=∑l，m∈Zαj+1l，mφ1j+1，l(x)ψ2j+1，m(y)

=∑l，m∈Zαj+1l，m∑p，q∈Zh1p-2lg2q-2mφ1j，p(x)φ2j，q(y)(593)

w2j+1(x，y)=∑l，m∈Z〈f(x，y)，ψ2j+1，l，m(x，y)〉ψ2j+1，l，m(x，y)=∑l，m∈Zβj+1l，mψ1j+1，l(x)φ2j+1，m(y)

=∑l，m∈Zβj+1l，m∑p，q∈Zg1p-2lh2q-2mφ1j，p(x)φ2j，q(y)(594)

w3j+1(x，y)=∑l，m∈Z〈f(x，y)，ψ3j+1，l，m(x，y)〉ψ3j+1，l，m(x，y)=∑l，m∈Zγj+1l，mψ1j+1，l(x)ψ2j+1，m(y)

=∑l，m∈Zγj+1l，m∑p，q∈Zg1p-2lg2q-2mφ1j，p(x)φ2j，q(y)(595)


所以，


fj+1(x，y)+w1j+1(x，y)+w2j+1(x，y)+w3j+1(x，y)

=∑p，q∈Z∑l，m∈Z［λj+1l，mh1p-2lh2q-2m+αj+1l，mh1p-2lg2q-2m+βj+1l，mg1p-2lh2q-2m+γj+1l，mg1p-2lg2q-2m］

φ1j，p(x)φ2j，q(y)


而


fj(x，y)=∑p，q∈Z〈f(x，y)，φj，p，q(x，y)〉φj，p，q(x，y)=∑p，q∈Zλjp，qφ1j，p(x)φ2j，q(y)


因此


λjp，q=∑l，m∈Z［λj+1l，mh1p-2lh2q-2m+αj+1l，mh1p-2lg2q-2m+βj+1l，mg1p-2lh2q-2m+γj+1l，mg1p-2lg2q-2m］(596)


式(596)所示的重构可以用图533所示的滤波过程来表示。


图533二维多分辨分析的重构



4. 二维正交小波的实现
利用MATLAB小波工具箱可以实现对图像的小波分解与重构。在命令窗口输入： 

>>wavemenu

打开小波分析工具箱，如图534(a)所示； 选择Wavelet 2D，如图534(b)所示。

在二维小波分析页面，从菜单中打开图像，在右侧选择Daubechies小波(N=4)进行分解与重构，如图535所示。可以尝试其他的分析和统计计算。


MATLAB提供了相应的二维小波函数，部分列举如下。
(1) wavefun2函数： 计算尺度、小波函数。
① ［S，W1，W2，W3，XYVAL］=wavefun2('wname'，ITER)： 返回wname指定的正交小波的尺度函数S与3个小波函数W1、W2、W3，ITER为迭代次数，XYVAL为取值点，支撑区间内有2ITER×2ITER个点。
② …=wavefun2(…，'plot')： 计算并绘制函数。
(2) dwt2函数： 实现一级二维离散小波变换。
① ［CA，CH，CV，CD］=dwt2(X，'wname')： 用wname指定的小波分解矩阵X为近似系数矩阵CA和细节系数矩阵CH、CV和CD。




图534MATLAB中的小波分析工具





图535lotus图像的小波分解与重构


② ［CA，CH，CV，CD］=dwt2(X，Lo_D，Hi_D)： 采用低通分解滤波器Lo_D和高通分解滤波器Hi_D分解矩阵X。
（3） idwt2函数： 一级二维离散小波逆变换。
① X=idwt2(CA，CH，CV，CD，'wname')： 基于近似系数矩阵CA和细节系数矩阵CH、CV和CD，使用wname指定的小波实现近似系数矩阵X的一级重构。
② X=idwt2(CA，CH，CV，CD，Lo_R，Hi_R)： 采用低通重构滤波器Lo_R和高通重构滤波器Hi_R重构矩阵X。
(4) wavedec2函数： 多级二维小波分解。
① ［C，S］=wavedec2(X，N，'wname')： 使用wname指定的小波实现矩阵X的N级分解。
② ［C，S］=wavedec2(X，N，Lo_D，Hi_D)： 使用分解滤波器分解矩阵X。

C=［A(N)|H(N)|V(N)|D(N)|H(N-1)|V(N-1)|D(N-1)|…|H(1)|V(1)|D(1)］，A、H、V、D分别为低频、水平高频、垂直高频、对角高频系数矩阵。N为正整数。S(1，:)是级数N的低频系数长度； S(i，:)是级数N-i+2的高频系数长度，i=2，…，N+1； S(N+2，:)=size(X)。
(5) waverec2函数： 多级二维小波重构。
① X=waverec2(C，S，'wname')： 基于多级小波分解结构［C，S］重构矩阵X。
② X=waverec2(C，S，Lo_R，Hi_R)： 采用重构滤波器实现矩阵X重构。
（6） appcoef2函数： 提取二维小波分解的低频系数。
① A=appcoef2(C，S，'wname'，N)： 从小波分解结构［C，L］中提取第N级近似系数矩阵A，N为正整数，且N≤size(S，1)-2。
② A=appcoef2(C，S，'wname')： 提取第size(S，1)-2级近似系数矩阵A。
③ A=appcoef2(C，S，Lo_R，Hi_R)或者A=appcoef2(C，S，Lo_R，Hi_R，N)： 基于分解滤波器提取近似系数A。
(7) detcoef2函数： 提取二维小波分解的高频系数。
① D=detcoef2(O，C，S，N)： 从小波分解结构［C，L］中提取第N级细节系数矩阵D。O可取'h' 、'v'或者'd'，对应水平、垂直和对角高频系数； N≥1且 N≤size(S，1)-2。
② ［H，V，D］=detcoef2('all'，C，S，N)： 提取第N级的水平H、垂直V、对角D细节系数矩阵。
（8） upcoef2函数： 二维小波系数的直接重构。
① Y=upcoef2(O，X，'wname'，N，S)： 计算矩阵X的向上N步的重构系数，并提取长度为S的中间部分。若O取'a'，则对近似系数重构； 若O取'h'、'v'、'd'，则对水平、垂直或对角细节系数重构。N为正整数。
② Y=upcoef2(O，X，Lo_R，Hi_R，N，S)： 基于重构滤波器实现重构。
（9） Y=wcodemat(X，NBCODES，OPT，ABSOL)： 返回对数据矩阵X编码的矩阵Y。NBCODES为伪编码的最大值，即编码范围为0~NBCODES-1，默认为16； 若ABSOL为0，返回编码矩阵； 若ABSOL非0，返回数据矩阵的绝对值abs(X)； 默认为1。OPT指定了编码方式，若为'row'或'r'，按行编码； 若为'col'或'c'，按列编码； 若为'mat' 或 'm'，按整个矩阵编码； 默认为'mat'。
【例523】对cameraman图像进行一级小波分解及重构。
程序如下： 

clear,clc,close all;

Image=imread('cameraman.jpg');

grayI=rgb2gray(Image);

［ca1,ch1,cv1,cd1］=dwt2(grayI,'db4');%用db4小波对图像进行一级小波分解

DWTI=［wcodemat(ca1,256),wcodemat(ch1,256);wcodemat(cv1,256),wcodemat(cd1,256)］;

%组成小波系数显示矩阵

result=idwt2(ca1,ch1,cv1,cd1,'db4');%一级重构

subplot(131),imshow(Image),title('原图');

subplot(132),imshow(DWTI/256),title('一级分解');%显示一级分解后的近似和细节图像

subplot(133),imshow(result,［］),title('一级重构');  %重构图像显示

程序运行结果如图536所示。


图536cameraman图像的一级分解及重构


【例524】对cameraman图像进行二级小波分解及重构。
程序如下： 

clc,clear,close all;

Image=imread('cameraman.jpg');


grayI=rgb2gray(Image);

［C,S］=wavedec2(grayI,2,'db4');%用db4小波对图像进行二级小波分解

siz = S(size(S,1),:); 

CA2=appcoef2(C,S,'db4',2);%提取二级小波分解低频变换系数

［CH2,CV2,CD2］ = detcoef2('all',C,S,2); %提取二级小波分解高频变换系数

［CH1,CV1,CD1］ = detcoef2('all',C,S,1);%提取一级小波分解高频变换系数

CA1=［wcodemat(CA2,256),wcodemat(CH2,256);wcodemat(CV2,256),wcodemat(CD2,256)］;

k=S(2,1)*2-S(3,1);%两级高频系数长度差

CH1=padarray(CH1,［k k］,1,'pre');

CV1=padarray(CV1,［k k］,1,'pre');

CD1=padarray(CD1,［k k］,1,'pre');%填充一级小波高频系数数组，使两级系数维数一致

DWTI=［CA1,wcodemat(CH1,256);wcodemat(CV1,256),wcodemat(CD1,256)］;

RecA=upcoef2('a',CA2,'db4',2,siz);%二级近似系数向上重构

RecV=upcoef2('v',CV2,'db4',2,siz); %垂直细节系数向上重构

result= waverec2(C,S,'db4');%二级重构

RecA=mat2gray(RecA);

RecV=mat2gray(RecV);

subplot(221),imshow(DWTI/256),title('二级分解');%显示二级分解后的近似和细节图像

subplot(222),imshow(RecA),title('二级近似系数重构'); 

subplot(223),imshow(RecV),title('垂直细节系数重构'); 

subplot(224),imshow(result,［］),title('二级重构');

程序运行结果如图537所示。


图537cameraman图像的二级分解及重构


5.5.4小波变换在图像处理中的应用

小波变换因其频率分解、多分辨分析等特性，应用于数字图像处理，可以出色地完成诸如图像滤波、图像增强、图像融合、图像压缩等多种处理，得到广泛的应用。
1. 基于小波变换的图像降噪

小波变换具有下述特点。
(1) 低熵性。图像变换后熵降低。
(2) 多分辨性。采用多分辨率的方法，可以非常好地刻画信号的非平稳特征，如边缘、尖峰、断点等，可在不同分辨率下根据信号和噪声分布的特点去噪。
(3) 小波变换可以灵活地选择不同的小波基。因此，小波去噪是小波变换在数字图像处理中的一个重要应用。

如前文所述，小波变换实际上是通过滤波器将图像信号分解为低频和高频的，噪声的大部分能量集中在高频部分，通过处理小波分解后的高频系数，实现噪声的降低。常见的基于小波变换的图像降噪方法有以下几种。
(1) 基于小波变换极大值原理的降噪方法。
该方法根据信号与噪声在小波变换各尺度上不同的传播特性，剔除由噪声产生的模极大值点，用剩余的模极大值点恢复信号。
(2) 基于相关性的降噪方法。
该方法对含噪声的信号进行变换后，计算相邻尺度间小波系数的相关性，根据相关性大小区别小波系数的类型，并进行取舍、重构。
(3) 基于阈值的降噪方法。
该方法按一定的规则(或阈值化)将小波系数划分成两类： 重要的、规则的小波系数和非重要的或受噪声干扰的小波系数，舍弃不重要的小波系数然后重构去噪后的图像。这种方法的关键是阈值的设计。常用的阈值函数有硬阈值和软阈值函数。硬阈值方法指的是设定阈值，小波系数绝对值大于阈值的保留，小于阈值的置0，可以很好地保留边缘等局部特征，但会出现振铃等失真现象； 软阈值方法指将较小的小波系数置0，较大的小波系数按一定的函数计算，向0收缩，处理结果较硬阈值方法平滑，但因绝对值较大的小波系数减小，会导致损失部分高频信息，造成图像边缘的失真模糊。
【例525】基于小波变换对图像进行硬阈值、软阈值去噪。
程序如下： 

clear,clc,close all;

Image=rgb2gray(imread('peppers.jpg'));

noiseI=imnoise(Image,'gaussian');%添加高斯噪声

［c,s］=wavedec2(noiseI,2,'sym5');%用sym5小波对图像进行二层小波分解

sigma=std(c);%小波系数标准差

thresh=2*sigma;%设定阈值

csize=size(c);

c(abs(c)<thresh)=0;%小波系数小于阈值则置0

denoiseI1=uint8(waverec2(c,s,'sym5'));

pos1=find(c>thresh);c(pos1)=c(pos1)-thresh;

pos2=find(c<-thresh);c(pos2)=c(pos2)+thresh;%大系数向0收缩

denoiseI2=uint8(waverec2(c,s,'sym5'));

subplot(221),imshow(Image),title('原图像');

subplot(222),imshow(noiseI),title('高斯噪声图像');

subplot(223),imshow(denoiseI1),title('硬阈值降噪');

subplot(224),imshow(denoiseI2),title('软阈值降噪');

程序运行结果如图538所示。


图538基于小波变换的图像降噪


2. 基于小波变换的边缘检测
图像边缘是指在图像平面中灰度值发生跳变的点连接所成的曲线段，包含了图像的重要信息。找出图像的边缘称为边缘检测，是图像处理中的重要内容。二维小波变换能检测二维函数f(x，y)的局部突变，因此是检测图像边缘的有力工具。
随着技术的发展，目前已经有很多新颖的基于小波变换的图像边缘检测技术和方法，如多尺度小波变换边缘提取算法、嵌入可信度的边缘检测方法、奇异点模极大值检测算法等。
3. 基于小波变换的图像压缩
小波变换特别适用于细节丰富、空间相关性差、冗余度低的图像数据压缩处理。同DCT类似，小波变换后使图像能量集中在少部分的小波系数上，可以通过简单的量化方法，将较小能量的小波系数省去，保留能量较大的小波系数，从而达到压缩的目的。所以，可以采用直接阈值方法实现基于小波变换的图像压缩，压缩效果好坏关键在于阈值的选择。考虑到人眼视觉系统对高频分量反应不敏感，而对低频分量反应敏感，所以，可以给低频区分配相对高的码率，给高频区分配相对低的码率，以降低数据量，如基于小波树结构的矢量量化法、嵌入式零树小波编码等。JPEG 2000压缩标准中采用基于小波变换的图像压缩技术。
4. 基于小波变换的图像增强

图像增强是指提高图像的对比度，以增加图像的视觉效果和可理解性，同时减少或抑制图像中的噪声，提高视觉质量。常用的
图像增强技术可以分为基于空间域和基于变换域两种，前者直接对像素点进行运算，后者通过将图像进行正交变换，对变换域内的系数进行调整以达到提高输出图像对比度的目的。小波变换将图像分解为大小、位置和方向不同的分量，根据需要改变某些分量系数，从而使得感兴趣的分量放大，不需要的分量减小，以达到图像增强的目的。
5. 基于小波变换的图像融合

图像融合是指将同一对象的两幅或更多的图像合成在一幅图像中，以便比原来任何一幅图像更容易为人所理解。基于小波变换的图像融合是指将原图像进行小波分解，在小波域通过一定的融合算子融合小波系数，再重构生成融合的图像，如图539所示。小波变换可以将图像分解到不同的频率域，在不同的频率域运用不同的融合算法，得到合成图像的多分辨分解，从而在合成图像中保留原图像在不同频率域的显著特征。


图539基于小波变换的图像融合过程


基于小波变换的图像融合的关键在于融合算法。例如，对于低频小波分解系数采用取平均的方法，高频分解系数的融合可以采用均值法、最大值法、基于区域的方法、基于边缘强度的方法等。

小波融合能够针对输入图像的不同特征来选择小波基及小波变换的级数，在融合时可以根据实际需要来引入双方的细节信息，表现出更强的针对性和实用性，融合效果更好。
【例526】采用DWT对图像进行融合。
程序如下： 

Image1=rgb2gray(imread('desert.jpg'));

Image2=rgb2gray(imread('car.jpg'));

［ca1,ch1,cv1,cd1］=dwt2(Image1,'db4');%用db4小波对背景图进行一级小波分解

［ca2,ch2,cv2,cd2］=dwt2(Image2,'db4');%用db4小波对前景图进行一级小波分解

ca=(ca1+ca2)/2;%DWT低频系数取平均融合

ch=max(ch1,ch2);cv=max(cv1,cv2);cd=max(cd1,cd2);  %DWT高频系数取最大值融合

result=idwt2(ca,ch,cv,cd,'db4')/256;

imshow(result),title('图像融合'); 

程序运行结果如图540所示。


图540综合实例结果图


以上是对小波变换在图像处理中的部分主要应用做了简要介绍，有兴趣的读者可以在学习过相关图像处理原理和概念后，结合小波变换的理论进行详细学习。
5.6本章小结
本章主要介绍了图像的常见正交变换，包括DFT、DCT、KL变换、Radon变换以及小波变换，详细介绍了各种正交变换的原理及MATLAB实现。正交变换在不同的处理算法中经常用到，应熟悉其基本原理、变换特点、常用函数及处理效果，以便灵活应用。