Matlab系列之小波分析基础

2020-12-22 15:15:24 浏览数 (1)

前言

原本想把MATLAB里关于概率论的相关进行记录，不过概率论学得不好，感觉在该部分的表达上还存在很大不足，就放弃了相关的篇章，直接开始了本篇，本篇主要是记录小波分析的一些东西，小波分析的原理就不细说了，所以还是老样子，主要介绍小波分析在MATLAB中的相关知识，不足之处请指出。

介绍

小波分析是数学分析方法里的一种，主要应用于信号处理、图像处理、语音分析以及其他的非线性科学领域，它被认为是继Fourier分析之后的又一有效的时频分析方法。小波变换与Fourier变换相比，是一个时间和频域的局域变换因而能有效地从信号中提取信息，通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析（Multiscale Analysis），解决了Fourier变换不能解决的许多困难问题。

MATLAB提供了小波分析工具箱，在主界面的命令窗口输入：wavemenu，就可以打开工具箱，如下所示。

常用的就是小波基函数、连续小波变换及其应用、离散小波变换及其应用、小波包变换、信号和图像的多尺度分解、基于小波变换的信号去噪、信号压缩，在上图也可以找到与这些对应的选项。常用的小波基函数如下表：

函数表示

函数描述

morl

Morlet小波

mexh

墨西哥草帽小波

meyr

Meyer小波

haar

Haar小波

dbN

紧支集正交小波

symN

近似对称的紧支集正交小波

coifN

Coifmant小波

biorNr.Nd

双正交样条小波

以下记录的是一些常用指令和语法使用，工具箱的操作就不弄了，自行根据指令进行对应和补充即可。

1、waveinfo函数

note：information on wavelets.

该语法的功能是提供工具箱中所有小波的信息查询，使用格式：waveinfo('wname')

wname指代的小波有

'haar' : Haar wavelet. 'db' : Daubechies wavelets. 'sym' : Symlets. 'coif' : Coiflets. 'bior' : Biorthogonal wavelets. 'rbio' : Reverse biorthogonal wavelets. 'meyr' : Meyer wavelet. 'dmey' : Discrete Meyer wavelet. 'gaus' : Gaussian wavelets. 'mexh' : Mexican hat wavelet. 'morl' : Morlet wavelet. 'cgau' : Complex Gaussian wavelets. 'cmor' : Complex Morlet wavelets. 'shan' : Complex Shannon wavelets. 'fbsp' : Complex Frequency B-spline wavelets. 'fk' : Fejer-Korovkin orthogonal wavelets

使用举例：waveinfo('haar')

结果：

查询小波包的信息，则使用：waveinfo('wsys')

2、wfilters函数

note：Wavelet filters.

这个函数有两种用法，结果也不太相同；

第一种：[LOD,HID,LOR,HIR] = wfilters('wname')

计算正交小波或双正交小波（wname）有关联的四个滤波器，分别为：

LO_D，分解低通滤波器 HI_D，分解高通滤波器 LO_R，重构低通滤波器 HI_R，重构高通滤波器

