图像处理-小波变换

　什么是小波？

“小”指的是具有衰减性；“波”指的是具有波动性，其振幅正负相间的振荡形式趋于不规则、不对称，均值为0的波形，类似于下面的：

 

小波的种类：

Matlab查看小波种类：

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

>> wavemngr('read',1)   ans =     71×44 char 数组       '===================================         '     'Haar              →→haar                    '     '===================================         '     'Daubechies        →→db                      '     '------------------------------              '     'db1→db2→db3→db4→                            '     'db5→db6→db7→db8→                            '     'db9→db10→db**→                              '     '===================================         '     'Symlets           →→sym                     '     '------------------------------              '     'sym2→sym3→sym4→sym5→                        '     'sym6→sym7→sym8→sym**→                       '     '===================================         '     'Coiflets          →→coif                    '     '------------------------------              '     'coif1→coif2→coif3→coif4→                    '     'coif5→                                      '     '===================================         '     'BiorSplines       →→bior                    '     '------------------------------              '     'bior1.1→bior1.3→bior1.5→bior2.2→            '     'bior2.4→bior2.6→bior2.8→bior3.1→            '     'bior3.3→bior3.5→bior3.7→bior3.9→            '     'bior4.4→bior5.5→bior6.8→                    '     '===================================         '     'ReverseBior       →→rbio                    '     '------------------------------              '     'rbio1.1→rbio1.3→rbio1.5→rbio2.2→            '     'rbio2.4→rbio2.6→rbio2.8→rbio3.1→            '     'rbio3.3→rbio3.5→rbio3.7→rbio3.9→            '     'rbio4.4→rbio5.5→rbio6.8→                    '     '===================================         '     'Meyer             →→meyr                    '     '===================================         '     'DMeyer            →→dmey                    '     '===================================         '     'Gaussian          →→gaus                    '     '------------------------------              '     'gaus1→gaus2→gaus3→gaus4→                    '     'gaus5→gaus6→gaus7→gaus8→                    '     '===================================         '     'Mexican_hat       →→mexh                    '     '===================================         '     'Morlet            →→morl                    '     '===================================         '     'Complex Gaussian  →→cgau                    '     '------------------------------              '     'cgau1→cgau2→cgau3→cgau4→                    '     'cgau5→cgau6→cgau7→cgau8→                    '     '===================================         '     'Shannon           →→shan                    '     '------------------------------              '     'shan1-1.5→shan1-1→shan1-0.5→shan1-0.1→      '     'shan2-3→shan**→                             '     '===================================         '     'Frequency B-Spline→→fbsp                    '     '------------------------------              '     'fbsp1-1-1.5→fbsp1-1-1→fbsp1-1-0.5→fbsp2-1-1→'     'fbsp2-1-0.5→fbsp2-1-0.1→fbsp**→             '     '===================================         '     'Complex Morlet    →→cmor                    '     '------------------------------              '     'cmor1-1.5→cmor1-1→cmor1-0.5→cmor1-1→        '     'cmor1-0.5→cmor1-0.1→cmor**→                 '     '===================================         '     'Fejer-Korovkin    →→fk                      '     '------------------------------              '     'fk4→fk6→fk8→fk14→                           '     'fk18→fk22→                                  '     '===================================         '

连续小波变换#

Continuous Wavelet Transform，CWT

$\mathrm{C}(scale,position)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}f(t)\psi (scale,position,t)dt$

小波变换是将信号f（t）与被缩放和平移的小波函数$\psi ()$之积在信号存在的整个期间里求和的结果，CWT变换结果是小波系数C

缩放因子scale越小，表示小波越窄，信号频率越高；

缩放因子scale越大，表示小波越宽，信息频率越低

二维离散小波变换#

将二维信号在不同尺度上的分解，得到原始信息的近似值和细节值。

分解的结果为近似分量、水平细节分量、垂直细节分量和对角细节分量

