第一篇:matlab小波变换函数的总结与程序
小波去噪举例
MATLAB中用wnoise函数测试去噪算法
sqrt_snr=3;init=231434;
[x,xn]=wnoise(3,11,sqrt_snr,init);% 加噪,信噪比为3 subplot(3,2,1),plot(x)
title('original test function')subplot(3,2,2),plot(xn)title('noised function')lev=5;
xd=wden(x,'heursure','s','one',lev,'sym8');%利用小波对一维信号进行降噪, XD为降噪后的%信号,CXD,LXD为XD的小波分解结构 % 's' or 'h'决定阈值的使用方式,SCAL决定阈值是%否随噪声变化:'one' 不调整,'sln'对第一层系%数的层噪声分别进行估计和调整; 'mln'对各层%系数的层噪声分别进行估计和调整;
subplot(3,2,3),plot(xd)
title('One de-noised function')xd=wden(x,'heursure','s','sln',lev,'sym8');subplot(3,2,4),plot(xd)
title('Sln de-noised function')xd=wden(x,'sqtwolog','s','sln',lev,'sym8');% 固定阈值选择算法去噪 subplot(3,2,5),plot(xd)
title('Sqtwolog de-noised function')
[c,l]=wavedec(x,lev,'sym8');subplot(3,2,6),plot(xd)
title('CL de-noised function')
MATLAB中图像噪声处理举例
load sinsin;
colormap('default');subplot(1,3,1),image(X);title('original image');axis('square');init=231434;
randn('seed',init);
X=X+18*randn(size(X));%产生噪声信号 subplot(1,3,2),image(x);title('noised image');axis('square');[thr,sorh,keepapp]=ddencmp('den','wv',x);%自动生成小波去躁或压缩的阈值选择方案,也 %就是寻找默认值
[xc,cxc,lxc,perf0,perfl2]=wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);%使用全局阈值进行
%图象降噪
subplot(1,3,3),image(xc);title('denoised image');axis('square')
可见,含躁图像的噪声含量很强,利用小波去躁,可以有效去除躁声,同时保留了边界。Wdencmp函数
[xc,cxc,lxc,perf0,perfl2]=wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp)是使用小波进行一维或二维小波压缩或降噪的函数。前面的语句是对于输入的一维或二维信号X,使用全局正阈值THR,由小波系数阈值得到降噪或压缩后的信号XC。附加的输出变量[cxc,lxc]是XC的小波分解结构; [perf0,perfl2]是恢复和压缩的L2范数百分比。使用小波'sym4'执行小波分解到第N=2层。Sorh是软阈值或硬阈值。若keepapp=1,低频系数不能进行阈值处理。
第二篇:小波变换 matlab 总结范文
小波变换matlab总结 目录
一、预置工具...........................................................................................................................4
1.预置信号......................................................................................................................4 2.预置小波......................................................................................................................4 3.滤波器函数..................................................................................................................6
wfilters函数.......................................................................................................6 4.量化编码......................................................................................................................6
wcodemat函数.......................................................................................................6 5.阈值获取......................................................................................................................6
ddencmp函数.........................................................................................................6 thselect函数.......................................................................................................7 wbmpen函数...........................................................................................................7 wdcbm函数..............................................................................................................7 6.阈值去噪......................................................................................................................8
wden函数................................................................................................................8 wdencmp函数.........................................................................................................8 wthresh函数.........................................................................................................9 wthcoef函数.........................................................................................................9 wpdencmp函数.......................................................................................................9
二、小波变换函数.................................................................................................................12
单尺度一维小波变换.....................................................................................................12
cwt一维连续小波变换.........................................................................................12 dwt一维离散小波变换.........................................................................................12 idwt一维离散小波逆变换..................................................................................13 upcoef 一维小波系数重构................................................................................13 多尺度一维小波变换.....................................................................................................14
wavedec多尺度一维分解...................................................................................14 waverec多尺度一维重构...................................................................................15 appcoef低频系数提取.......................................................................................16 detcoef高频系数提取.......................................................................................16 wrcoef多尺度小波系数重构.............................................................................17 一维静态(平稳)小波变换.........................................................................................18
swt一维平稳小波变换.........................................................................................18 iswt一维平稳小波逆变换..................................................................................18 实例.........................................................................................................................19 单尺度二维小波变换.....................................................................................................19
dwt2二维离散小波变换......................................................................................19 idwt2二维离散小波逆变换................................................................................20 upcoef2二维系数重构.......................................................................................20 多尺度二维小波变换.....................................................................................................21
wavedec2多尺度二维分解.................................................................................21 waverec2多尺度二维重构.................................................................................22 appcoef2低频系数提取.....................................................................................23 detcoef2高频系数提取.....................................................................................23 wrcoef2多尺度小波系数重构...........................................................................24 二维静态(平稳)小波变换.........................................................................................26
swt2二维静态小波变换......................................................................................26 iswt2二维静态小波逆变换................................................................................26 实例.........................................................................................................................26 直接调用的小波函数.....................................................................................................28
meyer函数............................................................................................................28 cgauwavf函数.....................................................................................................28 mexihat函数.......................................................................................................28 morlet函数.........................................................................................................29 symwavf函数.......................................................................................................29
三、图像接口调用.................................................................................................................30
使用图形接口做一维连续小波分析.....................................................................30 使用图形接口做一维离散小波分析.....................................................................33 使用图形接口分析复信号.....................................................................................36 使用图形接口做一维除噪分析.............................................................................36
四、小波变换在图像处理中的应用.....................................................................................40
4.1 小波分析用于图像压缩........................................................................................40
4.1.1 基于小波变换的图像局部压缩...............................................................40 4.1.2 小波变换用于图像压缩的一般方法.......................................................41 4.1.3 基于小波包变换的图像压缩...................................................................45 4.2 小波分析用于图像去噪........................................................................................47
小噪声阈值去噪.....................................................................................................48 大噪声滤波去噪.....................................................................................................49 少量噪声的小波分解系数阈值量化去噪.............................................................50 4.3 小波分析用于图像增强........................................................................................52
4.3.1 图像增强问题描述...................................................................................52 4.3.2 图像钝化...................................................................................................53 4.3.3 图像锐化...................................................................................................54 4.4 小波分析用于图像融合........................................................................................56 4.5 小波分析用于图像分解........................................................................................57
一、预置工具 1.预置信号
Matlab 内置了大量的信号实例,进行信号试验的时候可以调用。调用格式:load signal matlab安装目录的toolbox/wavelet/wavedemo noissin 包含噪声的正弦波
leleccum一维电压信号,有365560个采样点 wbarb 专指图片:
2.预置小波
Matlab预置了共计15种小波。查看语句:wavemngr('read',1)查看单个小波:waveinfo('wname')
1.Haar小波 小波名haar
2.Daubechies小波系 小波名db 调用名db1db2db3db4 db5 db6 db7 db8 db9 db10
3.对称小波系Symlets 小波名sym 调用名sym2sym3 sym4 sym5 sym6 sym7 sym8
4.Coiflets 小波系 小波名coif 调用名 coif1coif2coif3coif4coif5 5.Biorthogonal小波系 小波名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
6.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
7.Meyer小波 小波名meyr
8.Dmeyer小波 小波名dmey
9.Gaussian小波系 小波名gaus 调用名gaus1gaus2gaus3gaus4gaus5gaus6gaus7gaus8
10.Mexican hat小波 小波名mexh
11.Morlet小波 小波名morl
12.Complex Gaussian小波系 小波名cgau 调用名cgau1cgau2cgau3cgau4cgau5 cgau
13.Shannon小波系 小波名shan 调用名shan1-1.5shan1-1shan1-0.5shan1-0.1shan2-3
14.Frequency B小波系 小波名fbsp 调用名fbsp1-1-1.5fbsp1-1-1fbsp1-1-0.5fbsp2-1-1fbsp2-1-0.5fbsp2-1-0.1
15.Complex Morlet小波系 小波名cmor 调用名cmor1-1.5cmor1-1cmor1-0.5cmor1-1cmor1-0.5cmor1-0.1 3.滤波器函数 wfilters函数
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('wname')计算对应小波的滤波器。
The four output filters are Lo_D, 用于分解的低通滤波器 Hi_D, 用于分解的高通滤波器 Lo_R, 用于重构的低通滤波器 Hi_R, 用于重构的高通滤波器
4.量化编码 wcodemat函数
y=wcodemat(x,nb,opt,absol)y=wcodemat(x,nb,opt)y=wcodemat(x,nb)y=wcodemat(x)该函数是用来对矩阵X进行量化编码,它返回矩阵X的一个编码矩阵,在编码中,把矩阵X中元素绝对值最大的作为NB(NB是一个整数),绝对值最小的作为1,其他元素依其绝对值的大小在1与NB中排列.当OPT为'row'时,做行编码;当OPT为'col'时,做列编码, 当OPT为'mat'时,做全局编码,即把整个矩阵中的元素绝对值最大的元素作为NB,最小的作为1,其他元素依其绝对值的大小在整个矩阵中排列.当ABSOL为0时,该函数返回输入矩阵X的一个编码版本;为非0时,返回X的绝对值
5.阈值获取
MATLAB中实现阈值获取的函数有ddencmp、thselect、wbmpen和wwdcbm,下面对它们的用法进行简单的说明。ddencmp函数
调用格式:
[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,IN2,X)函数ddencmp用于获取信号在消噪或压缩过程中的默认阈值。输入参数X为一维或二维信号;
IN1取值为'den'或'cmp','den'表示进行去噪,'cmp'表示进行压缩; IN2取值为'wv'或'wp',wv表示选择小波,wp表示选择小波包。输出参数THR为函数选择的阈值,SORH为函数选择阈值使用方式。Sorh=s,为软阈值;Sorh=h,为硬阈值;输出参数KEEPAPP决定了是否对近似分量进行阈值处理。可选为0或1。CRIT为使用小波包进行分解时所选取的熵函数类型。例:自动生成信号小波处理的阈值选取方案。r=2055415866;randn('seed',r);x=randn(1,1000);%产生白噪声;%求取对信号进行小波消噪处理的默认阈值、软阈值,并且保留低频系数;[thr,sorh,keepapp]=ddencmp('den','wv',x);输出:
thr = 3.8593 sorh = s keepapp = 1 thselect函数
调用格式如下:
THR=thselect(X,TPTR);根据字符串TPTR定义的阈值选择规则来选择信号X的自适应阈值。自适应阈值的选择规则包括以下四种:
*TPTR='rigrsure',自适应阈值选择使用Stein的无偏风险估计原理。*TPTR='heursure',使用启发式阈值选择。
*TPTR='sqtwolog',阈值等于sqrt(2*log(length(X))).*TPTR='minimaxi',用极大极小原理选择阈值。
阈值选择规则基于模型 y = f(t)+ e,e是高斯白噪声N(0,1)。wbmpen函数
调用格式:
THR=wbmpen(C,L,SIGMA,ALPHA);返回去噪的全局阈值THR。THR通过给定的一种小波系数选择规则计算得到,小波系数选择规则使用Birge-Massart的处罚算法。
[C,L]是进行去噪的信号或图像的小波分解结构; SIGMA是零均值的高斯白噪声的标准偏差;
ALPHA是用于处罚的调整参数,它必须是一个大于1的实数,一般去ALPHA=2。wdcbm函数 调用格式:
(1)[THR,NKEEP]=wdcbm(C,L,ALPHA);(2)[THR,NKEEP]=wdcbm(C,L,ALPHA,M);函数wdcbm是使用Birge-Massart算法获取一维小波变换的阈值。返回值THR是与尺度无关的阈值,NKEEP是系数的个数。
[C,L]是要进行压缩或消噪的信号在j=length(L)-2层的分解结构; LAPHA和M必须是大于1的实数;
THR是关于j的向量,THR(i)是第i层的阈值;
NKEEP也是关于j的向量,NKEEP(i)是第i层的系数个数。一般压缩时ALPHA取1.5,去噪时ALPHA取3.6.阈值去噪
MATLAB中实现信号的阈值去噪的函数有wden、wdencmp、wthresh、wthcoef、wpthcoef以及wpdencmp。下面对它们的用法作简单的介绍。wden函数 调用格式:
(1)[XD,CXD,LXD]=wden(X,TPTR,SORH,SCAL,N,'wname')(2)[XD,CXD,LXD]=wden(C,L,TPTR,SORH,SCAL,N,'wname')函数wden用于一维信号的自动消噪。X为原始信号,[C,L]为信号的小波分解,N为小波分解的层数。
THR为阈值选择规则:
*TPTR='rigrsure',自适应阈值选择使用Stein的无偏风险估计原理。*TPTR='heursure',使用启发式阈值选择。
*TPTR='sqtwolog',阈值等于sqrt(2*log(length(X))).*TPTR='minimaxi',用极大极小原理选择阈值。SORH是软阈值或硬阈值的选择(分别对应's'和'h')。SCAL指所使用的阈值是否需要重新调整,包含下面三种: *SCAL='one' 不调整;
*SCAL='sln' 根据第一层的系数进行噪声层的估计来调整阈值。*SCAL='mln' 根据不同的噪声估计来调整阈值。
XD为消噪后的信号,[CXD,LXD]为消噪后信号的小波分解结构。格式(1)返回对信号X经过N层分解后的小波系数进行阈值处理后的消噪信号XD和信号XD的小波分解结构[CXD,LXD]。格式(2)返回参数与格式(1)相同,但其结构是由直接对信号的小波分解结构[C,L]进行阈值处理得到的。
wdencmp函数
调用格式有以下三种:
(1)[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('gbl',X,'wname',N,THTR,SORH,KEEPAPP);(2)[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',X,'wname',N,THTR,SORH);(3)[XC,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',C,L,'wname',N,THTR,SORH);函数wdencmp用于一维或二维信号的消噪或压缩。wname是所用的小波函数。
gbl(global的缩写)表示每一层都采用同一个阈值进行处理。lvd表示每层采用不同的阈值进行处理。N表示小波分解的层数。THR为阈值向量,对于格式(2)和(3)每层都要求有一个阈值,因此阈值向量THR的长度为N。
SORH表示选择软阈值或硬阈值(分别取值为's'和'h'),参数KEEPAPP取值为1时,则低频系数不进行阈值量化,反之,低频系数要进行阈值量化。XC是要进行消噪或压缩的信号,[CXC,LXC]是XC的小波分解结构,PERF0和PERFL2是恢复或压缩L^2的范数百分比。如果[C,L]是X的小波分解结构,则PERFL2=100*(CXC向量的范数/C向量的范数)^2;如果X是一维信号,小波wname是一个正交小波,则PERFL2=100||XC||^2/||X||^2。
wthresh函数
调用格式:
Y=wthresh(X,SORH,T)返回输入向量或矩阵X经过软阈值(如果SORH='s')或硬阈值(如果SORH='h')处理后的信号。T是阈值。
Y=wthresh(X,'s',T)返回的是Y=SIG(X)*(|X|-T)+,即把信号的绝对值与阈值进行比较,小于或等于阈值的点变为零,大于阈值的点为该点值与阈值的差值。
Y=wthresh(X,'h',T)返回的是Y=X*1(|X|>T),即把信号的绝对值和阈值进行比较,小于或等于阈值的点变为零,大于阈值的点保持不变。一般来说,用硬阈值处理后的信号比用软阈值处理后的信号更粗糙。wthcoef函数 调用格式:
(1)NC=wthcoef('d',C,L,N,P)(2)NC=wthcoef('d',C,L,N)(3)NC=wthcoef('a',C,L)(4)NC=wthcoef('t',C,L,N,T,SORH)用于一维信号小波系数的阈值处理。
格式(1)返回小波分解结构[C,L]经向量N和P定义的压缩率处理后的新的小波分解向量NC,[NC,L]构成一个新的小波分解结构。N包含被压缩的细节向量,P是把较小系数置0的百分比信息的向量。N和P的长度必须相同,向量N必须满足1<=N(i)<=length(L)-2。格式(2)返回小波分解结构[C,L]经过向量N中指定的细节系数置0后的小波分解向量NC。格式(3)返回小波分解结构[C,L]经过近似系数置0后的小波分解向量NC。
格式(4)返回小波分解结构[C,L]经过将向量N作阈值处理后的小波分解向量NC。如果SORH=‟s„,则为软阈值;如果SORH='h'则为硬阈值。N包含细节的尺度向量,T是N相对应的阈值向量。N和T的长度必须相等。
wpdencmp函数
调用格式:
[XD,TREED,PERF0,PERFL2]=wpdencmp(X,SORH,N,'wname',CRIT,PAR,KEEPAPP)[XD,TREED,PERF0,PERFL2]=wpdencmp(TREE,SORH,CRIT,PAR,KEEPAPP)函数wpdencmp用于使用小波包变换进行信号的压缩或去噪。格式(1)返回输入信号X(一维或二维)的去噪或压缩后的信号XD。输出参数TREED是XD的最佳小波包分解树;
PERFL2和PERF0是恢复和压缩L2的能量百分比。PERFL2=100*(X的小波包系数范数/X的小波包系数)^2;如果X是一维信号,小波wname是一个正交小波,则PERFL2=100*||XD||^2/||X||^2。
SORH的取值为's'或'h',表示的是软阈值或硬阈值。
输入参数N是小波包的分解层数,wname是包含小波名的字符串。
函数使用由字符串CRIT定义的熵和阈值参数PAR实现最佳分解。如果KEEPAPP=1,则近似信号的小波系数不进行阈值量化;否则,进行阈值量化。
格式(2)与格式(1)的输出参数相同,输入选项也相同,只是它从信号的小波包分解树TREE进行去噪或压缩。
实例:
load noisdopp;x=noisdopp;%给出全局阈值
[thr,sorh,keepapp]=ddencmp('den','wv',x);%根据全局阈值对信号降噪
xc=wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);[c,l]=wavedec(x,2,'sym4');[thr1,nkeep]=wdcbm(c,l,3);xc1=wdencmp('lvd',c,l,'sym4',2,thr1,'s');subplot(311);plot(x);title('原始信号')subplot(312);plot(xc);title('使用全局阈值降噪后的信号');subplot(313);plot(xc1);title('使用分层阈值降噪后的信号');
load woman %下面进行噪声的产生 init=3;rand('seed',init);Xn=X+18*(rand(size(X)));%下面进行图像的去噪处理
[thr,sorh,keepapp]=ddencmp('den','wv',Xn);Xd=wdencmp('gbl',Xn,'sym5',2,thr,sorh,keepapp);%用sym5小波对图像信号进行二层的小波分解 [c,s]=wavedec2(X,2,'sym5');[thr1,nkeep1]=wdcbm2(c,s,1.5);Xd1=wdencmp('lvd',c,s,'sym5',2,thr1,'s');%显示图像
subplot(221);image(X);colormap(map);title('原始图像X');axis square subplot(222);image(Xn);colormap(map);title('含噪声的图像');axis square subplot(223);image(Xd);colormap(map);title('去噪后的图像');axis square subplot(224);image(Xd1);colormap(map);title('分层去噪后的图像');axis square
二、小波变换函数 单尺度一维小波变换 cwt一维连续小波变换 调用格式:
coefs = cwt(s,scale,’wname’)【仅给出系数值】
coefs = cwt(s,scale,’wname’,’plot’)【给出系数值同时生成图像】 实例:
load noissin;figure subplot(121)plot(noissin);subplot(122)c=cwt(noissin,1:48,'db4','plot');
dwt一维离散小波变换 调用格式:
[cA1,cD1]=dwt(X,‟wname‟)【使用指定的小波基函数 'wname' 对信号 X 进行分解】 [cA1,cD1]=dwt(X,Lo_D,Hi_D)【使用指定的滤波器组 Lo_D、Hi_D 对信号进行分解】 实例:
load leleccum;s=leleccum(1:3920);ls=length(s);[cA,cD]=dwt(s,'db1');figure subplot(121)plot(cA);subplot(122)plot(cD)[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db1')[cA1,cD1]=dwt(s, Lo_D, Hi_D);figure subplot(121)plot(cA1);subplot(122)plot(cD1)
idwt一维离散小波逆变换 调用格式:
X=idwt(cA,cD,'wname')【由近似分量 cA 和细节分量 cD 经小波反变换重构原始信号 X】 X=idwt(cA,cD,Lo_R, Hi_R)【用指定的重构滤波器 Lo_R 和 Hi_R 经小波反变换重构原始信号 X 】 X=idwt(cA,cD,'wname',L)【指定返回信号 X 中心附近的 L 个点】 X=idwt(cA,cD, Lo_R, Hi_R, L)【与上边一致】
upcoef一维小波系数重构 调用格式:
Y = upcoef(O,X,'wname',N)Y = upcoef(O,X,'wname',N,L)Y = upcoef(O,X,Lo_R,Hi_R,N)Y = upcoef(O,X,Lo_R,Hi_R,N,L)这一函数是一维小波分解系数的直接重构
它重建的是原信号在指定层次的,高频或者低频分量。
也就是说,这个信号不是原本的信号,而且某个层次上的逼近。N=1条件下相当与缺省模式的idwt A1 = upcoef('a',cA1,'db1',1,l_s)相当于A1 = idwt(cA1,[],'db1',l_s);D1 = upcoef('d',cD1,'db1',1,l_s)相当于 D1 = idwt([],cD1,'db1',l_s);而O是一个标识符,O=„a‟ 低频,O=„d‟ 高频 N代表了小波系数的级别,是几级小波系数。L则指定返回信号 X 中心附近的 L 个点 例子:
load leleccum;s=leleccum(1:3920);ls=length(s);[cA,cD]=dwt(s,'db1');A1=upcoef('a',cA,'db1',1,ls);D1=upcoef('d',cD,'db1',1,ls);% A1 = idwt(cA,[],'db1',ls);D1 = idwt([],cD,'db1',ls);subplot(1,2,1);plot(A1);title('Approximation A1')subplot(1,2,2);plot(D1);title('Detail D1')
多尺度一维小波变换 wavedec多尺度一维分解
调用格式:
[C, L]=wavedec(X,N,’wname’)[C, L]=wavedec(X,N,Lo_D,Hi_D)函数返回N层分解的各组分系数C(连接在一个向量里),向量L里返回的是各组分的长度。分解的结构如下
例子:
load leleccum;s=leleccum(1:3920);[C,L]=wavedec(s,3,'db1');subplot(1,2,1);plot(C);title('C')subplot(1,2,2);plot(L);title('L')
waverec多尺度一维重构 调用格式:
X=waverec(C,L,’wname’)X=waverec(C,L,Lo_R,Hi_R)实例:
load leleccum;s=leleccum(1:3920);[C,L]=wavedec(s,3,'db1');A0=waverec(C,L,'db1');Err=abs(s-A0);figure(1)subplot(3,1,1);plot(s);title('original s');subplot(3,1,2);plot(A0);title('wavereconstruction s');subplot(3,1,3);plot(Err);title('wavereconstruction err');
appcoef低频系数提取 调用格式:
A=appcoef(C,L,‟wname‟,N)从C中抽取N层近似系数 A=appcoef(C,L,‟wname‟)A=appcoef(C,L,Lo_R,Hi_R,N)A=appcoef(C,L, Lo_R,Hi_R)detcoef高频系数提取 调用格式:
A=detcoef(C,L,N)A=detcoef(C,L)实例:
load leleccum;s=leleccum(1:3920);[C,L]=wavedec(s,3,'db1');cA3=appcoef(C,L,'db1',3);cD3=detcoef(C,L,3);cD2=detcoef(C,L,2);cD1=detcoef(C,L,1);subplot(2,2,1);plot(cA3);title('cA3')subplot(2,2,2);plot(cD3);title('cD3')subplot(2,2,3);plot(cD2);title('cD2')subplot(2,2,4);plot(cD1);title('cD1')
wrcoef多尺度小波系数重构
调用格式:
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)其中type=„a‟ 低频,type=„d‟ 高频 作用类似于upcoef,只不过是多尺度函数。实例:
load leleccum;s=leleccum(1:3920);[C,L]=wavedec(s,3,'db1');cA3=wrcoef('a',C,L,'db1',3);cD1=wrcoef('d',C,L,'db1',1);cD2=wrcoef('d',C,L,'db1',2);cD3=wrcoef('d',C,L,'db1',3);subplot(2,2,1);plot(cA3);title('cA3')subplot(2,2,2);plot(cD3);title('cD3')subplot(2,2,3);plot(cD2);title('cD2')subplot(2,2,4);plot(cD1);title('cD1')
一维静态(平稳)小波变换
离散平稳小波分析所用到的函数有swt小波分解和iswt小波重构 swt一维平稳小波变换 调用格式
SWC = swt(X,N,'wname')SWC = swt(X,N,Lo_D,Hi_D)[swa,swd] = swt(X,N,'wname')[swa,swd] = swt(X,N,Lo_D,Hi_D)对于SWT变换,如果需要在第k层分解信号,那么原始信号需要能够平分成2^k份。所以如果原始信号的长度不满足要求,需要使用Signal Extension GUI工具或使用wextend函数来扩展它。
函数执行将产生N层近似和细节的系数,两者和信号的长度是相等的,这也是平稳小波和普通小波不同的地方,从而使它在某些领域有好的效果。iswt一维平稳小波逆变换 调用格式
X = iswt(SWC,'wname')X = iswt(SWA,SWD,'wname')X = iswt(SWC,Lo_R,Hi_R)X = iswt(SWA,SWD,Lo_R,Hi_R)实例
load noisdopp s = noisdopp;[swa,swd] = swt(s,1,'db1');%一层平稳小波分解
subplot(1,2,1), plot(swa);title('Approximation cfs')subplot(1,2,2), plot(swd);title('Detail cfs')res = iswt(swa,swd,'db1');%一层平稳小波重构 %从系数中构建近似和细节
nulcfs = zeros(size(swa));A1 = iswt(swa,nulcfs,'db1');D1 = iswt(nulcfs,swd,'db1');subplot(4,1,1), plot(s);title('original s');subplot(4,1,2), plot(res);title('recovered s');subplot(4,1,3), plot(A1);title('Approximation A1');subplot(4,1,4), plot(D1);title('Detail D1');
单尺度二维小波变换 dwt2二维离散小波变换
[cA1,cH1,cV1,cD1]=dwt2(X,‟wname‟)[cA1,cH1,cV1,cD1]=dwt2(X,Lo_D,Hi_D)其中cA1, cH1水平;cV1垂直;cD1对角 二维小波变换有四个值 实例:
load wbarb;figure(1);image(X);colormap(map);colorbar;[cA1,cH1,cV1,cD1]=dwt2(X,'bior3.7');figure(2);colormap(map);subplot(2,2,1);image(wcodemat(cA1,192));title('cA1')subplot(2,2,2);image(wcodemat(cH1,192));title('cH1')subplot(2,2,3);image(wcodemat(cV1,192));title('cV1')subplot(2,2,4);image(wcodemat(cD1,192));title('cD1')
idwt2二维离散小波逆变换
X = idwt2(cA1,cH1,cV1,cD1,'bior3.7');
upcoef2二维系数重构
调用格式:
Y=upcoef2(O,X,‟wname‟,N,S)Y=upcoef2(O,X,Lo_R,Hi_R,N,S)Y=upcoef2(O,X,‟wname‟,N)Y=upcoef2(O,X,Lo_R,Hi_R,N)Y=upcoef2(O,X,‟wname‟)Y=upcoef2(O,X,Lo_R,Hi_R)O:„a‟低频;„h‟水平;„v‟垂直;„d‟对角 实例: load wbarb;[cA1,cH1,cV1,cD1]=dwt2(X,'bior3.7');A1 = upcoef2('a',cA1,'bior3.7',1);H1 = upcoef2('h',cH1,'bior3.7',1);V1 = upcoef2('v',cV1,'bior3.7',1);D1 = upcoef2('d',cD1,'bior3.7',1);figure(2);colormap(map);subplot(2,2,1);image(wcodemat(A1,192));title('Approximation A1')subplot(2,2,2);image(wcodemat(H1,192));title('Horizontal Detail H1')subplot(2,2,3);image(wcodemat(V1,192));title('Vertical Detail V1')subplot(2,2,4);image(wcodemat(D1,192));title('Diagonal Detail D1')
多尺度二维小波变换 wavedec2多尺度二维分解
[C, S]=wavedec2(X,N,’wname’)[C, S]=wavedec2(X,N,Lo_D,Hi_D)
c为各层分解系数,s为各层分解系数长度,也就是大小.c的结构: c=[A(N)|H(N)|V(N)|D(N)|H(N-1)|V(N-1)|D(N-1)|H(N-2)|V(N-2)|D(N-2)|...|H(1)|V(1)|D(1)] 可见,c是一个行向量,即:1*(size(X)),(e.g,X=256*256,then c
大小为:1*(256*256)=1*65536)A(N)代表第N层低频系数,H(N)|V(N)|D(N)代表第N层高频系数,分别是水平,垂直,对角高频,以此类推,到H(1)|V(1)|D(1).s的结构: 储存各层分解系数长度的,即第一行A(N)的长度,第二行是H(N)|V(N)|D(N)|的长度,第三行是
H(N-1)|V(N-1)|D(N-1)的长度,倒数第二行是H(1)|V(1)|D(1)长度,最后一行是X的长度(大小)实例:
load wbarb;[C,S] = wavedec2(X,2,'bior3.7');figure(2);colormap(map);subplot(1,2,1);plot(C);subplot(1,2,2);plot(S);
waverec2多尺度二维重构 调用格式:
X=waverec2(C,S,’wname’)X=waverec2(C,S,Lo_R,Hi_R)实例:
load wbarb;[C,S] = wavedec2(X,2,'bior3.7');X0=waverec2(C,S,'bior3.7');Err=abs(X0-X);figure(1);colormap(map);subplot(3,1,1);image(wcodemat(X,192));title('original X')subplot(3,1,2);image(wcodemat(X0,192));title('wavereconstruction X')subplot(3,1,3);imshow(Err,[])
appcoef2低频系数提取 调用格式:
A=appcoef2(C,S,‟wname‟,N)A=appcoef2(C,S,Lo_R,Hi_R,N)detcoef2高频系数提取
调用格式:
A= detcoef2(‘type’,C,S,’wname’,N)Type: ’h’ 水平;‘v’垂直;‘d’对角 实例:
load wbarb;[C,S] = wavedec2(X,2,'bior3.7');cA2 = appcoef2(C,S,'bior3.7',2);cH2 = detcoef2('h',C,S,2);cV2 = detcoef2('v',C,S,2);cD2 = detcoef2('d',C,S,2);cA1 = appcoef2(C,S,'bior3.7',1);cH1 = detcoef2('h',C,S,1);cV1 = detcoef2('v',C,S,1);cD1 = detcoef2('d',C,S,1);figure(1);colormap(map);subplot(2,2,1);image(wcodemat(cA1,192));title('Approximation A1')subplot(2,2,2);image(wcodemat(cH1,192));title('Horizontal Detail H1')subplot(2,2,3);image(wcodemat(cV1,192));title('Vertical Detail V1')subplot(2,2,4);image(wcodemat(cD1,192));title('Diagonal Detail D1')figure(2);colormap(map);subplot(2,2,1);image(wcodemat(cA2,192));title('Approximation A2')subplot(2,2,2);image(wcodemat(cH2,192));title('Horizontal Detail H2')subplot(2,2,3);image(wcodemat(cV2,192));title('Vertical Detail V2')subplot(2,2,4);image(wcodemat(cD2,192));title('Diagonal Detail D2')
wrcoef2多尺度小波系数重构
调用格式:
X= wrcoef2(„type‟,C,S,‟wname‟,N)X= wrcoef2(„type‟,C,S,Lo_R,Hi_R,N)X= wrcoef2(„type‟,C,S,‟wname‟)X= wrcoef2(„type‟,C,S, Lo_R,Hi_R,N)Type: „a‟低频;‟h‟ 水平;„v‟垂直;„d‟对角 实例:
load wbarb;[C,S] = wavedec2(X,2,'bior3.7');A1 = wrcoef2('a',C,S,'bior3.7',1);H1 = wrcoef2('h',C,S,'bior3.7',1);V1 = wrcoef2('v',C,S,'bior3.7',1);D1 = wrcoef2('d',C,S,'bior3.7',1);A2 = wrcoef2('a',C,S,'bior3.7',2);H2 = wrcoef2('h',C,S,'bior3.7',2);V2 = wrcoef2('v',C,S,'bior3.7',2);D2 = wrcoef2('d',C,S,'bior3.7',2);figure(1);colormap(map);subplot(2,2,1);image(wcodemat(A1,192));title('reApproximation A1')subplot(2,2,2);image(wcodemat(H1,192));title('reHorizontal Detail H1')subplot(2,2,3);image(wcodemat(V1,192));title('reVertical Detail V1')subplot(2,2,4);image(wcodemat(D1,192));title('reDiagonal Detail D1')figure(2);colormap(map);subplot(2,2,1);image(wcodemat(A2,192));title('reApproximation A2')subplot(2,2,2);image(wcodemat(H2,192));title('reHorizontal Detail H2')subplot(2,2,3);image(wcodemat(V2,192));title('reVertical Detail V2')subplot(2,2,4);image(wcodemat(D2,192));title('reDiagonal Detail D2')
二维静态(平稳)小波变换
回顾从一维离散小波变换到二维的扩展,二维静态小波变换采用相似的方式。对行和列分别采用高通和低通滤波器。这样分解的结果仍然是四组图像:近似图像、水平细节图像、竖直细节图像和对角图像,与离散小波变换不同的只是静态小波分解得到的四幅图像与原图像尺寸一致,道理与一维情况相同。
二维离散小波变换同样只提供了一个函数swt2,因为它不对分解系数进行下采样,所以单层分解和多层分解的结果是一样的。不需要另外提供多层分解的功能。swt2二维静态小波变换
调用格式:
[A,H,V,D] = swt2(X,N,'wname')[A,H,V,D] = swt2(X,N,Lo_D,Hi_D)iswt2二维静态小波逆变换 调用格式
X = iswt2(SWC,'wname')X = iswt2(A,H,V,D,wname)X = iswt2(SWC,Lo_R,Hi_R)X = iswt2(A,H,V,D,Lo_R,Hi_R)实例
clc;clear;close all;load woman;[cA,cH,cV,cD]=swt2(X,2,'haar');%用haar小波基进行2尺度平稳小波分解 cA1=cA(:,:,1);cH1=cH(:,:,1);cV1=cV(:,:,1);cD1=cD(:,:,1);%尺度1低、高频系数
cA2=cA(:,:,2);cH2=cH(:,:,2);cV2=cV(:,:,2);cD2=cD(:,:,2);%尺度2低、高频系数
reX = iswt2(cA,cH,cV,cD,'haar');X1 = iswt2(cA1,cH1,cV1,cD1,'haar');X2 = iswt2(cA2,cH2,cV2,cD2,'haar');figure;colormap(map);subplot(1,4,1),image(X);title('原始图像');subplot(1,4,2),image(reX);title('重构图像');subplot(1,4,3),image(X1);title('尺度1的重构图像');subplot(1,4,4),image(X2);title('尺度1的重构图像');figure;colormap(map);subplot(1,2,1),image(uint8(cA1));title('尺度1的低频系数图像');subplot(1,2,2),image(uint8(cA2));title('尺度2的低频系数图像');figure;colormap(map);subplot(2,3,1),image(uint8(cH1));title('尺度1水平方向高频系数');subplot(2,3,2),image(uint8(cV1));title('尺度1垂直方向高频系数');subplot(2,3,3),image(uint8(cD1));title('尺度1斜线方向高频系数');subplot(2,3,4),image(uint8(cH2));title('尺度2水平方向高频系数');subplot(2,3,5),image(uint8(cV2));title('尺度2垂直方向高频系数');subplot(2,3,6),image(uint8(cD2));title('尺度2斜线方向高频系数');
直接调用的小波函数 meyer函数 功能:Meyer小波
函数语法格式:
[PHI,PSI,T] = meyer(LB,UB,N)[PHI,T] = meyer(LB,UB,N,'phi')[PSI,T] = meyer(LB,UB,N,'psi')
cgauwavf函数
功能:Complex Gaussian小波
函数语法格式:
[PSI,X] = cgauwavf(LB,UB,N,P)
mexihat函数 功能:墨西哥帽小波
函数语法格式:
[PSI,X] = mexihat(LB,UB,N)
morlet函数 功能:Morlet小波
函数语法格式:
[PSI,X] = morlet(LB,UB,N)
symwavf函数
功能:Symlets小波滤波器 函数语法格式:F = symwavf(W)
三、图像接口调用
使用图形接口做一维连续小波分析
1.开启一维连续小波工具,只需输入如下命令 wavemenu
出现如下小波工具箱主菜单
选择Continuous Wavelet 1-D菜单项,出现如下一维信号分析连续小波分析工具
2.加载信号
选择菜单File->Load Signal,在Load Signal对话框里选择noissin.mat文件,它在matlab安装目录的toolbox/wavelet/wavedemo文件夹下,点击OK加载信号。一维连续小波工具开始加载信号,加载后默认采样频率为1s。3.执行连续小波变换
下面来测试使用db4小波对尺度1到48做小波分析,设置如下
4.点击Analyze按钮
在短暂的计算后,工具将绘制小波系数,并在Coefficients line坐标系中绘制尺度为24的小波系数,在local maxima坐标系中绘制各尺度的小波系数最大值。
5.查看小波Coefficients Line 在小波系数图中右键点击可以选择展示其他尺度的小波系数,选择后点击New Coefficients Line按钮,Coefficients Line会相应更新。
6.查看Maxima Line 点击Refresh Maxima Line按钮,可以显示从尺度1到所选尺度的小波系数的最大值。
注意当在系数图中按下鼠标右键并移动时,会在最下面的Info框中显示当前鼠标位于的X位置和尺度。
7.在尺度和伪频率之间切换
在右边选择Frequencies,当再在系数图中选择时,在Info中显示的将是Hz。
而关于尺度和频率的转换关系,可以看How to Connect Scale to Frequency? 8.选择要显示的坐标系
9.放大细节
在系数框中按鼠标左键可以选择放大的范围。
10.选择好放大范围后点击最下面的按钮可以实现指定的放大
11.显示普通系数或系数绝对值
两种显示方式的区别在于,普通模式下,颜色映射是在系数的最大最小之间;而绝对模式,颜色映射是在0和最大的系数绝对值之间。图形接口的导入导出信息 导入信号到一维连续小波工具
首先将要处理的信号保存到mat文件中,要求信号是一维的向量。然后使用工具的File->Load Signal菜单功能,选择此信号文件即可导入信号。
文件中第一个一维变量被认为是信号,变量在文件中顺序是按字母排序的。保存小波系数
小波分析完成后,点击File->Save->Coefficients,可以将分析结果保存到mat文件。保存后,可以使用load函数加载数据,会看到保存的变量有小波系数coeff、尺度scales、小波的名字wname。
使用图形接口做一维离散小波分析 1.开启一维小波分析工具 Wavemenu->Wavelet 1-D 2.加载信号
3.执行一层小波分解 使用db1小波执行一层分解
4.放大有关细节 5.执行多层小波分解 使用db1小波执行3层分解。
选择不同的显示方式:在Display mode下拉菜单下可以选择不同的显示方式,默认的显示方式为Full Decomposition Mode,其他的显示方式及其意义如下 Separate Mode:在不同的列中显示细节和近似;
Superimpose Mode:在一张图上以不同的颜色显示细节、近似;
Tree Mode:显示分解树、原始信号和选择的成分,在分解树上选择你想显示的成分; Show and Scroll Mode:显示3个窗口,第一个显示原始信号和选择的近似信号,第二个显示选择的细节,第三个显示小波系数; Show and Scroll Mode(Stem Cfs):和Show and Scroll Mode很接近,除了第三个窗口中以杆状图替代颜色条显示小波系数。
对于每个分析任务,可以改变默认的显示方式,只要在View->Default Display Mode子菜单下选择理想的方式即可;不同的显示方式会有额外的显示选项,在More Display Options中做选择,这些选项可以控制不同成分的显示、选择是否显示原始信号与细节、近似对比。6.从信号中移除噪声
图形接口提供了以预定义的阈值策略除噪的选项,这使得从信号中除噪非常容易,只需点击De-noise按钮就可以弹出除噪工具。
点击Close可以关闭除噪窗口,由于不能同时打开除噪和压缩窗口,所以需要关闭除噪窗口再进行信号压缩。关闭时会提示Update Synthesized Signal提示对话框,点击No,如果点击Yes,合成的信号会加载到主窗口。7.改善分析
图像工具可以在任何时候轻易的改善分析,只需要改变分析的方法就可以了,如使用db3做5层小波分析。8.压缩信号
图形接口提供了自动化或人工阈值的做压缩的功能。
默认使用的是全阈值方法,当然也可以使用人工阈值的方法,选择By Level thresholding选项即可,下面的滑动条提供了各级阈值独立调整的功能,相应的调整可以在左边的窗口中看到,在图形窗口中也可以直接拖动来改变阈值。
完成选择后,点击压缩按钮可以完成压缩。从压缩的结果可以看到,压缩过程去除了大多数噪声,但保存了信号99.74%的能量。自动化阈值是非常有效的,它使除3.2%的小波系数都归零化了。9.显示残差
点击Residuals按钮可以查看压缩的残差。显示的统计数据包括测量的趋势(平均值、众数、中值)和散布情况(极差、标准差)。另外,工具还提供了概率分布直方图和累计直方图以及时间序列图,如自相关函数、频谱。这些都是和去噪工具是一样的。10.显示统计分布
可以显示一系列有关信号及其组分的统计数据。
点击Statistics按钮可以查看统计数据信息,点击Histograms可以查看直方图。从图形接口中导入导出信息 保存信息
l 保存合成的信号
如加载如下信号File > Example Analysis > Basic Signals > with db3 at level 5 → Sum of sines,做除噪或压缩处理后,保存合成信号File > Save > Synthesized Signal,保存后加载文件,会得到如下变量: 如果使用除全阈值外的方法时,得到的变量结构如下
Synthsig是合成的信号,除噪或压缩的小波方法保存在wname中,相互依赖的各级阈值保存在thrParams中,小波分解的等级数和cell的长度相等,thrParams{i},i从1到5分别保存了阈值间距上下限的值和阈值(间距阈值是允许的,在自适应阈值方法中会用到,参见One-Dimensional Variance Adaptive Thresholding of Wavelet Coefficients)如果使用的全阈值方法,保存的数据结构如下
alTHR保存的是全阈值的值。l 保存离散小波变换的系数 一个例子的文件内容如下
Coefs包含了离散小波变换的系数,longs包含了各组分的长度,thrParams为空,因为合成信号不存在,wname是小波的名字。l 保存分解结果(即保存小波分析的全体数据)小波工具将保存为.wal文件,加载方式为 load wdecex1d.wa1 –mat 文件内容为
加载信息
加载的文件只要和保存的相应文件中的变量一样即可。
使用图形接口分析复信号
与实信号不同的是,选择Complex Continuous Wavelet 1-D,得到的结果如下
具体操作过程与实信号的相似,如下
使用图形接口做一维除噪分析 1.开启一维平稳小波除噪工具
输入wavemenu,选择SWT Denoising 1-D,出现如下GUI
2.加载信号 Load Signal 3.执行平稳小波分解
使用db1小波执行5层小波分解,得到的是非抽取系数(nondecimated coefficients),它们是使用相同的离散小波变换来得到的,只是省略了抽取的步骤。得到的结果如下
4.使用平稳小波变换除噪
可以使用GUI默认的参数做除噪处理。右边的滑动条可以控制各级系数的阈值大小,也可以直接在系数图中直接拖动来调整阈值的大小,注意近似系数中没有阈值。点击denoise进行除噪处理
得到的效果是非常好的,但似乎在信号不连续的地方出现了过平滑,这个可以从残差图中看出来,在800的位置出现了突降点。
选择hard阈值模式代替soft模式,再进行除噪,结果如下
可以看到这次效果非常好,而且残差图也看起来像白噪声序列。为了验证这一点,可以点击Residuals按钮查看残差图及相关统计数据来详细说明。
四、小波变换在图像处理中的应用 4.1 小波分析用于图像压缩
4.1.1 基于小波变换的图像局部压缩
基于离散余弦变换的图像压缩算法,其基本思想是在频域对信号进行分解,驱除信号点之间的相关性,并找出重要系数,滤掉次要系数,以达到压缩的效果,但该方法在处理过程中并不能提供时域的信息,在我们比较关心时域特性的时候显得无能为力。
但是这种应用的需求是很广泛的,比如遥感测控图像,要求在整幅图像有很高压缩比的同时,对热点部分的图像要有较高的分辨率,例如医疗图像,需要对某个局部的细节部分有很高的分辨率,单纯的频域分析的方法显然不能达到这个要求,虽然可以通过对图像进行分快分解,然后对每块作用不同的阈值或掩码来达到这个要求,但分块大小相对固定,有失灵活。
在这个方面,小波分析的就优越的多,由于小波分析固有的时频特性,我们可以在时频两个方向对系数进行处理,这样就可以对我们感兴趣的部分提供不同的压缩精度。
下面我们利用小波变化的时频局部化特性,举一个局部压缩的例子,大家可以通过这个例子看出小波变换在应用这类问题上的优越性。
我们把图像中部的细节系数都置零,从压缩图像中可以很明显地看出只有中间部分变得模糊(比如在原图中很清晰的围巾的条纹不能分辨),而其他部分的细节信息仍然可以分辨的很清楚。
实例:
load wbarb %使用sym4小波对信号进行一层小波分解 [ca1,ch1,cv1,cd1]=dwt2(X,'sym4');codca1=wcodemat(ca1,192);codch1=wcodemat(ch1,192);codcv1=wcodemat(cv1,192);codcd1=wcodemat(cd1,192);%将四个系数图像组合为一个图像
codx=[codca1,codch1,codcv1,codcd1];%复制原图像的小波系数 rca1=ca1;rch1=ch1;rcv1=cv1;rcd1=cd1;%将三个细节系数的中部置零
rch1(33:97,33:97)=zeros(65,65);rcv1(33:97,33:97)=zeros(65,65);rcd1(33:97,33:97)=zeros(65,65);codrca1=wcodemat(rca1,192);codrch1=wcodemat(rch1,192);codrcv1=wcodemat(rcv1,192);codrcd1=wcodemat(rcd1,192);%将处理后的系数图像组合为一个图像
codrx=[codrca1,codrch1,codrcv1,codrcd1];%重建处理后的系数
rx=idwt2(rca1,rch1,rcv1,rcd1,'sym4');subplot(221);image(wcodemat(X,192)),colormap(map);title('原始图像');subplot(222);image(codx),colormap(map);title('一层分解后各层系数图像');subplot(223);image(wcodemat(rx,192)),colormap(map);title('压缩图像');subplot(224);image(codrx),colormap(map);title('处理后各层系数图像');%求压缩信号的能量成分 per=norm(rx)/norm(X)%求压缩信号与原信号的标准差 err=norm(rx-X)运行结果
per = 1.0000 err = 586.4979 4.1.2 小波变换用于图像压缩的一般方法
二维小波分析用于图像压缩是小波分析应用的一个重要方面。
它的特点是压缩比高,压缩速度快,压缩后能保持图像的特征基本不变,且在传递过程中可以抗干扰。4.1.2.1 利用二维小波分析进行图像压缩
基于小波分析的图像压缩方法很多,比较成功的有小波包、小波变换零树压缩、小波变换矢量量化压缩等。
利用二维小波变换进行图像压缩时,小波变换将图像从空间域变换到时间域,它的作用与以前在图像压缩中所用到的离散余弦(DCT)、傅立叶变换(FFT)等的作用类似。但是要很好的进行图像的压缩,需要综合的利用多种其他技术,特别是数据的编码与解码算法等,所以利用小波分析进行图像压缩通常需要利用小波分析和许多其他相关技术共同完成。小波高频滤波压缩 %%%%%%%%装入图像、变换 load wbarb;%对图像用bior3.7小波进行2层小波分解 [c,s]=wavedec2(X,2,'bior3.7');%提取小波分解结构中第一层低频系数和高频系数 ca1=appcoef2(c,s,'bior3.7',1);ch1=detcoef2('h',c,s,1);cv1=detcoef2('v',c,s,1);cd1=detcoef2('d',c,s,1);%分别对各频率成分进行重构
a1=wrcoef2('a',c,s,'bior3.7',1);h1=wrcoef2('h',c,s,'bior3.7',1);v1=wrcoef2('v',c,s,'bior3.7',1);d1=wrcoef2('d',c,s,'bior3.7',1);c1=[a1,h1;v1,d1];%%%%%%%%%下面进行图像压缩处理
%保留小波分解第一层低频信息,进行图像的压缩 %第一层的低频信息即为ca1,显示第一层的低频信息 %首先对第一层信息进行量化编码
ca1=appcoef2(c,s,'bior3.7',1);ca1=wcodemat(ca1,440,'mat',0);%改变图像的高度 ca1=0.5*ca1;%保留小波分解第二层低频信息,进行图像的压缩,此时压缩比更大 %第二层的低频信息即为ca2,显示第二层的低频信息 ca2=appcoef2(c,s,'bior3.7',2);%首先对第二层信息进行量化编码
ca2=wcodemat(ca2,440,'mat',0);%改变图像的高度 ca2=0.25*ca2;%%%%%%%显示图像、分解后各频率成分的信息
subplot(221);image(X);colormap(map);title('原始图像');axis square subplot(222);image(c1);axis square;title('分解后低频和高频信息');subplot(223);image(ca1);colormap(map);axis square;title('第一次压缩');subplot(224);image(ca2);colormap(map);axis square;title('第二次压缩');disp('压缩前图像X的大小:');whos('X')disp('第一次压缩图像的大小为:');whos('ca1')disp('第二次压缩图像的大小为:');whos('ca2')输出结果如下所示:
压缩前图像X的大小: Name SizeBytesClass X256x256524288double array Grand total is 65536 elements using 524288 bytes 第一次压缩图像的大小为: Name SizeBytesClass ca1135x135145800double array Grand total is 18225 elements using 145800 bytes 第二次压缩图像的大小为: Name SizeBytesClass ca2 75x7545000double array Grand total is 5625 elements using 45000 bytes
可以看出,第一次压缩提取的是原始图像中小波分解第一层的低频信息,此时压缩效果较好,压缩比较小(约为1/3):第二次压缩是提取第一层分解低频部分的低频部分(即小波分解第二层的低频部分),其压缩比较大(约为1/12),压缩效果在视觉上也基本过的去。这是一种最简单的压缩方法,只保留原始图像中低频信息,不经过其他处理即可获得较好的压缩效果。
小波高频阈值压缩 %装入一个二维信号 load tire;
%下面进行图像压缩
[c,s]=wavedec2(X,2,'db3');%对图像用db3小波进行2层小波分解, [thr,sorh,keepapp]=ddencmp('cmp','wv',X);%使用ddencmp函数来求取阈值
[Xcomp,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'db3',2,thr,sorh,keepapp);%输入参数中选择了全局阈值选项„gbl‟,用来对所有高频系数进行相同的阈值量化处理 %将压缩后的图像与原始图像相比较,并显示出来
subplot(121);image(X);colormap(map);title('原始图像');axis square subplot(122);image(Xcomp);colormap(map);title('压缩图像');axis square disp('小波分解系数中置0的系数个数百分比:');perf0 disp('压缩后图像剩余能量百分比:');perfl2
输出结果如下所示:
小波分解系数中置0的系数个数百分比: perf0 =49.1935 压缩后图像剩余能量百分比: perfl2 =99.9928 4.1.2.2 二维信号压缩中的阈值的确定与作用命令
由于阈值处理只关心系数的绝对值,并不关心系数的位置,所以二维小波变换系数的阈值化方法同一维情况大同小异,为了方便用户使用小波工具箱对某些阈值化方法提供了专门的二维处理命令
阈值化压缩方法 %装入一个二维信号 load tire;%下面进行图像压缩
[c,s]=wavedec2(X,2,'db3');%对图像用db3小波进行2层小波分解, [thr,sorh,keepapp]=ddencmp('cmp','wv',X);%使用ddencmp函数来求取阈值 [Xcomp,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'db3',2,thr,sorh,keepapp);%将压缩后的图像与原始图像相比较,并显示出来 subplot(121);image(X);colormap(map);title('原始图像');axis square subplot(122);image(Xcomp);colormap(map);title('压缩图像');axis square disp('小波分解系数中置0的系数个数百分比:');perf0 disp('压缩后图像剩余能量百分比:');perfl2
显示结果如图所示
可见分层阈值化压缩方法同全局阈值化方法相比,在能量损失不是很大的情况下可以获得最高的压缩化,这主要是因为层数和方向相关的阈值化方法能利用更精细的细节信息进行阈值化处理。
4.1.3 基于小波包变换的图像压缩
小波分析之所以在信号处理中有着强大的功能,是基于其分离信息的思想,分离到各个小波域的信息除了与其他小波域的关联,使得处理的时候更为灵活。全局阈值化方法作用的信息粒度太大,不够精细,所以很难同时获得高的压缩比和能量保留成分,在作用的分层阈值以后,性能明显提高,因为分层阈值更能体现信号固有的时频局部特性。
但是小波分解仍然不够灵活,分解出来的小波树只有一种模式,不能完全地体现时频局部化信息。而压缩的核心思想既是尽可能去处各小波域系数之间的信息关联,最大限度体现时频局部化的信息,因此,实际的压缩算法多采用小波包算法,而小波树的确定则是根据不同的信息论准则,以达到分解系数表达的信息密度最高。
下面我通过一个例子来说明小波包分析在图像压缩中的应用,并给出性能参数以便于同基于小波分析的压缩进行比较。
小波包优化 %读入信号 load julia %求颜色索引表长度 nbc=size(map,1);%得到信号的阈值,保留层数,小波树优化标准
[thr,sorh,keepapp,crit]=ddencmp('cmp','wp',X);%通过以上得到的参数对信号进行压缩
[xd,treed,perf0,perfl2]=wpdencmp(X,sorh,4,'sym4',crit,thr*2,keepapp);%更改索引表为pink索引表 colormap(pink(nbc));subplot(121);image(wcodemat(X,nbc));title('原始图像');subplot(122);image(wcodemat(xd,nbc));title('全局阈值化压缩图像');xlabel(['能量成分',num2str(perfl2),'%','零系数成分',num2str(perf0),'%']);plot(treed);
得到的压缩结果如图所示
图6 基于小波包分析的图像压缩
压缩过程中使用的最优小波数如图所示
图7 最优小波树
ddencmp、wpdencmp这两个命令是Matlab小波工具箱提供的自动获取阈值和自动使用小波包压缩的命令,后者将分解阈值化和重建综合起来。在将小波包用于信号压缩的过程中,ddencmp命令返回的最优小波树标准都是阈值化标准。根据这个标准确定的最优小波树可以使得压缩过程的零系数成分最高,并且自动降低计算量。
图像压缩是应用非常广泛的一类问题,所以其机器实现效率是至关重要的,在实际的应用中,如JPEG2000,一般不采用通常的mallat算法做小波分解,而是应用特定的双正交小波,利用其滤波器分布规则的特性,用移位操作来实现滤波操作。4.2 小波分析用于图像去噪
噪声可以理解为妨碍人的视觉器官或系统传感器对所接收图像源进行理解或分析的各种因素。一般噪声是不可预测的随机信号,它只能用概率统计的方法去认识。噪声对图像处理十分重要,它影响图像处理的输入、采集、处理的各个环节以及输出结果的全过程。特别是图像的输入、采集的噪声是个十分关键的问题,若输入伴有较大噪声,必然影响处理全过程及输出结果。因此一个良好的图像处理系统,不论是模拟处理还是计算机处理无不把减少最前一级的噪声作为主攻目标。去噪已成为图像处理中极其重要的步骤。
对二维图像信号的去噪方法同样适用于一维信号,尤其是对于几何图像更适合。二维模型可以表述为
s(i,j)=f(i,j)+δ·e(i,j)i,j=0,1,…,m-1 其中,e是标准偏差不变的高斯白噪声。二维信号用二维小波分析的去噪步骤有3步:(1)二维信号的小波分解。选择一个小波和小波分解的层次N,然后计算信号s到第N层的分解。
(2)对高频系数进行阈值量化。对于从1到N的每一层,选择一个阈值,并对这一层的高频系数进行软阈值量化处理。
(3)二维小波的重构。根据小波分解的第N层的低频系数和经过修改的从第一层到第N层的各层高频系数计算二维信号的小波重构。
在这3个步骤中,重点是如何选取阈值和阈值的量化 小噪声阈值去噪
下面给出一个二维信号(文件名为detfinger.mat),并利用小波分析对信号进行去噪处理。
%装入图像 load tire %下面进行噪声的产生 init=3718025452;rand('seed',init);Xnoise=X+18*(rand(size(X)));%用sym5小波对图像信号进行二层的小波分解 [c,s]=wavedec2(X,2,'sym5');%下面进行图像的去噪处理
%使用ddencmp函数来计算去噪的默认阈值和熵标准
%使用wdencmp函数来实现图像的压缩
[thr,sorh,keepapp]=ddencmp('den','wv',Xnoise);[Xdenoise,cxc,lxc,perf0,perfl2]=wdencmp('gbl',c,s,'sym5',2,thr,sorh,keepapp);%显示去噪后的图像 colormap(map);subplot(131);image(wcodemat(X,192));title('原始图像X');axis square subplot(132);image(wcodemat(X,192));title('含噪声的图像Xnoise');axis square subplot(133);image(Xdenoise);title('去噪后的图像');axis square
输出结果
从图中3个图像的比较可以看出,Matlab中的ddencmp和wdencmp函数可以有效地进行去噪处理。大噪声滤波去噪
再给定一个有较大白噪声的delmontl.mat图像。由于图像所含的噪声主要是白噪声,而且主要集中在图像的高频部分,所以我们可以通过全部滤掉图像中的高频部分实现图像的去噪。
load wmandril;%产生含噪图像
init=2055615866;randn('seed',init)x=X+38*randn(size(X));
%下面进行图像的去噪处理
%用小波函数sym4对x进行2层小波分解 [c,s]=wavedec2(x,2,'sym4');%提取小波分解中第一层的低频图像,即实现了低通滤波去噪 a1=wrcoef2('a',c,s,'sym4',1);%提取小波分解中第二层的低频图像,即实现了低通滤波去噪 a2=wrcoef2('a',c,s,'sym4',2);
%画出图像
subplot(221);image(X);colormap(map);title('原始图像');axis square subplot(222);image(x);colormap(map);title('含噪声图像');axis square;subplot(223);image(a1);title('第一次去噪图像');axis square;subplot(224);image(a2);title('第二次去噪图像');axis square;输出结果如图:
从上面的输出结果可以看出,第一次去噪已经滤去了大部分的高频噪声,但从去噪图像与原始图像相比可以看出,第一次去噪后的图像中还是含有不少的高频噪声;第二次去噪是在第一次去噪的基础上,再次滤去其中的高频噪声。从去噪的结果可以看出,它具有较好的去噪效果。
少量噪声的小波分解系数阈值量化去噪
下面再给出定一个喊有较少噪声的facets.mat图像。由于原始图像中只喊有较少的高频噪声,如果按照上一个例子把高频噪声全部滤掉的方法将损坏图像中固有的高频有用信号。因此这幅图像适合采用小波分解系数阈值量化方法进行去噪处理。
%下面装入原始图像,X中含有被装载的图像 load facets;%产生含噪声图像
init=2055615866;randn('seed',init)x=X+10*randn(size(X));%下面进行图像的去噪处理
[c,s]=wavedec2(x,2,'coif3');%用小波画数coif3对x进行2层小波分解 n=[1,2];%设置尺度向量n p=[10.12,23.28];%设置阈值向量p %对三个方向高频系数进行阈值处理 nc=wthcoef2('h',c,s,n,p,'s');nc=wthcoef2('v',nc,s,n,p,'s');nc=wthcoef2('d',nc,s,n,p,'s');%对新的小波分解结构[nc,s]进行重构 xx=waverec2(nc,s,'coif3');%画出图像
subplot(131);image(X);colormap(map);title('原始图像');axis square subplot(132);image(x);colormap(map);title('含噪声图像');axis square subplot(133);image(xx);colormap(map);title('去噪后的图像');axis square
输出结果如图
更加标准的做法是下面这样:
第三篇:MATLAB函数总结
MATLAB函数总结
Matlab有没有求矩阵行数/列数/维数的函数? ndims(A)返回A的维数
size(A)返回A各个维的最大元素个数 length(A)返回max(size(A))[m,n]=size(A)如果A是二维数组,返回行数和列数 nnz(A)返回A中非0元素的个数
MATLAB的取整函数:fix(x), floor(x):,ceil(x), round(x)(1)fix(x): 截尾取整.>> fix([3.12-3.12])ans = 3-3(2)floor(x):不超过x 的最大整数.(高斯取整)
>> floor([3.12-3.12])ans = 3-4
(3)ceil(x): 大于x 的最小整数
>> ceil([3.12-3.12])ans = 4-3
(4)四舍五入取整
>> round(3.12-3.12)ans = 0
>> round([3.12-3.12])ans = 3-3 >>
如何用matlab生成随机数函数 rand(1)rand(n):生成0到1之间的n阶随机数方阵 rand(m,n):生成0到1之间的m×n的随机数矩阵(现成的函数)另外:
Matlab随机数生成函数
betarnd 贝塔分布的随机数生成器 binornd 二项分布的随机数生成器 chi2rnd 卡方分布的随机数生成器 exprnd 指数分布的随机数生成器 frnd f分布的随机数生成器 gamrnd 伽玛分布的随机数生成器 geornd 几何分布的随机数生成器 hygernd 超几何分布的随机数生成器 lognrnd 对数正态分布的随机数生成器 nbinrnd 负二项分布的随机数生成器 ncfrnd 非中心f分布的随机数生成器 nctrnd 非中心t分布的随机数生成器 ncx2rnd 非中心卡方分布的随机数生成器 normrnd 正态(高斯)分布的随机数生成器 poissrnd 泊松分布的随机数生成器 raylrnd 瑞利分布的随机数生成器 trnd 学生氏t分布的随机数生成器 unidrnd 离散均匀分布的随机数生成器 unifrnd 连续均匀分布的随机数生成器 weibrnd 威布尔分布的随机数生成器
一、MATLAB常用的基本数学函数
abs(x):纯量的绝对值或向量的长度
angle(z):复数z的相角(Phase angle)
sqrt(x):开平方
real(z):复数z的实部
imag(z):复数z的虚部
conj(z):复数z的共轭复数
round(x):四舍五入至最近整数
fix(x):无论正负,舍去小数至最近整数
floor(x):地板函数,即舍去正小数至最近整数
ceil(x):天花板函数,即加入正小数至最近整数
rat(x):将实数x化为分数表示
rats(x):将实数x化为多项分数展开
sign(x):符号函数(Signum function)。
当x<0时,sign(x)=-1;
当x=0时,sign(x)=0;
当x>0时,sign(x)=1。
rem(x,y):求x除以y的馀数
gcd(x,y):整数x和y的最大公因数
lcm(x,y):整数x和y的最小公倍数
exp(x):自然指数
pow2(x):2的指数
log(x):以e为底的对数,即自然对数或
log2(x):以2为底的对数
log10(x):以10为底的对数
二、MATLAB常用的三角函数
sin(x):正弦函数
cos(x):馀弦函数
tan(x):正切函数
asin(x):反正弦函数
acos(x):反馀弦函数
atan(x):反正切函数
atan2(x,y):四象限的反正切函数
sinh(x):超越正弦函数
cosh(x):超越馀弦函数
tanh(x):超越正切函数
asinh(x):反超越正弦函数
acosh(x):反超越馀弦函数
atanh(x):反超越正切函数
三、适用於向量的常用函数有:
min(x): 向量x的元素的最小值
max(x): 向量x的元素的最大值
mean(x): 向量x的元素的平均值
median(x): 向量x的元素的中位数
std(x): 向量x的元素的标准差
diff(x): 向量x的相邻元素的差
sort(x): 对向量x的元素进行排序(Sorting)
length(x): 向量x的元素个数
norm(x): 向量x的欧氏(Euclidean)长度
sum(x): 向量x的元素总和
prod(x): 向量x的元素总乘积
cumsum(x): 向量x的累计元素总和
cumprod(x): 向量x的累计元素总乘积
dot(x, y): 向量x和y的内积
cross(x, y): 向量x和y的外积
四、MATLAB的永久常数
i或j:基本虚数单位(即)
eps:系统的浮点(Floating-point)精确度
inf:无限大,例如1/0
nan或NaN:非数值(Not a number),例如0/0
pi:圆周率 p(= 3.1415926...)
realmax:系统所能表示的最大数值
realmin:系统所能表示的最小数值
nargin: 函数的输入引数个数
nargin: 函数的输出引数个数
五、MATLAB基本绘图函数
plot: x轴和y轴均为线性刻度(Linear scale)
loglog: x轴和y轴均为对数刻度(Logarithmic scale)
semilogx: x轴为对数刻度,y轴为线性刻度
semilogy: x轴为线性刻度,y轴为对数刻度
六、plot绘图函数的叁数
字元颜色字元图线型态
y 黄色.点
k 黑色 o 圆
w 白色 x x
b 蓝色 + +
g 绿色 * *
r 红色-实线
c 亮青色 : 点线
m 锰紫色-.点虚线
--虚线
七、注解
xlabel('Input Value');% x轴注解
ylabel('Function Value');% y轴注解
title('Two Trigonometric Functions');% 图形标题
legend('y = sin(x)','y = cos(x)');% 图形注解
grid on;% 显示格线八、二维绘图函数
bar 长条图
errorbar 图形加上误差范围
fplot 较精确的函数图形
polar 极座标图
hist 累计图
rose 极座标累计图
stairs 阶梯图
stem 针状图
fill 实心图
feather 羽毛图
compass 罗盘图
quiver 向量场图
---------------------------- 附录1 常用命令
附录1.1 管理用命令函数名功能描述函数名功能描述
addpath 增加一条搜索路径 rmpath 删除一条搜索路径
demo 运行Matlab演示程序 type 列出.M文件
doc 装入超文本文档 version 显示Matlab的版本号
help 启动联机帮助 what 列出当前目录下的有关文件
lasterr 显示最后一条信息 whatsnew 显示Matlab的新特性
lookfor 搜索关键词的帮助 which 造出函数与文件所在的目录
path 设置或查询Matlab路径
附录1.2管理变量与工作空间用命令函数名功能描述函数名功能描述
clear 删除内存中的变量与函数 pack 整理工作空间内存
disp 显示矩阵与文本 save 将工作空间中的变量存盘
length 查询向量的维数 size 查询矩阵的维数
load 从文件中装入数据 who,whos 列出工作空间中的变量名
附录1.3文件与操作系统处理命令函数名功能描述函数名功能描述
cd 改变当前工作目录 edit 编辑.M文件
delete 删除文件 matlabroot 获得Matlab的安装根目录
diary 将Matlab运行命令存盘 tempdir 获得系统的缓存目录
dir 列出当前目录的内容 tempname 获得一个缓存(temp)文件
!执行操作系统命令
附录1.4窗口控制命令函数名功能描述函数名功能描述
echo 显示文件中的Matlab中的命令 more 控制命令窗口的输出页面
format 设置输出格式
附录1.5启动与退出命令函数名功能描述函数名功能描述
matlabrc 启动主程序 quit 退出Matlab环境
startup
Matlab自启动程序
附录2 运算符号与特殊字符附录
2.1运算符号与特殊字符函数名功能描述函数名功能描述
+ 加...续行标志
-减 , 分行符(该行结果不显示)
* 矩阵乘;分行符(该行结果显示)
.* 向量乘 % 注释标志
^ 矩阵乘方!操作系统命令提示符
.^ 向量乘方 ' 矩阵转置
kron 矩阵kron积.向量转置
矩阵左除 = 赋值运算
/ 矩阵右除 == 关系运算之相等
.向量左除 ~= 关系运算之不等
./ 向量右除<关系运算之小于
: 向量生成或子阵提取<= 关系运算之小于等于
()下标运算或参数定义>关系运算之大于
[] 矩阵生成>= 关系运算之大于等于
{} &逻辑运算之与
.结构字段获取符 | 逻辑运算之或
.点乘运算,常与其他运算符联合使用(如.)~ 逻辑运算之非
xor 逻辑运算之异成附录2.2逻辑函数函数名功能描述函数名功能描述
all 测试向量中所用元素是否为真 is*(一类函数)
检测向量状态.其中*表示一个确定的函数(isinf)
any 测试向量中是否有真元素 *isa 检测对象是否为某一个类的对象
exist 检验变量或文件是否定义 logical 将数字量转化为逻辑量
find 查找非零元素的下标
附录3 语言结构与调试
附录3.1编程语言函数名功能描述函数名功能描述
builtin 执行Matlab内建的函数 global 定义全局变量
eval 执行Matlab语句构成的字符串 nargchk 函数输入输出参数个数检验
feval 执行字符串指定的文件 script Matlab语句及文件信息
function Matlab函数定义关键词
附录3.2控制流程函数名功能描述函数名功能描述
break 中断循环执行的语句 if 条件转移语句
case 与switch结合实现多路转移 otherwise 多路转移中的缺省执行部分
else 与if一起使用的转移语句 return 返回调用函数
elseif 与if一起使用的转移语句 switch 与case结合实现多路转移
end 结束控制语句块 warning 显示警告信息
error 显示错误信息 while 循环语句
for 循环语句
附录3.3交互输入函数名功能描述函数名功能描述
input 请求输入 menu 菜单生成
keyboard 启动键盘管理 pause 暂停执行
附录3.4面向对象编程函数名功能描述函数名功能描述
class 生成对象 isa 判断对象是否属于某一类
double 转换成双精度型 superiorto 建立类的层次关系
inferiorto 建立类的层次关系 unit8 转换成8字节的无符号整数
inline 建立一个内嵌对象
附录3.5调试函数名功能描述函数名功能描述
dbclear 清除调试断点 dbstatus 列出所有断点情况
dbcont 调试继续执行 dbstep 单步执行
dbdown 改变局部工作空间内存 dbstop 设置调试断点
dbmex 启动对Mex文件的调试 sbtype 列出带命令行标号的.M文件
dbquit 退出调试模式 dbup 改变局部工作空间内容
dbstack 列出函数调用关系
附录4 基本矩阵与矩阵处理
附录4.1基本矩阵函数名功能描述函数名功能描述
eye 产生单位阵 rand 产生随机分布矩阵
linspace 构造线性分布的向量 randn 产生正态分布矩阵
logspace 构造等对数分布的向量 zeros 产生零矩阵
ones 产生元素全部为1的矩阵 : 产生向量
附录4.2特殊向量与常量函数名功能描述函数名功能描述
ans 缺省的计算结果变量 non 非数值常量常由0/0或Inf/Inf获得
computer 运行Matlab的机器类型 nargin 函数中参数输入个数
eps 精度容许误差(无穷小)nargout 函数中输出变量个数
flops 浮点运算计数 pi 圆周率
i 复数单元 realmax 最大浮点数值
inf 无穷大 realmin 最小浮点数值
inputname 输入参数名 varargin 函数中输入的可选参数
j 复数单元 varargout 函数中输出的可选参数
附录4.3时间与日期函数名功能描述函数名功能描述
calender 日历 eomday 计算月末
clock 时钟 etime 所用时间函数
cputime 所用的CPU时间 now 当前日期与时间
date 日期 tic 启动秒表计时器
datenum 日期(数字串格式)toc 读取秒表计时器
datestr 日期(字符串格式)weekday 星期函数
datevoc 日期(年月日分立格式)
附录4.4矩阵处理函数名功能描述函数名功能描述
cat 向量连接 reshape 改变矩阵行列个数
diag 建立对角矩阵或获取对角向量 rot90 将矩阵旋转90度
fliplr 按左右方向翻转矩阵元素 tril 取矩阵的下三角部分
flipud 按上下方向翻转矩阵元素 triu 取矩阵的上三角部分
repmat 复制并排列矩阵函数
附录5 特殊矩阵函数名功能描述函数名功能描述
compan 生成伴随矩阵 invhilb 生成逆hilbert矩阵
gallery 生成一些小的测试矩阵 magic 生成magic矩阵
hadamard 生成hadamard矩阵 pascal 生成pascal矩阵
hankel 生成hankel矩阵 toeplitz 生成toeplitz矩阵
hilb 生成hilbert矩阵 wilkinson 生成wilkinson特征值测试矩阵
附录6 数学函数
附录6.1三角函数函数名功能描述函数名功能描述
sin/asin 正弦/反正弦函数 sec/asec 正割/反正割函数
sinh/asinh 双曲正弦/反双曲正弦函数 sech/asech 双曲正割/反双曲正割函数
cos/acos 余弦/反余弦函数 csc/acsc 余割/反余割函数
cosh/acosh 双曲余弦/反双曲余弦函数 csch/acsch 双曲余割/反双曲余割函数
tan/atan 正切/反正切函数 cot/acot 余切/反余切函数
tanh/atanh 双曲正切/反双曲正切函数 coth/acoth 双曲余切/反双曲余切函数
atan2 四个象限内反正切函数
附录6.2指数函数函数名功能描述函数名功能描述
exp 指数函数 log10 常用对数函数
log 自然对数函数 sqrt平方根函数
附录6.3复数函数函数名功能描述函数名功能描述
abs 绝对值函数 imag 求虚部函数
angle 角相位函数 real 求实部函数
conj 共轭复数函数
附录6.4数值处理函数名功能描述函数名功能描述
fix 沿零方向取整 round 舍入取整
floor 沿-∞方向取整 rem 求除法的余数
ceil 沿+∞方向取整 sign 符号函数
附录6.5其他特殊数学函数函数名功能描述函数名功能描述
airy airy函数 erfcx 比例互补误差函数
besselh bessel函数(hankel函数)erfinv 逆误差函数
bessili 改进的第一类bessel函数 expint 指数积分函数
besselk 改进的第二类bessel函数 gamma gamma函数
besselj 第一类bessel函数 gammainc 非完全gamma函数
bessely 第二类bessel函数 gammaln gamma对数函数
beta beta函数 gcd 最大公约数
betainc 非完全的beta函数 lcm 最小公倍数
betaln beta对数函数 log2 分割浮点数
elipj Jacobi椭圆函数 legendre legendre伴随函数
ellipke 完全椭圆积分 pow2 基2标量浮点数
erf 误差函数 rat 有理逼近
erfc 互补误差函数 rats 有理输出
----------------------------- A a abs 绝对值、模、字符的ASCII码值 acos 反余弦 acosh 反双曲余弦 acot 反余切 acoth 反双曲余切 acsc 反余割 acsch 反双曲余割
align 启动图形对象几何位置排列工具 all 所有元素非零为真 angle 相角
ans 表达式计算结果的缺省变量名 any 所有元素非全零为真 area 面域图
argnames 函数M文件宗量名 asec 反正割 asech 反双曲正割 asin 反正弦 asinh 反双曲正弦 assignin 向变量赋值 atan 反正切 atan2 四象限反正切 atanh 反双曲正切 autumn 红黄调秋色图阵 axes 创建轴对象的低层指令 axis 控制轴刻度和风格的高层指令 B b
bar 二维直方图 bar3 三维直方图 bar3h 三维水平直方图 barh 二维水平直方图
base2dec X进制转换为十进制 bin2dec 二进制转换为十进制 blanks 创建空格串 bone 蓝色调黑白色图阵 box 框状坐标轴
break while 或for 环中断指令 brighten 亮度控制
C c
capture(3版以前)捕获当前图形 cart2pol 直角坐标变为极或柱坐标 cart2sph 直角坐标变为球坐标 cat 串接成高维数组 caxis 色标尺刻度 cd 指定当前目录 cdedit 启动用户菜单、控件回调函数设计工具 cdf2rdf 复数特征值对角阵转为实数块对角阵 ceil 向正无穷取整 cell 创建元胞数组
cell2struct 元胞数组转换为构架数组 celldisp 显示元胞数组内容 cellplot 元胞数组内部结构图示
char 把数值、符号、内联类转换为字符对象 chi2cdf 分布累计概率函数 chi2inv 分布逆累计概率函数 chi2pdf 分布概率密度函数 chi2rnd 分布随机数发生器 chol Cholesky分解 clabel 等位线标识 cla 清除当前轴
class 获知对象类别或创建对象 clc 清除指令窗
clear 清除内存变量和函数 clf 清除图对象 clock 时钟
colorcube 三浓淡多彩交叉色图矩阵 colordef 设置色彩缺省值 colormap 色图 colspace 列空间的基 close 关闭指定窗口 colperm 列排序置换向量 comet 彗星状轨迹图 comet3 三维彗星轨迹图 compass 射线图 compose 求复合函数 cond(逆)条件数
condeig 计算特征值、特征向量同时给出条件数 condest 范-1条件数估计 conj 复数共轭 contour 等位线 contourf 填色等位线 contour3 三维等位线
contourslice 四维切片等位线图 conv 多项式乘、卷积 cool 青紫调冷色图 copper 古铜调色图 cos 余弦 cosh 双曲余弦 cot 余切 coth 双曲余切
cplxpair 复数共轭成对排列 csc 余割 csch 双曲余割 cumsum 元素累计和 cumtrapz 累计梯形积分 cylinder 创建圆柱
D d
dblquad 二重数值积分 deal 分配宗量
deblank 删去串尾部的空格符 dec2base 十进制转换为X进制 dec2bin 十进制转换为二进制 dec2hex 十进制转换为十六进制 deconv 多项式除、解卷 delaunay Delaunay 三角剖分 del2 离散Laplacian差分 demo Matlab演示 det 行列式
diag 矩阵对角元素提取、创建对角阵 diary Matlab指令窗文本内容记录 diff 数值差分、符号微分
digits 符号计算中设置符号数值的精度 dir 目录列表 disp 显示数组
display 显示对象内容的重载函数 dlinmod 离散系统的线性化模型
dmperm 矩阵Dulmage-Mendelsohn 分解 dos 执行DOS 指令并返回结果
double 把其他类型对象转换为双精度数值 drawnow 更新事件队列强迫Matlab刷新屏幕 dsolve 符号计算解微分方程
E e
echo M文件被执行指令的显示 edit 启动M文件编辑器 eig 求特征值和特征向量 eigs 求指定的几个特征值
end 控制流FOR等结构体的结尾元素下标 eps 浮点相对精度
error 显示出错信息并中断执行
errortrap 错误发生后程序是否继续执行的控制 erf 误差函数 erfc 误差补函数 erfcx 刻度误差补函数 erfinv 逆误差函数
errorbar 带误差限的曲线图 etreeplot 画消去树 eval 串演算指令 evalin 跨空间串演算指令 exist 检查变量或函数是否已定义 exit 退出Matlab环境 exp 指数函数
expand 符号计算中的展开操作 expint 指数积分函数 expm 常用矩阵指数函数 expm1 Pade法求矩阵指数 expm2 Taylor法求矩阵指数 expm3 特征值分解法求矩阵指数 eye 单位阵
ezcontour 画等位线的简捷指令 ezcontourf 画填色等位线的简捷指令 ezgraph3 画表面图的通用简捷指令 ezmesh 画网线图的简捷指令
ezmeshc 画带等位线的网线图的简捷指令 ezplot 画二维曲线的简捷指令 ezplot3 画三维曲线的简捷指令 ezpolar 画极坐标图的简捷指令 ezsurf 画表面图的简捷指令
ezsurfc 画带等位线的表面图的简捷指令
F f
factor 符号计算的因式分解 feather 羽毛图 feedback 反馈连接 feval 执行由串指定的函数 fft 离散Fourier变换 fft2 二维离散Fourier变换 fftn 高维离散Fourier变换 fftshift 直流分量对中的谱 fieldnames 构架域名 figure 创建图形窗 fill3 三维多边形填色图 find 寻找非零元素下标
findobj 寻找具有指定属性的对象图柄 findstr 寻找短串的起始字符下标 findsym 机器确定内存中的符号变量 finverse 符号计算中求反函数 fix 向零取整
flag 红白蓝黑交错色图阵 fliplr 矩阵的左右翻转 flipud 矩阵的上下翻转 flipdim 矩阵沿指定维翻转 floor 向负无穷取整 flops 浮点运算次数 flow Matlab提供的演示数据
fmin 求单变量非线性函数极小值点(旧版)fminbnd 求单变量非线性函数极小值点 fmins 单纯形法求多变量函数极小值点(旧版)fminunc 拟牛顿法求多变量函数极小值点 fminsearch 单纯形法求多变量函数极小值点 fnder 对样条函数求导 fnint 利用样条函数求积分
fnval 计算样条函数区间内任意一点的值 fnplt 绘制样条函数图形 fopen 打开外部文件 for 构成for环用 format 设置输出格式 fourier Fourier 变换 fplot 返函绘图指令 fprintf 设置显示格式 fread 从文件读二进制数据 fsolve 求多元函数的零点 full 把稀疏矩阵转换为非稀疏阵 funm 计算一般矩阵函数 funtool 函数计算器图形用户界面 fzero 求单变量非线性函数的零点
G g
gamma 函数
gammainc 不完全函数 gammaln 函数的对数 gca 获得当前轴句柄
gcbo 获得正执行“回调”的对象句柄 gcf 获得当前图对象句柄 gco 获得当前对象句柄 geomean 几何平均值 get 获知对象属性 getfield 获知构架数组的域 getframe 获取影片的帧画面 ginput 从图形窗获取数据 global 定义全局变量 gplot 依图论法则画图 gradient近似梯度 gray 黑白灰度 grid 画分格线
griddata 规则化数据和曲面拟合 gtext 由鼠标放置注释文字
guide 启动图形用户界面交互设计工具
H h
harmmean 调和平均值 help 在线帮助
helpwin 交互式在线帮助
helpdesk 打开超文本形式用户指南 hex2dec 十六进制转换为十进制 hex2num 十六进制转换为浮点数 hidden 透视和消隐开关 hilb Hilbert矩阵
hist 频数计算或频数直方图 histc 端点定位频数直方图 histfit 带正态拟合的频数直方图 hold 当前图上重画的切换开关 horner 分解成嵌套形式 hot 黑红黄白色图 hsv 饱和色图
I i
if-else-elseif 条件分支结构 ifft 离散Fourier反变换 ifft2 二维离散Fourier反变换 ifftn 高维离散Fourier反变换 ifftshift 直流分量对中的谱的反操作 ifourier Fourier反变换 i, j 缺省的“虚单元”变量 ilaplace Laplace反变换 imag 复数虚部 image 显示图象 imagesc 显示亮度图象 imfinfo 获取图形文件信息 imread 从文件读取图象 imwrite 把
imwrite 把图象写成文件 ind2sub 单下标转变为多下标 inf 无穷大
info MathWorks公司网点地址 inline 构造内联函数对象 inmem 列出内存中的函数名 input 提示用户输入 inputname 输入宗量名 int 符号积分
int2str 把整数数组转换为串数组 interp1 一维插值 interp2 二维插值 interp3 三维插值 interpn N维插值 interpft 利用FFT插值 intro Matlab自带的入门引导 inv 求矩阵逆
invhilb Hilbert矩阵的准确逆 ipermute 广义反转置 isa 检测是否给定类的对象 ischar 若是字符串则为真 isequal 若两数组相同则为真 isempty 若是空阵则为真 isfinite 若全部元素都有限则为真 isfield 若是构架域则为真 isglobal 若是全局变量则为真 ishandle 若是图形句柄则为真
ishold 若当前图形处于保留状态则为真 isieee 若计算机执行IEEE规则则为真 isinf 若是无穷数据则为真 isletter 若是英文字母则为真 islogical 若是逻辑数组则为真 ismember 检查是否属于指定集 isnan 若是非数则为真 isnumeric 若是数值数组则为真 isobject 若是对象则为真 isprime 若是质数则为真 isreal 若是实数则为真 isspace 若是空格则为真 issparse 若是稀疏矩阵则为真 isstruct 若是构架则为真
isstudent 若是Matlab学生版则为真 iztrans 符号计算Z反变换
J j , K k
jacobian 符号计算中求Jacobian 矩阵 jet 蓝头红尾饱和色
jordan 符号计算中获得 Jordan标准型 keyboard 键盘获得控制权
kron Kronecker乘法规则产生的数组
L l
laplace Laplace变换 lasterr 显示最新出错信息 lastwarn 显示最新警告信息
leastsq 解非线性最小二乘问题(旧版)legend 图形图例 lighting 照明模式 line 创建线对象 lines 采用plot 画线色
linmod 获连续系统的线性化模型 linmod2 获连续系统的线性化精良模型 linspace 线性等分向量 ln 矩阵自然对数
load 从MAT文件读取变量 log 自然对数 log10 常用对数 log2 底为2的对数 loglog 双对数刻度图形 logm 矩阵对数 logspace 对数分度向量 lookfor 按关键字搜索M文件 lower 转换为小写字母
lsqnonlin 解非线性最小二乘问题 lu LU分解
M m
mad平均绝对值偏差 magic 魔方阵
maple &nb, sp;运作 Maple格式指令 mat2str 把数值数组转换成输入形态串数组 material 材料反射模式 max 找向量中最大元素
mbuild 产生EXE文件编译环境的预设置指令 mcc 创建MEX或EXE文件的编译指令 mean 求向量元素的平均值 median 求中位数
menuedit 启动设计用户菜单的交互式编辑工具 mesh 网线图 meshz 垂帘网线图 meshgrid 产生“格点”矩阵
methods 获知对指定类定义的所有方法函数 mex 产生MEX文件编译环境的预设置指令 mfunlis 能被mfun计算的MAPLE经典函数列表 mhelp 引出 Maple的在线帮助 min 找向量中最小元素 mkdir 创建目录
mkpp 逐段多项式数据的明晰化 mod 模运算
more 指令窗中内容的分页显示 movie 放映影片动画
moviein 影片帧画面的内存预置
mtaylor 符号计算多变量Taylor级数展开
N n
ndims 求数组维数 NaN 非数(预定义)变量 nargchk 输入宗量数验证 nargin 函数输入宗量数 nargout 函数输出宗量数 ndgrid 产生高维格点矩阵 newplot 准备新的缺省图、轴 nextpow2 取最接近的较大2次幂 nnz 矩阵的非零元素总数 nonzeros 矩阵的非零元素 norm 矩阵或向量范数
normcdf 正态分布累计概率密度函数 normest 估计矩阵2范数
norminv 正态分布逆累计概率密度函数 normpdf 正态分布概率密度函数 normrnd 正态随机数发生器
notebook 启动Matlab和Word的集成环境 null 零空间
num2str 把非整数数组转换为串
numden 获取最小公分母和相应的分子表达式 nzmax 指定存放非零元素所需内存
O o
ode1 非Stiff 微分方程变步长解算器 ode15s Stiff 微分方程变步长解算器 ode23t 适度Stiff 微分方程解算器 ode23tb Stiff 微分方程解算器 ode45 非Stiff 微分方程变步长解算器 odefile ODE 文件模板
odeget 获知ODE 选项设置参数
odephas2 ODE 输出函数的二维相平面图 odephas3 ODE 输出函数的三维相空间图 odeplot ODE 输出函数的时间轨迹图 odeprint 在Matlab指令窗显示结果 odeset 创建或改写 ODE选项构架参数值 ones 全1数组
optimset 创建或改写优化泛函指令的选项参数值 orient 设定图形的排放方式 orth 值空间正交化
P p
pack 收集Matlab内存碎块扩大内存 pagedlg 调出图形排版对话框 patch 创建块对象
path 设置Matlab搜索路径的指令 pathtool 搜索路径管理器 pause 暂停
pcode 创建预解译P码文件 pcolor 伪彩图 peaks Matlab提供的典型三维曲面 permute 广义转置 pi(预定义变量)圆周率 pie 二维饼图 pie3 三维饼图 pink 粉红色图矩阵 pinv 伪逆 plot平面线图 plot3 三维线图
plotmatrix 矩阵的散点图 plotyy 双纵坐标图
poissinv 泊松分布逆累计概率分布函数 poissrnd 泊松分布随机数发生器 pol2cart 极或柱坐标变为直角坐标 polar 极坐标图
poly 矩阵的特征多项式、根集对应的多项式 poly2str 以习惯方式显示多项式
poly2sym 双精度多项式系数转变为向量符号多项式 polyder 多项式导数 polyfit 数据的多项式拟合 polyval 计算多项式的值 polyvalm 计算矩阵多项式 pow2 2的幂
ppval 计算分段多项式
pretty 以习惯方式显示符号表达式 print 打印图形或SIMULINK模型 printsys 以习惯方式显示有理分式 prism 光谱色图矩阵
procread 向MAPLE输送计算程序 profile 函数文件性能评估器 propedit 图形对象属性编辑器 pwd 显示当前工作目录
Q q
quad 低阶法计算数值积分
quad8 高阶法计算数值积分(QUADL)quit 推出Matlab 环境 quiver 二维方向箭头图 quiver3 三维方向箭头图
R r
rand 产生均匀分布随机数 randn 产生正态分布随机数 randperm 随机置换向量 range 样本极差 rank 矩阵的秩 rats 有理输出
rcond 矩阵倒条件数估计 real 复数的实部
reallog 在实数域内计算自然对数 realpow 在实数域内计算乘方 realsqrt 在实数域内计算平方根 realmax 最大正浮点数 realmin 最小正浮点数 rectangle 画“长方框” rem 求余数
repmat 铺放模块数组 reshape 改变数组维数、大小 residue 部分分式展开 return 返回
ribbon 把二维曲线画成三维彩带图 rmfield 删去构架的域 roots 求多项式的根 rose 数扇形图 rot90 矩阵旋转90度 rotate 指定的原点和方向旋转
rotate3d 启动三维图形视角的交互设置功能 round 向最近整数圆整 rref 简化矩阵为梯形形式
rsf2csf 实数块对角阵转为复数特征值对角阵 rsums Riemann和
S s
save 把内存变量保存为文件 scatter 散点图 scatter3 三维散点图 sec 正割 sech 双曲正割
semilogx X轴对数刻度坐标图 semilogy Y轴对数刻度坐标图 series 串联连接 set 设置图形对象属性 setfield 设置构架数组的域 setstr 将ASCII码转换为字符的旧版指令 sign 根据符号取值函数
signum 符号计算中的符号取值函数 sim 运行SIMULINK模型
simget 获取SIMULINK模型设置的仿真参数 simple 寻找最短形式的符号解 simplify 符号计算中进行简化操作
simset 对SIMULINK模型的仿真参数进行设置 simulink 启动SIMULINK模块库浏览器 sin 正弦 sinh 双曲正弦 size 矩阵的大小 slice 立体切片图
solve 求代数方程的符号解 spalloc 为非零元素配置内存 sparse 创建稀疏矩阵
spconvert 把外部数据转换为稀疏矩阵 spdiags 稀疏对角阵 spfun 求非零元素的函数值 sph2cart 球坐标变为直角坐标 sphere 产生球面
spinmap 色图彩色的周期变化 spline 样条插值
spones 用1置换非零元素 sprandsym 稀疏随机对称阵 sprank 结构秩 spring 紫黄调春色图 sprintf 把格式数据写成串 spy 画稀疏结构图 sqrt平方根 sqrtm 方根矩阵
squeeze 删去大小为1的“孤维” sscanf 按指定格式读串 stairs 阶梯图 std 标准差 stem 二维杆图 step 阶跃响应指令
str2double 串转换为双精度值 str2mat 创建多行串数组 str2num 串转换为数 strcat 接成长串 strcmp 串比较 strjust 串对齐 strmatch 搜索指定串 strncmp 串中前若干字符比较 strrep 串替换
strtok 寻找第一间隔符前的内容 struct 创建构架数组
struct2cell 把构架转换为元胞数组 strvcat 创建多行串数组 sub2ind 多下标转换为单下标 subexpr 通过子表达式重写符号对象 subplot 创建子图
subs 符号计算中的符号变量置换 subspace 两子空间夹角 sum 元素和
summer 绿黄调夏色图 superiorto 设定优先级 surf 三维着色表面图 surface 创建面对象 surfc 带等位线的表面图 surfl 带光照的三维表面图 surfnorm 空间表面的法线 svd 奇异值分解
svds 求指定的若干奇异值 switch-case-otherwise 多分支结构
sym2poly 符号多项式转变为双精度多项式系数向量 symmmd 对称最小度排序 symrcm 反向Cuthill-McKee排序 syms 创建多个符号对象 T t
tan 正切 tanh 双曲正切
taylortool 进行Taylor逼近分析的交互界面 text 文字注释 tf 创建传递函数对象 tic 启动计时器 title 图名 toc 关闭计时器 trapz 梯形法数值积分 treelayout 展开树、林 treeplot 画树图 tril 下三角阵 trim 求系统平衡点 trimesh 不规则格点网线图
trisurf 不规则格点表面图 triu 上三角阵 try-catch 控制流中的Try-catch结构 type 显示M文件 U u uicontextmenu 创建现场菜单 uicontrol 创建用户控件 uimenu 创建用户菜单
unmkpp 逐段多项式数据的反明晰化 unwrap 自然态相角 upper 转换为大写字母
V v
var 方差
varargin 变长度输入宗量 varargout 变长度输出宗量
vectorize 使串表达式或内联函数适于数组运算 ver 版本信息的获取 view 三维图形的视角控制 voronoi Voronoi多边形 vpa 任意精度(符号类)数值
W w
warning 显示警告信息 what 列出当前目录上的文件
whatsnew 显示Matlab中 Readme文件的内容 which 确定函数、文件的位置 while 控制流中的While环结构 white 全白色图矩阵 whitebg 指定轴的背景色 who 列出内存中的变量名 whos 列出内存中变量的详细信息 winter 蓝绿调冬色图 workspace 启动内存浏览器
X x , Y y , Z z
xlabel X轴名 xor 或非逻辑
yesinput 智能输入指令 ylabel Y轴名 zeros 全零数组 zlabel Z轴名
zoom 图形的变焦放大和缩小 ztrans 符号计算Z变换
第四篇:matlab制图函数总结
Subplot(a,b,c)图像位置函数,a表示分成的行数,b表示当前行的列数,c为位置序号。
Plot(x1,y1,’s1’,x2,y2,’s2’,……)二维绘图函数,绘制一般曲线,参数x表示x轴量,y表示y轴量,s为曲线颜色及形状参数。
Axis([x1,x2,y1,y2])二维绘图函数,参数x1和x2为x轴初始及末尾值,y1和y2为y轴值。
Stair(x1,y1)二维绘图函数,绘制台阶型曲线,参数x表示x轴量,y表示y轴量。
Hold on 保持之前的图形,同时显示之后的图形。
[x,y,z]=cylinder(f(x),s)三维制图函数,绘制柱状立体图,f(x)为边界曲线函数,s为边界曲线条数。
[x,y,z]=sphere(s)三围制图函数,绘制球状立体图,s为球体各个侧面的图块数,默认为30。
Figure('toolbar','none',...%是否显示工具栏,否
'name','制图',...%对话框名称,制图
'NumberTitle','off',...%是否显示对话框编号,否
'color','w',...%背景颜色,白色
'Resize','on',...%是否可调对话框大小,是
'Position',[300,100,700,600]);
%默认对话框大小及位置 对话框属性编辑函数
第五篇:MATLAB程序总结
%******************************读取数据************************************* MATLAB读取数据
Xlsread(‘lujing’,’mingcheng’)%******************************输出数据************************************* MATLAB输出Excel数据
Xlswrite(‘eular.xla’,[A],’sheet 1’)%******************************三维曲面************************************* 三维曲面
x1=[0.05 0.05 0.05 0.1 0.05 0.1 0.2 0 0 0 0 0 0 0 0];x2=[0.1 0.25 0.12 0.12 0 0 0 0 0.05 0.08 0.1 0.12 0.15 0.2 0.25];y1=[153 130 146 133 160 154 140 133 165 169 171 186 175 156 152];y2=[7.22 5.57 6.66 6 6.6 8.21 5.1 4.66 7.58 8.26 8.99 8.73 7.71 7 6.49];y3=[112 100 131 88 80 72 72 80 116 83 80 80 64 68 74];scatter3(x1,x2,y2)figure [X,Y,Z]=griddata(x1,x2,y2,linspace(min(x1),max(x1))',linspace(min(x2),max(x2)),'v4');%插值 pcolor(X,Y,Z);shading interp%伪彩色图 figure,contourf(X,Y,Z)%等高线图 figure,surf(X,Y,Z);%三维曲面
%******************************各种距离************************************* pdist(X)— 样本X中各n维向量的欧氏距离 pdist(X, 'cityblock')— 各n维向量的绝对距离 pdist(X, ‘minkowski',r)— 闵可夫斯基距离 pdist(X, 'mahal')— 各n维向量的马氏距离 mean(A)— A中各列向量的均值
Var(A)— A中各列向量的方差
Std(A)— A中各列向量的标准差
Cov(A)— A中各列向量的协方差矩阵 Corrcoef(A)— A中各列向量的相关矩阵 命令:[x,fval] = fminbnd(F,a,b)%******************************泰勒公式************************************* 计算F在a,b之间取极小值时的x与y(即fval).格式:taylor(F,a)功能:F在x=a处的泰勒级数前5项 格式:taylor(F,v)功能:F对变量v的泰勒展式前5项
格式:taylor(F,v,n)功能:求F的n 阶泰勒展式,且
(n缺省时默认为 5)
%******************************一维插值************************************* x=[-2*pi:0.5*pi:2*pi];y=cos(x);xi=-2*pi:0.3*pi:2*pi;y_nearest=interp1(x,y, xi,‘nearest’)%最近邻插值 y_linear= interp1(x,y,xi)%双线性差值
y_spline= interp1(x,y,xi, ‘spline’)%三次样条插值 y_cubic= interp1(x,y,xi, ‘cubic’)%双三次插值
plot(x,y,'o',xi,y_nearest,'-',xi,y_linear, 'r* ', xi,y_spline,'k:',xi,y_cubic,'k-');legend('original data','nearest','linear','spline','cubic')%******************************曲线拟合************************************* x1=[2:16];y1=[6.42,8.2,9.58,9.5,9.7,10,9.93,9.99,10.49,10.59,10.6,10.8,10.6,10.9,10.76];b01=[0.1435,0.084];
%初始参数值
fun1=inline('x./(b(1)+b(2)*x)','b','x');% 定义函数 [b1,r1,j1]=nlinfit(x1,y1,fun1,b01);y=x1./(0.1152+0.0845*x1);%根据b1写出具体函数
plot(x1,y1,'*',x1,y,'-or');R2=1-sum((y-y1).^2)/sum((y1-mean(y1)).^2)%*****************************残差与置信区间******************************** X=[ones(30,1), x1‘,x2’,x3‘];
%输入自变量(注意1与自变量组成的矩阵)[b,bint,r,rint,s]=regress(y',X);%多元线性回归 b,bint,s rcoplot(r,rint)%残差及其置信区间作图
a=[1,3:9,11:30];X1=X(a,:);y1=y(a);%剔除异常点
[b1,bint1,r1,rint1,s1]=regress(y1',X1);b1,bint1,s1 %再次线性回归 rcoplot(r1,rint1)
%------------%
Copula理论及其在matlab中的实现程序应用实例 %------------
%******************************读取数据************************************* % 从文件hushi.xls中读取数据 hushi = xlsread('hushi.xls');% 提取矩阵hushi的第5列数据,即沪市的日收益率数据 X = hushi(:,5);% 从文件shenshi.xls中读取数据 shenshi = xlsread('shenshi.xls');% 提取矩阵shenshi的第5列数据,即深市的日收益率数据 Y = shenshi(:,5);
%****************************绘制频率直方图*********************************
% 调用ecdf函数和ecdfhist函数绘制沪、深两市日收益率的频率直方图 [fx, xc] = ecdf(X);figure;ecdfhist(fx, xc, 30);xlabel('沪市日收益率');% 为X轴加标签 ylabel('f(x)');% 为Y轴加标签 [fy, yc] = ecdf(Y);figure;ecdfhist(fy, yc, 30);xlabel('深市日收益率');% 为X轴加标签 ylabel('f(y)');% 为Y轴加标签
%****************************计算偏度和峰度********************************* % 计算X和Y的偏度 xs = skewness(X)ys = skewness(Y)
% 计算X和Y的峰度 kx = kurtosis(X)ky = kurtosis(Y)
%******************************正态性检验*********************************** % 分别调用jbtest、kstest和lillietest函数对X进行正态性检验 [h,p] = jbtest(X)% Jarque-Bera检验
[h,p] = kstest(X,[X,normcdf(X,mean(X),std(X))])% Kolmogorov-Smirnov检验 [h, p] = lillietest(X)% Lilliefors检验
% 分别调用jbtest、kstest和lillietest函数对Y进行正态性检验 [h,p] = jbtest(Y)% Jarque-Bera检验
[h,p] = kstest(Y,[Y,normcdf(Y,mean(Y),std(Y))])% Kolmogorov-Smirnov检验 [h, p] = lillietest(Y)% Lilliefors检验
%****************************求经验分布函数值******************************* % 调用ecdf函数求X和Y的经验分布函数 [fx, Xsort] = ecdf(X);[fy, Ysort] = ecdf(Y);% 调用spline函数,利用样条插值法求原始样本点处的经验分布函数值 U1 = spline(Xsort(2:end),fx(2:end),X);V1 = spline(Ysort(2:end),fy(2:end),Y);
% 调用ecdf函数求X和Y的经验分布函数 [fx, Xsort] = ecdf(X);[fy, Ysort] = ecdf(Y);% 提取fx和fy的第2个至最后一个元素,即排序后样本点处的经验分布函数值 fx = fx(2:end);fy = fy(2:end);
% 通过排序和反排序恢复原始样本点处的经验分布函数值U1和V1 [Xsort,id] = sort(X);[idsort,id] = sort(id);U1 = fx(id);[Ysort,id] = sort(Y);[idsort,id] = sort(id);V1 = fy(id);
%*******************************核分布估计********************************** % 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值 U2 = ksdensity(X,X,'function','cdf');V2 = ksdensity(Y,Y,'function','cdf');
% **********************绘制经验分布函数图和核分布估计图********************** [Xsort,id] = sort(X);% 为了作图的需要,对X进行排序 figure;% 新建一个图形窗口
plot(Xsort,U1(id),'c','LineWidth',5);% 绘制沪市日收益率的经验分布函数图 hold on plot(Xsort,U2(id),'k-.','LineWidth',2);% 绘制沪市日收益率的核分布估计图 legend('经验分布函数','核分布估计', 'Location','NorthWest');% 加标注框 xlabel('沪市日收益率');% 为X轴加标签 ylabel('F(x)');% 为Y轴加标签
[Ysort,id] = sort(Y);% 为了作图的需要,对Y进行排序 figure;% 新建一个图形窗口
plot(Ysort,V1(id),'c','LineWidth',5);% 绘制深市日收益率的经验分布函数图 hold on plot(Ysort,V2(id),'k-.','LineWidth',2);% 绘制深市日收益率的核分布估计图 legend('经验分布函数','核分布估计', 'Location','NorthWest');% 加标注框 xlabel('深市日收益率');% 为X轴加标签 ylabel('F(x)');% 为Y轴加标签
%****************************绘制二元频数直方图***************************** % 调用ksdensity函数分别计算原始样本X和Y处的核分布估计值 U = ksdensity(X,X,'function','cdf');V = ksdensity(Y,Y,'function','cdf');figure;% 新建一个图形窗口
% 绘制边缘分布的二元频数直方图,hist3([U(:)V(:)],[30,30])xlabel('U(沪市)');% 为X轴加标签 ylabel('V(深市)');% 为Y轴加标签 zlabel('频数');% 为z轴加标签
%****************************绘制二元频率直方图***************************** figure;% 新建一个图形窗口
% 绘制边缘分布的二元频数直方图,hist3([U(:)V(:)],[30,30])h = get(gca, 'Children');% 获取频数直方图的句柄值 cuv = get(h, 'ZData');% 获取频数直方图的Z轴坐标
set(h,'ZData',cuv*30*30/length(X));% 对频数直方图的Z轴坐标作变换 xlabel('U(沪市)');% 为X轴加标签 ylabel('V(深市)');% 为Y轴加标签 zlabel('c(u,v)');% 为z轴加标签
%***********************求Copula中参数的估计值****************************** % 调用copulafit函数估计二元正态Copula中的线性相关参数 rho_norm = copulafit('Gaussian',[U(:), V(:)])% 调用copulafit函数估计二元t-Copula中的线性相关参数和自由度 [rho_t,nuhat,nuci] = copulafit('t',[U(:), V(:)])
%********************绘制Copula的密度函数和分布函数图************************ [Udata,Vdata] = meshgrid(linspace(0,1,31));% 为绘图需要,产生新的网格数据 % 调用copulapdf函数计算网格点上的二元正态Copula密度函数值 Cpdf_norm = copulapdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);% 调用copulacdf函数计算网格点上的二元正态Copula分布函数值 Ccdf_norm = copulacdf('Gaussian',[Udata(:), Vdata(:)],rho_norm);% 调用copulapdf函数计算网格点上的二元t-Copula密度函数值 Cpdf_t = copulapdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);% 调用copulacdf函数计算网格点上的二元t-Copula分布函数值 Ccdf_t = copulacdf('t',[Udata(:), Vdata(:)],rho_t,nuhat);% 绘制二元正态Copula的密度函数和分布函数图 figure;% 新建图形窗口
surf(Udata,Vdata,reshape(Cpdf_norm,size(Udata)));% 绘制二元正态Copula密度函数图 xlabel('U');% 为X轴加标签 ylabel('V');% 为Y轴加标签 zlabel('c(u,v)');% 为z轴加标签 figure;% 新建图形窗口
surf(Udata,Vdata,reshape(Ccdf_norm,size(Udata)));% 绘制二元正态Copula分布函数图 xlabel('U');% 为X轴加标签 ylabel('V');% 为Y轴加标签 zlabel('C(u,v)');% 为z轴加标签
% 绘制二元t-Copula的密度函数和分布函数图 figure;% 新建图形窗口
surf(Udata,Vdata,reshape(Cpdf_t,size(Udata)));% 绘制二元t-Copula密度函数图 xlabel('U');% 为X轴加标签 ylabel('V');% 为Y轴加标签 zlabel('c(u,v)');% 为z轴加标签 figure;% 新建图形窗口
surf(Udata,Vdata,reshape(Ccdf_t,size(Udata)));% 绘制二元t-Copula分布函数图 xlabel('U');% 为X轴加标签 ylabel('V');% 为Y轴加标签 zlabel('C(u,v)');% 为z轴加标签
%**************求Kendall秩相关系数和Spearman秩相关系数*********************** % 调用copulastat函数求二元正态Copula对应的Kendall秩相关系数 Kendall_norm = copulastat('Gaussian',rho_norm)% 调用copulastat函数求二元正态Copula对应的Spearman秩相关系数 Spearman_norm = copulastat('Gaussian',rho_norm,'type','Spearman')% 调用copulastat函数求二元t-Copula对应的Kendall秩相关系数 Kendall_t = copulastat('t',rho_t)% 调用copulastat函数求二元t-Copula对应的Spearman秩相关系数 Spearman_t = copulastat('t',rho_t,'type','Spearman')
% 直接根据沪、深两市日收益率的原始观测数据,调用corr函数求Kendall秩相关系数 Kendall = corr([X,Y],'type','Kendall')% 直接根据沪、深两市日收益率的原始观测数据,调用corr函数求Spearman秩相关系数 Spearman = corr([X,Y],'type','Spearman')
%******************************模型评价************************************* % 调用ecdf函数求X和Y的经验分布函数 [fx, Xsort] = ecdf(X);[fy, Ysort] = ecdf(Y);% 调用spline函数,利用样条插值法求原始样本点处的经验分布函数值 U = spline(Xsort(2:end),fx(2:end),X);V = spline(Ysort(2:end),fy(2:end),Y);% 定义经验Copula函数C(u,v)C = @(u,v)mean((U <= u).*(V <= v));% 为作图的需要,产生新的网格数据 [Udata,Vdata] = meshgrid(linspace(0,1,31));% 通过循环计算经验Copula函数在新产生的网格点处的函数值 for i=1:numel(Udata)
CopulaEmpirical(i)= C(Udata(i),Vdata(i));end
figure;% 新建图形窗口
% 绘制经验Copula分布函数图像
surf(Udata,Vdata,reshape(CopulaEmpirical,size(Udata)))xlabel('U');% 为X轴加标签 ylabel('V');% 为Y轴加标签
zlabel('Empirical Copula C(u,v)');% 为z轴加标签
% 通过循环计算经验Copula函数在原始样本点处的函数值 CUV = zeros(size(U(:)));for i=1:numel(U)
CUV(i)= C(U(i),V(i));end
% 计算线性相关参数为0.9264的二元正态Copula函数在原始样本点处的函数值 rho_norm = 0.9264;Cgau = copulacdf('Gaussian',[U(:), V(:)],rho_norm);% 计算线性相关参数为0.9325,自由度为4的二元t-Copula函数在原始样本点处的函数值 rho_t = 0.9325;k = 4.0089;Ct = copulacdf('t',[U(:), V(:)],rho_t,k);% 计算平方欧氏距离
dgau2 =(CUV-Cgau)'*(CUV-Cgau)dt2 =(CUV-Ct)'*(CUV-Ct)
%******************************灰色预测*************************************
A=[ 96 144 194 276 383 466 554 652 747 832 972
169 235 400 459 565 695 805 881 1011 1139
151 238 335 425 541 641 739 866 975 1087 1238
164 263 376 531 600 711 913 1038 1173 1296 1497
182 318 445 576 708 856 1000 1145 1292 1435 1667
216 361 504 642 818 979 1142 1305 1479 1644 1920] m=size(A,2);x0=mean(A,2);x1=cumsum(x0);alpha=0.4;n=length(x0);z1=alpha*x1(2:n)+(1-alpha)*x1(1:n-1);Y=x0(2:n);B=[-z1,ones(n-1,1)];ab=BY
x_A=(x0(1)-ab(2)/ab(1))*(exp(-ab(1)*n)-exp(-ab(1)*(n-1)))
图标参数************************************* y 黄-实线.点< 小于号 m 紫: 点线o 圆s 正方形 c 青-.点划线x 叉号d 菱形 r 红--虚线+ 加号h 六角星 g 绿* 星号p 五角星 b 蓝v 向下三角形 w 白^ 向上三角形 k 黑> 大于号
%******************************图形修饰*************************************
函数 含义
grid on(/off)给当前图形标记添加(取消)网络 xlable(‘string’)标记横坐标 ylabel(‘string’)标记纵坐标 title(‘string’)给图形添加标题
text(x,y,’string’)在图形的任意位置增加说明性文本信息 gtext(‘string’)利用鼠标添加说明性文本信息
axis([xmin xmax ymin ymax])设置坐标轴的最小最大值
例5.1.2 给例5.1.1 的图形中加入网络和标记。(见图5.1.3 和5.1.4)>> x=0:pi/10:2*pi;>> y1=sin(x);>> y2=cos(x);>> plot(x,y1,x,y2)>> grid on >> xlabel('independent variable X')>> ylabel('Dependent Variable Y1 & Y2')>> title('Sine and Cosine Curve')>> text(1.5,0.3,'cos(x)')>> gtext('sin(x)')>> axis([0 2*pi-0.9 0.9])
图5.1.3 使用了图形修饰的plot 函数绘制的正弦曲线
例5.2.1 绘制方程x=t y=sin(t)z=cos(t)
在t=[0,2*pi]上的空间方程。(见图5.2.1)
>> clf >> x=0:pi/10:2*pi;>> y1=sin(x);>> y2=cos(x);>> plot3(y1,y2,x,'m:p')>> grid on >> xlabel('Dependent Variable Y1')>> ylabel('Dependent Variable Y2')>> zlabel('Independent Variable X')>> title('Sine and Cosine Curve')
图5.2.1 函数plot 绘制的三维曲线图
%******************************三维图形*************************************
例5.2.3 绘制方程
sin((x^2+y^2)^(1/2))z =---------------------(x^2+y^2)^(1/2)
在x∈[-7.5,7.5];y∈[-7.5,7.5] 的图形。
>> x=-7.5:0.5:7.5;y=x;>> [X,Y]=meshgrid(x,y);>> R=sqrt(X.^2+Y.^2)+eps;>> Z=sin(R)./R;>> surf(X,Y,Z)>> xlabel('X 轴方向')>> ylabel('Y 轴方向')>> zlabel('Z 轴方向')(见图5.2.4)
图5.2.4
例5.2.4 绘制由方程形成的立体图。(见图5.2.5)z=x*exp(-(x^2+y^2))
>> clear >> x=-2:0.1:2;y=x;>> [X,Y]=meshgrid(x,y);>> Z=X.*exp(-X.^2-Y.^2);>> surf(X,Y,Z)
图5.2.5
%******************************三维多角度观察*********************************
MTALAB 允许用户设置观察点,其指令是: view(azimuth,elevation)其中方位角azimuth 是观察点和坐标原点连线在x-y平面的投影和y 轴负方向的夹角,仰
角elevation 是观察点与坐标原点的连线和x-y平面的夹角。对于这两个角度,三维图形的
默认值分别是-37.5 和30,二维图形的默认值是0 和90。
例5.2.5 从不同的角度观察高斯矩阵的曲面。
>> z=peaks(40);>> subplot(2,2,1);>> mesh(z);>> subplot(2,2,2);>> mesh(z);>> view(-37.5,-30);>> subplot(2,2,3);>> mesh(z);>> view(180,0);>> subplot(2,2,4);>> mesh(z);>> view(0,90);
图5.2.6 对应不同观察点的三维曲面图
%******************************其他图形*************************************
除了plot 绘图函数以外,在有些场合对绘制的曲线会有一些特殊要求,这就要其他函
数来实现,常用的几种函数如下(见表5.3.1)
表5.3.1 其他图形函数表 函数 含义
loglog
使用对数坐标系绘图
semilogx
横坐标为对数坐标轴,纵坐标为线性坐标轴 semilogy
横坐标为线性坐标轴,纵坐标为对数坐标轴 polar
绘制极坐标图 fill
绘制实心图 bar
绘制直方图 pie
绘制饼图 area
绘制面积图 quiver 绘制向量场图 stairs 绘制阶梯图 sterm 绘制火柴杆图
例5.3.1 >> x=0:pi/10:2*pi;>> y1=sin(x);>> subplot(2,2,1);>> plot(x,y1);>> subplot(2,2,2);>> bar(x,y1);>> subplot(2,2,3);>> fill(x,y1,'g');>> subplot(2,2,4);>> stairs(x,y1,'k');
%******************************直方图*************************************
函数bar(x)可以绘制直方图,这对统计或者数据采集非常直观实用。它共有四种形式: bar,bar3,barh 和bar3h,其中bar 和bar3 分别用来绘制二维和三维竖直方图,barh 和b ar3h 分别用来绘制二维和三维水平直方图,调用格式是:
bar(x,y)其中x 必须单调递增或递减,y 为n m× 矩阵,可视化结果为m 组,每
组n 个垂直柱,也就是把y 的行画在一起,同一列的数据用相同的颜色表示; bar(x,y,width)(或bar(y,width))指定每个直方条的宽度,如width>1,则直方条会重
叠,默认值为width=0.8;
bar(…,’grouped’)使同一组直方条紧紧靠在一起; bar(…,’stack’)把同一组数据描述在一个直方条上。
例5.3.2
>> y=[5 3 2 9;4 7 2 7;1 5 7 3];>> subplot(2,2,1),bar(y)>> x=[5 9 11];>> subplot(2,2,2),bar3(x,y)>> subplot(2,2,3),bar(x,y,'grouped')>> subplot(2,2,4),bar(rand(2,3),.75,'stack')
%******************************面积图*************************************
5.3.2 面积图
函数area 用来绘制面积图,面积图在plot 的基础上填充x 轴和曲线之间的面积,该图
用于查看某个数在该列所有数的总和中所占的比例。
例5.3.3
>> x=-3:3;>> y=[3 2 5;6 1 8;7 4 9;6 3 7;8 2 9;4 2 9;3 1 7];>> area(x,y)
%******************************饼图*************************************
5.3.3 饼图
函数pie 用来绘制饼图,它可以形象地表示出向量中各元素所占比例。其调用格式是:
pie(x)x 中的元素通过x/sum(x)进行归一化,以确定饼图中的份额; pie(x,explode)向量explode 和x 元素数相同,用来指出需要分开的饼片,explode 中
不为零的部分会被分开。
例5.3.4 设某班的某课程的考试成绩如下:90 分以上有32 人,81 至90 有58 人,71 至80 分有27 人,60 至70 分为21 人,60 分以下有16 人,画出饼图。(见图5.3.4)
>> x=[32 58 27 21 16];>> explode0=[1 0 0 0 0];>> subplot(1,2,1)>> pie(x,explode0)>> explode1=[0 0 0 0 1];>> subplot(1,2,2)>> pie(x,explode1)
5.3.4 不同坐标系中的绘图
Semilogx,semilogy,loglo,polar(theta,rho)的使用方法和plot 完全类似,不同的只是绘
制到不同的图形坐标上。函数semilogx 绘制x 轴为对数标度的图形,在半对数坐标系中绘图; 函数semilogy 绘制y 轴为对数标度的图形;函数loglog 绘制两个轴都为对数间隔的图形;
函数polar(theta,rho)绘制极坐标图形,其中theta 为相角,rho 为其对应的半径。
例5.3.5 绘制ρ=acos(3θ),a=2 的图形。(见图5.3.5)
>> theta=-pi:pi/80:pi;>> polar(theta,2*cos(3*theta))
图5.3.5 极坐标图
5.4 符号表达式绘图
MATLAB 软件提供了将表达式进行图形显示的功能。完成此功能需调用fplot 函数和
ezplot 函数。
函数fplot 用来绘制数学函数,其调用格式为: fplot(fun,lims)其中fun 就是所要绘制的函数,可以是定义函数的M 文件名,也可以是以x 为变量的可计
算字符串。例如’diric(x,10)’或’[sin(x),cos(x)]’,对于向量x 的每个元素,函数 fun(x)必须返回一个行向量。如果fun 返回[f1(x),f2(x),f3(x)],输入[x1;x2],就会返回矩阵
f1(x1)f2(x1)f3(x1)f1(x2)f2(x2)f3(x2)lims=[XMIN XMAX YMIN YMAX]限定了x,y 轴上的绘图空间。
例5.4.1 >> subplot(2,2,1),fplot('humps',[0 1])>> subplot(2,2,2),fplot('abs(exp(-j*x*(0:9))*ones(10,1))',[0 2*pi])>> subplot(2,2,3),fplot('[tan(x),sin(x),cos(x)]',2*pi*[-1 1-1 1])>> subplot(2,2,4),fplot('sin(1./x)',[0.01 0.1],1e-3)
图5.4.1 fplot 函数绘制表达式图形
ezplot 函数是简捷绘图指令之一,它无需数据准备,直接画出函数图形,基本调用格式
为ezplot(f)其中f 是字符串或代表数学函数的符号表达式,只有一个符号变量,可以是x,缺省情况下
x 轴的绘图区域为[-π, π ],但我们可以用ezplot(f,xmin,xmax)或ezplot(f,[xmin,xmax])来指定x 的范围。
例5.4.2
>> y='x^2';>> subplot(1,2,1)>> ezplot(y)>> subplot(1,2,2)>> y='sin(x)';>> ezplot(y,[-pi,pi])
图5.4.2 ezplot 函数绘制表达式图形
5.5 plot 函数
MATLAB 对数据是按列存储和计算的,运用plot(x)时,当x 为一个向量时,以其元
素为纵坐标,其序号为横坐标值绘制曲线。当x 为实矩阵时,则以其序号为横坐标,按列
绘制每列元素相对于序号的曲线,当x 为n m× 矩阵时,就有n 条曲线。
如果x,y 是同维向量,plot(x,y)指令以x 元素为横坐标值,y 元素为纵坐标值绘制曲线。
如x 是向量,y 是有一维与x 元素数量相等的矩阵,则以x 为共同横坐标,按列绘制y 每
列元素值,曲线数为y 的另一维的元素数。如果x,y 是同维矩阵,则以x,y 对应列元素为、纵坐标分别绘制曲线,数目等于矩阵的列数。
例5.5.1
>> x=[3 5 10 8];>> subplot(2,2,1)>> plot(x)>> x=[3 5 10 8;7 2 9 4;2 7 2 7]';>> subplot(2,2,2)>> plot(x)>> x=[3 5 6 8];>> y=[1 5 10 4];>> subplot(2,2,3)>> plot(x,y)>> x=[1 3 5 7;2 4 6 8]';>> y=[6 2 5 10;3 5 2 6]';>> subplot(2,2,4)>> plot(x,y,'k:*')