第二种：[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（高通滤波器）

应用实例

代码语言：javascript复制
clc;
clear all;
[LO_D,HI_D,LO_R,HI_R]=wfilters('db45');%多贝西小波
%stem产生离散序列的图形，xlim限制x轴的长度，另外三个就是对图形做描述
subplot(221);stem(LO_D);xlim([0 95]);title('分解低通滤波器');
xlabel('x');ylabel('y');
subplot(222);stem(HI_D);xlim([0 95]);title('分解高通滤波器');
xlabel('x');ylabel('y');
subplot(223);stem(LO_R);xlim([0 95]);title('重构低通滤波器');
xlabel('x');ylabel('y');
subplot(224);stem(HI_R);xlim([0 95]);title('重构高通滤波器');
xlabel('x');ylabel('y');

运行结果

3、dwt函数

note：Single-level discrete 1-D wavelet transform.

该函数是通过指定小波‘wname’或者指定小波的分解滤波器LOD和HID执行单层一维小波分解。

使用方法：

[CA,CD] = dwt(X,'wname')或[CA,CD] = dwt(X,LO_D,HI_D)

还有一种是指定模式计算小波分解

[cA,cD] = dwt(...,'mode',MODE)

该用法的举例：[cA,cD] = dwt(x,'db1','mode','sym')；

mode对应值及其英语表述如下：

'mode'

DWT Extension Mode

'sym' or 'symh'

Symmetric-padding (half-point): boundary value symmetric replication — default mode

'symw'

Symmetric-padding (whole-point): boundary value symmetric replication

'asym' or 'asymh'

Antisymmetric-padding (half-point): boundary value antisymmetric replication

'asymw'

Antisymmetric-padding (whole-point): boundary value antisymmetric replication

'zpd'

Zero-padding

'spd' or 'sp1'

Smooth-padding of order 1 (first derivative interpolation at the edges)

'sp0'

Smooth-padding of order 0 (constant extension at the edges)

'ppd'

Periodic-padding (periodic extension at the edges)

二维的类似，函数为dwt2，使用语法为：

代码语言：javascript复制
[CA,CH,CV,CD] = dwt2(X,'wname')
 [CA,CH,CV,CD] = dwt2(X,Lo_D,Hi_D) 
 [CA,CH,CV,CD] = dwt2(...,'mode',MODE)

小波分解的逆过程就是小波重构，类似FFT和IFFT，很多时候傅里叶变换也被人拿来和小波变换作一些比对。

iwdt和iwdt2的使用语法分别如下所示：

代码语言：javascript复制
%iwdt使用语法
X = idwt(CA,CD,'wname',L) 
X = idwt(CA,CD,LO_R,HI_R,L)
X = idwt(...,'mode',MODE)
%iwdt2使用语法
X = idwt2(CA,CH,CV,CD,'wname') 
X = idwt2(CA,CH,CV,CD,Lo_R,Hi_R) 
X = = idwt2(...,'mode',MODE)

除了以上的语法以外，在这些基础上还有一些扩展的语法，这些就留给“help”了。

应用实例

代码语言：javascript复制
clc;
clear all;
close all;
n=randn(1,256);
z=1.5*sin(1:256);
s=n z;
[ca1,cd1]=dwt(s,'haar');%%等同于db1
subplot(311);plot(s,'b-');title('原始信号');
xlabel('x');ylabel('y');
subplot(323);plot(ca1,'b-');title('haar低频系数1');
xlabel('x');ylabel('y');
subplot(324);plot(cd1,'b-');title('haar高频系数1');
xlabel('x');ylabel('y');
%计算滤波器的相关值，再直接指定分解滤波器进行系数的计算
[LO_D,HI_D]=wfilters('haar','d');
[ca2,cd2]=dwt(s,LO_D,HI_D);
subplot(325);plot(ca2,'b-');title('haar低频系数2');
xlabel('x');ylabel('y');
subplot(326);plot(cd2,'b-');title('haar高频系数2');
xlabel('x');ylabel('y');

运行结果

4、wavedec

note：Multi-level 1-D wavelet decomposition.

该函数的功能依然是小波分解，只是其层级变多了，所以可以猜到其语法和dwt会有点相似，如下：

代码语言：javascript复制
[C,L] = wavedec(X,N,'wname') 
[C,L] = wavedec(X,N,LO_D,HI_D),

N必须是正整数，其用意就是，返回信号X在N层的小波分解，C代表的是分解向量，L代表一个记录向量。

应用实例

代码语言：javascript复制
%对音频进行小波分解
clc;
clear all;
close all;
%产生声音文件
load handel.mat
filename = 'dzkr_SoundTest.wav';
audiowrite(filename,y,Fs);
clear y Fs
%将生成的音频导入
s=audioread('dzkr_SoundTest.wav');%旧版本可能用的是wavread()，需要自行比对下版本差别
n=length(s);
figure;
subplot(211);plot(4000:n,s(4000:n));
xlim([4000 6200]);ylim([-1 1]);
xlabel('信号序列');ylabel('信号值');
title('原始声音图像');
[C,L]=wavedec(s(4000:n),2,'db2');
subplot(212);plot(C);xlim([0 2200]);ylim([-2 2]);
title('小波分解结构');
xlabel('低频系数 and 第二层及第一层高频的信号序列');ylabel('信号值');

运行结果

wavedec也有二维的形式，即wavedec2，语法如下：

代码语言：javascript复制
[C,S] = wavedec2(X,N,'wname') 
[C,S] = wavedec2(X,N,LO_D,HI_D)

其中S则变成了记录矩阵，C依然还是分解向量，其他的参数与以上的一致。

很理所当然的就可以想到逆过程了，不过可不再是加个i了，其逆过程变成了waverec，二维则是waverec2，两者的语法分别为：

代码语言：javascript复制
%waverec
X = waverec(C,L,'wname') 
X = waverec(C,L,LO_R,HI_R)
%waverec2
X = waverec2(C,S,'wname') 
X = waverec2(C,S,LO_R,HI_R)

5、wrcoef

这个函数也是重构的功能，和waverec有点类似，用法如下：

代码语言：javascript复制
X = wrcoef('type',C,L,'wname',N) 
X = wrcoef('type',C,L,LO_R,HI_R,N) 
X = wrcoef('type',C,L,'wname') 
X = wrcoef('type',C,L,LO_R,HI_R)

参数上和waverec的也基本一致，不过多了一个‘type’，而且N

有N的就是相当于进行N层重构，且N的值和‘type’有关系，type可取低频‘a’或高频'd'，取低频‘a’的时候，N最小可以为0；两种类型下N的值都得满足于：N <= length(L)-2，若没有指定N，则默认在N = length(L)-2进行重构。

应用实例

代码语言：javascript复制
clc;
clear all;
close all;
N=1000;
t=1:N;
%产生信号
sig1=sin(0.3*t);%正弦波
%三角波
sig2(1:500)=((1:N/2)-1)/500;
sig2(501:N)=(N-((N/2 1):1000))/500;
sig=sig1 sig2;%叠加信号
[C,L]=wavedec(sig,2,'db6');%进行小波分解
a1=wrcoef('a',C,L,'db6',1);%重构第1层逼近系数
a2=wrcoef('a',C,L,'db6',2);%重构第2层逼近系数
d1=wrcoef('d',C,L,'db6',1);%重构第1层细节系数
d2=wrcoef('d',C,L,'db6',2);%重构第2层细节系数
subplot(511);plot(sig);xlabel('信号序列');ylabel('原始信号');
subplot(512);plot(a1);xlabel('信号序列');ylabel('a1');
subplot(513);plot(a2);xlabel('信号序列');ylabel('a2');
subplot(514);plot(d1);xlabel('信号序列');ylabel('d1');
subplot(515);plot(d2);xlabel('信号序列');ylabel('d2');

运行结果

结语

本篇暂告一段落，仔细看完的话，你会发现本篇介绍到的小波分析展示了其”选取滤波器“的功能，之后还会写一篇用小波分析的知识做一些图像处理，比如图像去噪和图像压缩，音频的话，本篇已经略微涉及到了音频信号简单分解，深入了解就不折腾了，还是等该系列的下篇关于小波分析在图像处理的部分应用吧~

matlab padding point replication

0 人点赞