每通过一次小波变换后，图像被分为四分之一的子频带区域，分别包含了相应频带的小波系数，LL频带是图像内容的缩略图，保持了原图的内容信息，LH频带含有水平方向的高频边缘信息，HL频带含有竖直方向的高频边缘信息，HH频带含有对角方向的高频边缘信息，反映了水平和竖直方向上图像灰度的综合变换

小波变换#

小波一次变换实际就是对图像既进行行变换又进行列变换，其结果是左上角显示缩小为原来的四分之一的图像，右上角是图像的行变换结果，左下角是图像的列变结果，而右下角则是对图像进行45度边缘检测的处理结果

百度百科：

离散小波变换是对基本小波的尺度和平移进行离散化。在图像处理时，常采用二进小波作为小波变换函数，即用2的整数次幂进行划分。

余弦变换是最经典的谱分解工具，是指整个时域过程的频域特征或整个频域过程的时域特征，故对于平稳过程，他有很好的效果，但对于非平稳过程，他却有诸多不足。在JPEG中，离散余弦变换将图像压缩为8×8 的小块，然后依次放入文件中，这种算法靠丢弃频率信息实现压缩，因而图像的压缩率越高，频率信息被丢弃的越多。在极端情况下，JPEG图像只保留了反映图像外貌的基本信息，精细的图像细节都损失了。

小波变换是现代谱分析工具，他既能考察局部时域过程的频域特征，又能考察局部频域过程的时域特征，因此即使对于非平稳过程，处理起来也得心应手。他能将图像变换为一系列小波系数，这些系数可以被高效压缩和存储，此外，小波的粗略边缘可以更好地表现图像，因为他消除了DCT压缩普遍具有的方块效应。

在数字图像处理中，需要将连续的小波及其小波变换离散化。一般计算机实现中使用二进制离散处理，将经过这种离散化的小波及其相应的小波变换成为离散小波变换（简称DWT）。实际上，离散小波变换是对连续小波变换的尺度、位移按照2的幂次进行离散化得到的，所以也称之为二进制小波变换。

虽然经典的傅里叶变换可以反映出信号的整体内涵，但表现形式往往不够直观，并且噪声会使得信号频谱复杂化。在信号处理领域一直都是使用一族带通滤波器将信号分解为不同频率分量，即将信号f(x)送到带通滤波器族Hi(x)中。

小波分解的意义就在于能够在不同尺度上对信号进行分解，而且对不同尺度的选择可以根据不同的目标来确定。

对于许多信号，低频成分相当重要，它常常蕴含着信号的特征，而高频成分则给出信号的细节或差别。人的话音如果去掉高频成分，听起来与以前可能不同，但仍能知道所说的内容；如果去掉足够的低频成分，则听到的是一些没有意义的声音。在小波分析中经常用到近似与细节。近似表示信号的高尺度，即低频信息；细节表示信号的高尺度，即高频信息。因此，原始信号通过两个相互滤波器产生两个信号。

通过不断的分解过程，将近似信号连续分解，就可以将信号分解成许多低分辨率成分。理论上分解可以无限制的进行下去，但事实上，分解可以进行到细节（高频）只包含单个样本为止。因此，在实际应用中，一般依据信号的特征或者合适的标准来选择适当的分解层数。

基本定义#

x[n]：离散的输入信号，长度为N。

g[n]：low pass filter低通滤波器，可以将输入信号的高频部份滤掉而输出低频部份。

h[n]：high pass filter高通滤波器，与低通滤波器相反，滤掉低频部份而输出高频部份。

 Q：downsampling filter降采样滤波器，如果以x[n]作为输入，则输出y[n]=x[Qn]。此处举例Q=2。

降噪#

在小波分析中经常用到近似于细节，近似表示信号的高尺度，即低频信息；细节表示信号的高尺度，即高频信息。

对含有噪声的信号，噪声分量的主要能量集中在小波解的细节分量中。

程序：

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

clc; clear all; close all;   load leleccum; % 载入信号数据   s = leleccum;   Len = length(s);   [ca1, cd1] = dwt(s, 'db1'); % 采用db1小波基分解   a1 = upcoef('a', ca1, 'db1', 1, Len); % 从系数得到近似信号   d1 = upcoef('d', cd1, 'db1', 1, Len); % 从系数得到细节信号   s1 = a1+d1; % 重构信号   figure;   subplot(2, 2, 1); plot(s); title('初始电源信号');   subplot(2, 2, 2); plot(ca1); title('一层小波分解的低频信息');   subplot(2, 2, 3); plot(cd1); title('一层小波分解的高频信息');   subplot(2, 2, 4); plot(s1, 'r-'); title('一层小波分解的重构信号');

总的来说，

小波分解可以使人们在任意尺度观察信号，只需所采用的小波函数的尺度合适。小波分解将信号分解为近似分量【低频】和细节分量【高频】，它们在应用中分别有不同的特点。比如，对含有噪声的信号，噪声分量的主要能量集中在小波分解的细节分量中，对细节分量做进一步处理，比如阈值处理，可以过滤噪声。

Matlab实现#

常用函数#

实验1#

连续小波分析只需要一个cwt函数

下面的例子是一个包含噪声的正弦波：

1、加载信号

?

1	load noissin

可以使用whos显示信号信息

?

1 2 3 4 5 6 7 8 9 10 11 12

>> whos   Name          Size              Bytes  Class     Attributes     Len           1x1                   8  double                a1            1x4320            34560  double                ca1           1x2160            17280  double                cd1           1x2160            17280  double                d1            1x4320            34560  double                leleccum      1x4320            34560  double                noissin       1x1000             8000  double                s             1x4320            34560  double                s1            1x4320            34560  double

2、执行连续小波变换

?

1	c=cwt(noissin,1:48,'db4');

函数cwt 的参数分别为分析的信号，分析的尺度和使用的小波

返回值c包含了在各个尺度下的小波系数，这里c是一个48*1000的矩阵，每一行与一个尺度相关

3、绘制小波系数

cwt函数可以接受第四个参数，来指定函数在执行结束后是否会绘制连续小波变换系数的绝对值

下面语句就是绘制系数结果：

?

1	b=cwt(noissin,1:48,'db4','plot');

4、选择分析的尺度

cwt函数的第二个参数可以设定任意小波分析的尺度，只要满足以下要求即可：

（1）所有尺度必须是正实数

（2）尺度的增量必须为正

（3）最高尺度不能超过由信号决定的一个最大值

如下面代码可以执行从2开始的偶数尺度计算：

?

1	b=cwt(noissin,2:2:128,'db4','plot');

上面图像显示信号的周期性

实验2#

1、加载信号

?

1	load leleccum

截取信号：

?

1 2	s=leleccum(1:3920); ls=length(s);

2、对信号进行一层小波分解

使用db1小波执行一层小波分解，执行下面语句产生近似低频CA1，高频系数CD1

?

1	[CA1,CD1]=dwt(s,'db1');

3、从系数中构建低频和高频

从系数CA1和CD1中构建一层低频A1和高频D1：

?

1 2 3 4 5 6 7 8 9

A1 = upcoef('a',cA1,'db1',1,ls);   D1 = upcoef('d',cD1,'db1',1,ls);   或   A1 = idwt(cA1,[],'db1',ls);   D1 = idwt([],cD1,'db1',ls);

系数重构：

4、显示低频和高频

?

1 2	subplot(1,2,1); plot(A1); title('Approximation A1') subplot(1,2,2); plot(D1); title('Detail D1')

5、使用逆小波变换恢复信号

?

1 2	A0=idwt(CA1,CD1,'db1',ls); subplot(2,2,3); plot(A0); title('逆变换 ')

6、多尺度一维分解

执行3层信号分解：

?

1	[C,L] = wavedec(s,3,'db1');

函数返回3层分解的各组分系数C（连接在一个向量里），向量L里返回的是各组分的长度

7、抽取低频和高频系数

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

从C中抽取3层低频系数   CA3 = appcoef(C,L,'db1',3);   从C中抽取3、2、1层高频系数   CD3 = detcoef(C,L,3);   CD2 = detcoef(C,L,2);   CD1 = detcoef(C,L,1);   或者   [CD1,CD2,CD3] = detcoef(C,L,[1,2,3]);   figure; subplot(2,2,1); plot(CA3); title('低频系数') subplot(2,2,2); plot(CD3); title('高频系数三层') subplot(2,2,3); plot(CD2); title('高频系数二层') subplot(2,2,4); plot(CD1); title('高频系数一层')

8、重构三层低频和1，2，3层高频

?

1 2 3 4 5 6 7 8 9 10 11

从C中重建3层近似   A3 = wrcoef('a',C,L,'db1',3);   从C中重建1、2、3层细节   D1 = wrcoef('d',C,L,'db1',1);   D2 = wrcoef('d',C,L,'db1',2);   D3 = wrcoef('d',C,L,'db1',3);

9、显示分层分解的结果

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

显示3层分解的结果   subplot(2,2,1); plot(A3);   title('Approximation A3')   subplot(2,2,2); plot(D1);   title('Detail D1')   subplot(2,2,3); plot(D2);   title('Detail D2')   subplot(2,2,4); plot(D3);   title('Detail D3')

10、从3层分解中重建原始信号

?

1 2	A0 = waverec(C,L,'db1'); err = max(abs(s-A0))

11、粗糙的去燥信号

使用小波从信号中移除噪声需要辨识哪个或哪些组分包含噪音，然后重建没有这些组分的信号

下面的例子中，连续的低频随着越来越多的高频信息从信号中滤除，噪声变得越来越少，3层低频与原始信号对比会更干净：

?

1 2 3

figure; subplot(2,1,1);plot(s);title('Original'); axis off subplot(2,1,2);plot(A3);title('Level 3 Approximation');axis off

 

当然，除去所有高频信息，会失去原始信号中的很多最尖锐的特征，最佳去噪音的方法是通过一种更精细的“阈值”方法，只丢弃部分超过一定范围的细节

12、通过阈值去除噪音

 先看3层分析的高频：

?

1 2 3 4

figure; subplot(3,1,1); plot(D1); title('Detail Level 1'); axis off subplot(3,1,2); plot(D2); title('Detail Level 2'); axis off subplot(3,1,3); plot(D3); title('Detail Level 3'); axis off

从图中可以看到，大多数噪声发生在信号的后面部分，表现在细节上就是出现大波动的地方。如果我们通过设定最大值来限定细节强度，会怎么样呢？这会有降低噪声效果，同时保留不影响的必要细节。这是一种很好的方法。

注意到CD1,CD2,CD3是向量，那么我们就可以通过直接操纵这些向量来达到目的，即设置这些向量小于其峰值或平均值的一部分，然后就可以由这些设定了阈值的系数重建新的细节信号D1、D2、D3。

实际去噪过程中，可以使用ddencmp函数来计算默认的阈值参数，然后用wdencmp函数来执行实际的去噪过程，代码如下

?

1 2	[thr,sorh,keepapp] = ddencmp('den','wv',s); clean = wdencmp('gbl',C,L,'db1',3,thr,sorh,keepapp);

注意wdencmp使用了第6步中小波分解的结果C、L，另外指定了db1小波来做分析，指定全局阈值选项gb1.

详细参见ddencmp函数和wdencmp函数。

显示原始信号和去噪信号如下：

?

1 2 3

figure; subplot(2,1,1); plot(s(2000:3920)); title('Original') subplot(2,1,2); plot(clean(2000:3920)); title('denoised')

这里我们只绘制了原始信号中包含噪声的部分，特别注意我们是如何在移除了噪声的情况下仍保持原有的尖锐细节的，这也是小波分析强大的地方。

参考#

1、知乎1   知乎2

2、博客园1

3、CSDN1

4、Matlab实现1   Matlab实现2

5、数字图像处理（杨淑莹）