当前位置:文档之家› 小波分解矩阵 Matlab

小波分解矩阵 Matlab

小波分解矩阵 Matlab
小波分解矩阵 Matlab

%----------------------------------------------------------%

小波图像分解Matlab 程序- V2.0版

https://www.doczj.com/doc/d53364852.html,/chenyusiyuan/archive/2008/06/05/2513865.aspx

小波图像重构Matlab 程序- V2.0版

https://www.doczj.com/doc/d53364852.html,/chenyusiyuan/archive/2008/06/05/2514119.aspx

%----------------------------------------------------------%

%----------------------------------------------------------%

小波分解矩阵Matlab 程序- V3.0版

%----------------------------------------------------------%

function [coef,scf]=mywavedec2(x,N,wname)

%----------------------------------------------------------%

% 函数MYWAVEDEC2() 对输入矩阵x 进行dim 层分解,得到相应的分解系数矩阵y % 输入参数:x ——输入矩阵

% N ——分解级数

% wname ——分解所用的小波函数

% 输出参数:scf ——存储各级分解系数矩阵的大小以及原始(图像)矩阵的大小

% coef ——分解系数矩阵,其结构如下:

% coef = {cA_N;cV_N;cH_N;cD_N;cV_N-1;cH_N-1;cD_N-1;……;cV_1;cH_1;cD_1}

%

% Copyright by Zou Yuhua ( chenyusiyuan )

% Version: 3.0, Date: 2008-07-08

%----------------------------------------------------------%

% 求出小波函数的滤波器组系数向量

[Lo_D,Hi_D] = wfilters(wname,'d');

lf=length(Lo_D);

% 画出原始图像

imshow(x);title('Original Image');

% 标明图像大小

[r,c]=size(x);

xlabel(['Size : ',num2str(r),'*',num2str(c)]);

% 将矩阵x的数据格式转换为适合数值处理的double格式

xd=double(x);

[rx,cx]=size(x);

[o2sa,f1sa,rsx]=sizcoef([rx,cx],lf,N); a=[o2sa,f1sa,rsx]

coef=[];

scf=[rx,cx];

for i=1:N

[cA,cV,cH,cD]=mydwt2(xd,Lo_D,Hi_D); % 第i 级小波分解

xd=cA; % 将第i 级分解得到的低频系数矩阵作为第i+1 级分解的源矩阵

outmp={cV;cH;cD}; % 将第i 级分解得到的高频系数矩阵cV,cH,cD存入细胞矩阵outmp

scf=[size(cV);scf]; % 将各级分解系数矩阵的大小存入矩阵scf

coef=[outmp;coef]; % 将细胞矩阵outmp 存入输出矩阵coef,coef将由空矩阵变为细胞矩阵

end

% 迭代结束后,矩阵coef 中保存的是各级分解中的高频系数矩阵

% 故需将迭代后得到的矩阵cA,即第dim 级低频矩阵存入矩阵coef

coef=[cA;coef];

scf=[size(cA);scf];

plotcoef(N,wname,coef);

plotcoef2one(N,wname,coef,rsx);

function plotcoef2one(N,wname,coef,rsx)

% 画出小波分解系数的塔式结构图

rsx=rsx(end:-1:1,:);

tmpcoef=[];

tA=wkeep(coef{1},rsx(1,:),'c');

tmpcoef=tA;

tA=uint8(tA); tA(end,:)=255; tA(:,end)=255;

for j=1:N

tV=wkeep(coef{(j-1)*3+2},rsx(j,:),'c');

tH=wkeep(coef{(j-1)*3+3},rsx(j,:),'c');

tD=wkeep(coef{(j-1)*3+4},rsx(j,:),'c');

tV=uint8(tV); tH=uint8(tH); tD=uint8(tD);

if j

tV(end,:)=255; tV(:,end)=255;

tH(end,:)=255; tH(:,end)=255;

tD(end,:)=255; tD(:,end)=255;

else

tV(end,:)=255; tH(:,end)=255;

end

tmpcoef=[tA,tV;tH,tD];stc=size(tmpcoef);

if stc>=rsx(j+1,:)

tA=tmpcoef(1:rsx(j+1,1),1:rsx(j+1,2));

else

tmpcoef=tmpcoef([1:end-1,end-2:end-1],[1:end-1,end-2:end-1]);

tA=tmpcoef(1:rsx(j+1,1),1:rsx(j+1,2));

end

tA=uint8(tA); tA(end,:)=255; tA(:,end)=255;

end

figure;

sc=tA;

[rx,cx]=size(sc);

imshow(sc);

xlabel(['Size : ',num2str(rx),'*',num2str(cx)]);

title(['Wavelet Decomposition -- Wavelet Type: ',wname,' , Levels: ',num2str(N)]);

function plotcoef(N,wname,coef)

% 画出各级低频、高频系数矩阵的层次结构图

% 首先建立一个名为“Wavelet Decomposition -- Wavelet Type: , Levels: ”的图像窗口figure('Name',['Wavelet Decomposition -- Wavelet Type: ',wname,' , Levels: ',num2str(N)]); % 图像的第1行显示低频系数,置中,左右两个subplot为空

subplot(N+1,3,2);

yt=uint8(coef{1});

[yrow,ycol]=size(yt);

imshow(yt);title( ['Approximation A',num2str(N)]);

xlabel(['Size : ',num2str(yrow),'*',num2str(ycol)]);

% 第2-(N+1)行显示各级高频系数

titllist={['Vertical Detail V'];['Horizontal Detail H'];['Diagonal Detail D']};

pn=2; % pn 是子图的显示序号

for pr=1:N

for pc=1:3

subplot(N+1,3,pn+2);

yt=[];

yt=uint8(coef{pn});

[yrow,ycol]=size(yt);

imshow(yt);title([ titllist{pc},num2str(N-pr+1)]);

xlabel(['Size : ',num2str(yrow),'*',num2str(ycol)]);

% 每行的第一个图像的Y轴,显示该行高频系数对应的分解级别

if mod(pn+2,3)==1

ylabel(['Level ',num2str(N-pr+1)]);

end

pn=pn+1;

end

end

function [o2sa,f1sa,rsx]=sizcoef(sx,lf,N)

% 函数sizcoef 计算小波分解的卷积过程中各分解级的矩阵大小

% 输入参数:sx ——原始(图像)矩阵的大小

% lf ——滤波器的长度

% N ——分解级数

% 输出参数:o2sa ——根据公式sigLen+filLen-1 得到的卷积后矩阵大小% f1sa ——经过下抽样后本级分解系数矩阵的大小

% rsx ——通过f1sa 求出对应于原始矩阵的有效区域大小

% Copyright by Zou Yuhua ( chenyusiyuan ), Version: 1.0, Date: 2008-07-07

o2sa=sx;

f1sa=[0 0];

osx=sx;

for i=1:N

ot=osx+lf-1;

sa=floor(ot/2);

o2sa=[o2sa;ot];

f1sa=[f1sa;sa];

osx=sa;

end

rsx=zeros(N+1,2);

rsx(1,:)=sx;

addln=0;

for j=2:N+1 % 在每级卷积后,序列的长度会增长floor(lf-1/2),aln=floor((addln+lf-1)/2); % 据此可由卷积后的序列长度求出原来的输入序列长度

tsx=f1sa(j,:)-aln;

rsx(j,:)=tsx;

addln=aln;

end

function [cA,cV,cH,cD]=mydwt2(x,Lo_D,Hi_D) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%

% 函数MYDWT2() 对输入的r*c维矩阵x 进行二维小波分解,输出四个分解系数子矩阵[LL,HL,LH,HH]

% 输入参数:x ——输入矩阵,为r*c维矩阵。

% Lo_D,Hi_D ——小波分解的滤波器组系数向量

% 输出参数:cA,cV,cH,cD ——是小波分解系数矩阵的四个相等大小的子矩阵

% cA:低频部分分解系数;cV:垂直方向分解系数;

% cH:水平方向分解系数;cD:对角线方向分解系数。

%

% Copyright by Zou Yuhua ( chenyusiyuan ), Version: 3.0, Date: 2008-07-07 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%

[row,col]=size(x); % 读取输入矩阵的大小

for j=1:row % 首先对输入矩阵的每一行序列进行一维离散小波分解tmp1=x(j,:);

[ca1,cd1]=mydwt(tmp1,Lo_D,Hi_D,1);

% tmp1 长度为row ,滤波器长度为lf

% 则[ca1,cd1] 的总长为2*floor(( row + lf -1 )/2)

x1(j,:)=[ca1,cd1]; % 将分解系数序列存入缓存矩阵x1 中

end

[row1,col1]=size(x1); % row1=row + lnf -1, col1=col+lnf-1

for k=1:col1 % 再对缓存矩阵x1 的每一列序列进行一维离散小波分解tmp2=x1(:,k);

[ca2,cd2]=mydwt(tmp2,Lo_D,Hi_D,1);

x2(:,k)=[ca2,cd2]'; % 将分解所得系数存入缓存矩阵x2 中

% 注意不要遗漏了上一行代码中的转置符号“ ’”。Matlab 6.5 及以下较低的版本不支% 持行、列向量的相互赋值,故要把行向量[ca2,cd2]转置为列向量,再存入x2 的相应列

end

[row2,col2]=size(x2);

cA=x2(1:row2/2,1:col2/2); % cA是矩阵x2的左上角部分

cV=x2(1:row2/2,col2/2+1:col2); % cV是矩阵x2的右上角部分

cH=x2(row2/2+1:row2,1:col2/2); % cH是矩阵x2的左下角部分

cD=x2(row2/2+1:row2,col2/2+1:col2); % cD是矩阵x2的右下角部分

function [cA,cD] = mydwt(x,lpd,hpd,dim) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%

% 函数[cA,cD]=MYDWT(X,LPD,HPD,DIM) 对输入序列x进行一维离散小波分解,输出分解序列[cA,cD]

% 输入参数:x——输入序列;

% lpd——低通滤波器;

% hpd——高通滤波器;

% dim——小波分解层数。

% 输出参数:cA——平均部分的小波分解系数;

% cD——细节部分的小波分解系数。

%

% Copyright by Zou Yuhua ( chenyusiyuan ), Version: 1.0, Date: 2007-11-10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%

cA=x; % 初始化cA,cD

cD=[];

for i=1:dim

cvl=conv(cA,lpd); % 低通滤波,为了提高运行速度,调用MATLAB提供的卷积函数conv() dnl=downspl(cvl); % 通过下抽样求出平均部分的分解系数

cvh=conv(cA,hpd); % 高通滤波

dnh=downspl(cvh); % 通过下抽样求出本层分解后的细节部分系数

cA=dnl; % 下抽样后的平均部分系数进入下一层分解

cD=[cD,dnh]; % 将本层分解所得的细节部分系数存入序列cD

end

function y=downspl(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%

% 函数Y=DOWMSPL(X) 对输入序列进行下抽样,输出序列Y。

% 下抽样是对输入序列取其偶数位,舍弃奇数位。例如x=[x1,x2,x3,x4,x5],则y=[x2,x4].

%

% Copyright by Zou Yuhua ( chenyusiyuan ), Version: 1.0, Date: 2007-11-10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%

N=length(x); % 读取输入序列长度

M=floor(N/2); % 输出序列的长度是输入序列长度的一半(带小数时取整数部分)i=1:M;

y(i)=x(2*i);

原始图像:

小波分解系数塔式结构图:

函数sizcoef 返回的信息:

lf =

16

[o2sa,f1sa,rsx] =

192 196 0 0 192 196 207 211 103 105 96 98 118 120 59 60 48 49

74 75 37 37 24 24

最新小波去噪matlab程序.优选

[转帖]小波去噪matlab程序 ****************************************** clear clc %在噪声环境下语音信号的增强 %语音信号为读入的声音文件 %噪声为正态随机噪声 sound=wavread('c12345.wav'); count1=length(sound); noise=0.05*randn(1,count1); for i=1:count1 signal(i)=sound(i); end for i=1:count1 y(i)=signal(i)+noise(i); end %在小波基'db3'下进行一维离散小波变换[coefs1,coefs2]=dwt(y,'db3'); %[低频高频] count2=length(coefs1); count3=length(coefs2); energy1=sum((abs(coefs1)).^2); energy2=sum((abs(coefs2)).^2);

energy3=energy1+energy2; for i=1:count2 recoefs1(i)=coefs1(i)/energy3; end for i=1:count3 recoefs2(i)=coefs2(i)/energy3; end %低频系数进行语音信号清浊音的判别 zhen=160; count4=fix(count2/zhen); for i=1:count4 n=160*(i-1)+1:160+160*(i-1); s=sound(n); w=hamming(160); sw=s.*w; a=aryule(sw,10); sw=filter(a,1,sw); sw=sw/sum(sw); r=xcorr(sw,'biased'); corr=max(r); %为清音(unvoice)时,输出为1;为浊音(voice)时,输出为0 if corr>=0.8

matlab小波去噪详解

小波去噪 [xd,cxd,lxd]=wden(x,tptr,sorh,scal,n,'wname') 式中: 输入参数x 为需要去噪的信号; 1.tptr :阈值选择标准. 1)无偏似然估计(rigrsure)原则。它是一种基于史坦无偏似然估计(二次方程)原理的自适应阈值选择。对于一个给定的阈值t,得到它的似然估计,再将似然t 最小化,就得到了所选的阈值,它是一种软件阈值估计器。 2)固定阈值(sqtwolog)原则。固定阈值thr2 的计算公式为:thr 2log(n) 2 = (6)式中,n 为信号x(k)的长度。 3)启发式阈值(heursure)原则。它是rigrsure原则和sqtwolog 原则的折中。如果信噪比很小,按rigrsure 原则处理的信号噪声较大,这时采用sqtwolog原则。 4)极值阈值(minimaxi)原则。它采用极大极小原理选择阈值,产生一个最小均方误差的极值,而不是没有误差。 2.sorh :阈值函数选择方式,即软阈值(s) 或硬阈值(h). 3.scal :阈值处理随噪声水平的变化,scal=one 表示不随噪声水平变化,scal=sln 表示根据第一层小波分解的噪声水平估计进行调整,scal=mln 表示根据每一层小波分解的噪声水平估计进行调整. 4.n 和wname 表示利用名为wname 的小波对信号进行n 层分解。输出去噪后的数据xd 及xd 的附加小波分解结构[cxd,lxd]. 常见的几种小波:haar,db,sym,coif,bior haar db db1 db2 db3 db4 db5 db6 db7 db8 db9 db10 sym sym2 sym3 sym4 sym5 sym6 sym7 sym8 coif coif1 coif2 coif3 coif4 coif5 coif6 coif7 coif8 coif9 coif10 bior bior1.1 bior1.3 bior1.5 bior2.2 bior2.4 bior2.6 bior2.8 bior3.5 bior3.7 bior3.9 bior4.4

小波去噪matlab程序

小波去噪matlab程序 ****************************************** clear clc %在噪声环境下语音信号的增强 %语音信号为读入的声音文件 %噪声为正态随机噪声 sound=wavread('c12345.wav'); count1=length(sound); noise=0.05*randn(1,count1); for i=1:count1 signal(i)=sound(i); end for i=1:count1 y(i)=signal(i)+noise(i); end %在小波基'db3'下进行一维离散小波变换 [coefs1,coefs2]=dwt(y,'db3');%[低频高频] count2=length(coefs1); count3=length(coefs2); energy1=sum((abs(coefs1)).^2); energy2=sum((abs(coefs2)).^2); energy3=energy1+energy2; for i=1:count2 recoefs1(i)=coefs1(i)/energy3; end for i=1:count3 recoefs2(i)=coefs2(i)/energy3; end %低频系数进行语音信号清浊音的判别 zhen=160; count4=fix(count2/zhen); for i=1:count4 n=160*(i-1)+1:160+160*(i-1); s=sound(n); w=hamming(160); sw=s.*w; a=aryule(sw,10); sw=filter(a,1,sw);

关于matlab矩阵分解

(1) LU分解 A是非奇异的,LU分解总是可以进行的。 [L,U]=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换),矩阵X必须是方阵。 [L,U,P]=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PX=LU。矩阵X必须是方阵。 实现LU分解后,线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb),这样可以大大提高运算速度。 例7-2 用LU分解求解例7-1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [L,U]=lu(A); x=U\(L\b) 或采用LU分解的第2种格式,命令如下: [L,U ,P]=lu(A); x=U\(L\P*b) (2) QR分解 对矩阵X进行QR分解,就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进 行QR分解,其调用格式为: [Q,R]=qr(X):产生一个一个正交矩阵Q和一个上三角矩阵R,使之满足X=QR。[Q,R,E]=qr(X):产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E,使之满足XE=QR。 实现QR分解后,线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。 例7-3 用QR分解求解例7-1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [Q,R]=qr(A); x=R\(Q\b) 或采用QR分解的第2种格式,命令如下: [Q,R,E]=qr(A); x=E*(R\(Q\b)) (3) Cholesky分解

小波图像去噪及matlab分析

小波图像去噪及matlab实例 图像去噪 图像去噪是信号处理的一个经典问题,传统的去噪方法多采用平均或线性方法进行,常用的是维纳滤波,但是去噪效果不太好(维纳滤波在图像复原中的作用)。 小波去噪 随着小波理论的日益完善,其以自身良好的时频特性在图像去噪领域受到越来越多的关注,开辟了用非线性方法去噪的先河。具体来说,小波能够去噪主要得益于小波变换有如下特点: (1)低熵性。小波系数的稀疏分布,使图像变换后的熵降低。意思是对信号(即图像)进行分解后,有 更多小波基系数趋于0(噪声),而信号主要部分多集中于某些小波基,采用阈值去噪可以更好的保留原 始信号。 (2)多分辨率特性。由于采用了多分辨方法,所以可以非常好地刻画信号的非平稳性,如突变和断点等(例如0-1突变是傅里叶变化无法合理表示的),可以在不同分辨率下根据信号和噪声的分布来消除噪声。(3)去相关性。小波变换可对信号去相关,且噪声在变换后有白化趋势,所以小波域比时域更利于去噪。(4)基函数选择灵活。小波变换可灵活选择基函数,也可根据信号特点和去噪要求选择多带小波和小波 包等(小波包对高频信号再次分解,可提高时频分辨率),对不同场合,选择不同小波基函数。 根据基于小波系数处理方式的不同,常见去噪方法可分为三类: (1)基于小波变换模极大值去噪(信号与噪声模极大值在小波变换下会呈现不同变化趋势)

(2)基于相邻尺度小波系数相关性去噪(噪声在小波变换的各尺度间无明显相关性,信号则相反)(3)基于小波变换阈值去噪 小波阈值去噪是一种简单而实用的方法,应用广泛,因此重点介绍。 阈值函数选择 阈值处理函数分为软阈值和硬阈值,设w是小波系数的大小,wλ是施加阈值后小波系数大小,λ为阈值。(1)硬阈值 当小波系数的绝对值小于给定阈值时,令其为0,而大于阈值时,保持其不变,即: (2)软阈值 当小波系数的绝对值小于给定阈值时,令其为0,大于阈值时,令其都减去阈值,即: 如下图,分别是原始信号,硬阈值处理结果,软阈值处理结果。硬阈值函数在|w| = λ处是不连续的,容易造成去噪后图像在奇异点附近出现明显的伪吉布斯现象。 阈值大小的选取 阈值的选择是离散小波去噪中最关键的一部。在去噪过程中,小波阈值λ起到了决定性作用:如果阈值太小,则施加阈值后的小波系数将包含过多的噪声分量,达不到去噪的效果;反之,阈值太大,则去除了有用的成分,造成失真。小波阈值估计方法很多,这里暂不介绍。 小波去噪实现步骤 (1)二维信号的小波分解。选择一个小波和小波分解的层次N,然后计算信号s到第N层的分解。

Matlab_中的矩阵分解函数

Matlab 中的矩阵分解函数 矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解(三角分解)、QR分解(正交变换)、Cholesky分解,以及Schur分解、Hessenberg分解、奇异分解等。 (1) LU分解(三角分解)lu函数 [L,U]=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换),使之满足X=LU。注意,这里的矩阵X必须是方阵。 实现LU分解后,线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb),这样可以大大提高运算速度。 [L,U,P]=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PX=LU。当然矩阵X同样必须是方阵。(设P 是一个m×n的(0,1) 矩阵,如m≤n且P*P′=E,则称P为一个m×n的置换矩阵。) 例1用LU分解求解例1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]';

x=U\(L\b) 或采用LU分解的第2种格式,命令如下: [L,U ,P]=lu(A); x=U\(L\P*b) (2) QR分解(正交变换) 对矩阵X进行QR分解,就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解,其调用格式为: [Q,R]=qr(X):产生一个一个正交矩阵Q和一个上三角矩阵R,使之满足X=QR。 [Q,R,E]=qr(X):产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E,使之满足XE=QR。 实现QR分解后,线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。例2用QR分解求解例1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]';

matlab小波函数

Matlab小波函数 一、Matlab小波去噪基本原理 1、带噪声的信号一般是由含有噪声的高频信号和原始信号所在的低频 信号。利用多层小波,将高频噪声信号从混合信号中分解出来。 2、选择合适的阈值对图像的高频信号进行量化处理 3、重构小波图像:依据图像小波分解的低频信号与处理之后的高频信 号来重构图像的信息。 二、第二代小波变换 1、构造方法特点: (1)继承了第一代小波的多分辨率的特性。 (2)不依赖fourior变换,直接在时域完成小波变换。 (3)变换之后的系数可以是整数。 (4)图像恢复质量与变换是边界采用何种延拓方式无关。 2、优点:算法简单,速度快,适合并行处理。对内存需求量小,便于DSP 芯片实现、可用于本位操作运算。 3、提升原理:构造紧支集双正交小波 (1)步骤:分裂—预测—更新 (2)分解与重构 三、matlab小波函数库 1、matlab小波通用函数: (1)wavemngr函数【小波管理器(用于小波管理,添加、删除、储存、读取小波)】 wavemngr(‘add’,FN,FSN,WT,NUMS,FILE) wavemngr(‘add’,FN,FSN,WT,NUMS,FILE,B) % 添加小波函数,FN为family name,FSN为family short name WT为小波类型:WT=1表示正交小波,=2表示非正交小波,=3表示带尺度函数的小波,=4表示无尺度函数的小波,=5表示 无尺度函数的复小波。 小波族只有一个小波,则NUMS=“,否则NUMS表示小波参数的字符串 FILE表示文件名 B=[lb ub]指定小波有效支撑的上下界 wavemngr(‘del’,N) %删除小波 wavemngr(‘restore’)/ wavemngr(‘restore’,IN2) %保存原始小波 OUT1= wavemngr(‘read’) %返回小波族的名称 OUT1= wavemngr(‘read’,IN2) %返回所有小波的名称 OUT1= wavemngr(‘read_asc’) %读取wavelets.asc文件并返回小波信息 (2)scal2frq函数【尺度转换频率】 F=scal2frq(A,’wname’,DELTA) %返回由尺度A,小波函数“wname”和采样周期DELTA决定的准 频率。 (3)orthfilt函数【正交小波滤波器组】

第4讲(4)Matlab中的矩阵分解命令

1 4—6 矩阵分解的Matlab 命令 2 (1)矩阵的LU 分解 (2)矩阵的QR 分解(3 )矩阵的Cholesky 分解(4) 矩阵的奇异值分解(5)矩阵的特征值分解(6) 矩阵的Schur 分解(7) 矩阵的Jordan 标准型分解 3(1)矩阵的LU 分解 在Matlab 中用函数lu 来实现矩阵的LU 分解,其命令格式为:[L, U]=lu(X) 说明:U 矩阵为上三角矩阵,满足X=L*U. 4 [L,U,P]=lu(X) 说明:返回的P 矩阵是置换矩阵,矩阵U 是上三角矩阵,矩阵L 满秩矩阵,满足L*U=P*X. 5例4.1 >>a=[1,2,4;2,-1,5;-1,10,4];>> [b,c,p]=lu(a)运行结果:b = 1.0000 0 0-0.5000 1.0000 00.5000 0.2632 1.0000 6 c = 2.0000 -1.0000 5.00000 9.5000 6.50000 0 -0.2105p =0 1 00 0 11 0 0

7 (2)矩阵的QR 分解 在Matlab 中,矩阵的QR 分解可由函数qr 来实现,其常用的调用格式如下:①[B,C]=qr(A) 说明:返回的C 矩阵为上三角矩阵,矩阵B 为满秩矩阵。 [Q,R,E]=qr(A) 说明:返回的矩阵E 是置换矩阵,矩阵R 是上三角矩阵,矩阵Q 是满秩矩阵,上述矩阵满足关系A*E=Q*R. 8 例4.2 >> a=[1,2,4;2,-1,5;-1,10,4];>> [b,c,e]=qr(a)运行结果:b = -0.1952 -0.5068 -0.83970.0976 -0.8619 0.4976-0.9759 0.0152 0.2177 9 c = -10.2470 -4.1964 0.9759 0 -6.2762 -2.24580 0 -0.0622e =0 0 11 0 00 1 0 10 (3 )矩阵的Cholesky 分解 在Matlab 中用函数chol 对矩阵进行Cholesky 分解,函数chol 的调用格式为:1.R=chol(X) 说明:矩阵X 必须是正定矩阵,否则会返回错误信息,返回的矩阵R 是上三角矩阵。2.[R,p]=chol(X) 说明:此调用格式不管矩阵X 是否正定,都不会返回错误信息。如果矩阵X 正定,则返回上三角矩阵R, p 为零;如果矩阵X 非正定,则返回的矩阵R 也是上三角矩阵,但p 为正数。 11例4.4 >> a=[3,-1,1;-1,5,2;1,2,4];>> b=chol(a)运行结果:b = 1.7321 -0.5774 0.57740 2.1602 1.08010 0 1.5811 12 (4) 矩阵的奇异值分解 在Matlab 中,矩阵的奇异值分解则由函数svd 来实现,其调用格式为:[b,c,d]=svd(A) 说明:返回的矩阵b 为左奇异矩阵,矩阵d 为右奇异矩阵,矩阵c 为奇异值矩阵。例4.5 >> a=[1,2,4;2,-1,5;-1,10,4];>> [b,c,d]=svd(a)

五种常用小波基含MATLAB实现

1.给出五种常用小波基的时域和频域波形图。 与标准的傅里叶变换相比,小波分析中使用到的小波函数具有不唯一性,即小波函数(t)ψ 具有多样性。小波分析在工程应用中,一个十分重要的问题就是最优小波基的选择问题,因为用不同的小波基分析同一个问题会产生不同的结果。目前我们主要是通过用小波分析方法处理信号的结果与理论结果的误差来判定小波基的好坏,由此决定小波基。常用小波基有Haar 小波、Daubechies(dbN)小波、Mexican Hat(mexh)小波、Morlet 小波、Meyer 小波等5种。 (1)Haar 小波 Haar 函数是小波分析中最早用到的一个具有紧支撑的正交小波函数,也是最简答的一个小波函数,它是支撑域在[0,1]∈t 围的单个矩形波。 Haar 函数的 定义如下:其他 1212 1 001-1(t)≤≤≤≤?????=ψt t Haar 小波在时域上是不连续的,所以作为基本小波性能不是特别好。但它也有自己的优点,如: 计算简单; (t)ψ不但与t)2(j ψz][j ∈正交,而且与自己的整数位移正交。 因此,在2j a =的多分辨率系统中Haar 小波构成一组最简单的正交归一的小波 族。 ()t ψ的傅里叶变换是: 2/24=sin ()j e a ψ-ΩΩ ΩΩ()j

Haar 小波的时域和频域波形图 -1.5 -1 -0.5 0.5 1 1.5 t haar 时域 x 10 5 1 2 3 4 5 6 75 f haar 频域 i=20; wav = 'haar'; [phi,g1,xval] = wavefun(wav,i); subplot(1,2,1); plot(xval,g1,'-r','LineWidth',1.5); xlabel('t') title('haar 时域'); g2=fft(g1); g3=abs(g2); subplot(1,2,2);plot(g3); xlabel('f') title('haar 频域')

MATLAB 矩阵分解算法大全

(1)LU 分解法程序:function x=solvebyLU(A,b) % 该函数利用LU分解法求线性方程组Ax=b的解 flag=isexist(A,b); %调用第一小节中的isexist函数判断方程组解的情况if flag==0 disp('该方程组无解!'); x=[]; return; else r=rank(A); [m,n]=size(A); [L,U,P]=lu(A); y(1)=b(1); if m>1 for i=2:m y(i)=b(i)-L(i,1:i-1)*y(1:i-1)'; end end y=y'; % 解Ux=y得原方程组的一个特解 x0(r)=y(r)/U(r,r); if r>1 for i=r-1:-1:1 x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i); end end x0=x0'; if flag==1 %若方程组有唯一解 x=x0; return; else %若方程组有无穷多解 format rat; Z=null(A,'r'); %求出对应齐次方程组的基础解系 [mZ,nZ]=size(Z); x0(r+1:n)=0; for i=1:nZ t=sym(char([107 48+i])); k(i)=t; %取k=[k1,k2...,]; end x=x0; for i=1:nZ x=x+k(i)*Z(:,i); %将方程组的通解表示为特解加对应齐次通解形式 end end end (2)矩阵的QR分解法(c语言):

void QR(double a[N][N],double q[N][N],double r1[N][N],int n) /*QR分解*/ { int i,j,k,r,m; double temp,sum,dr,cr,hr; double ur[N],pr[N],wr[N]; double q1[N][N],emp[N][N]; for(i=1;i=ZERO) { sum=0; for(k=r;kZERO)m=-1; else m=1; cr=m*dr; hr=cr*(cr-a[r][r]); for(i=1;ir) ur[i]=a[i][r]; }; for(i=1;i

基于小波去噪matlab程序示例

clear all clc %在噪声环境下语音信号的增强 %语音信号为读入的声音文件 %噪声为正态随机噪声 sound=wavread('c12345.wav'); count1=length(sound); noise=0.05*randn(1,count1); for i=1:count1 signal(i)=sound(i); end for i=1:count1 y(i)=signal(i)+noise(i); end %在小波基'db3'下进行一维离散小波变换 [coefs1,coefs2]=dwt(y,'db3'); %[低频高频] count2=length(coefs1); count3=length(coefs2); energy1=sum((abs(coefs1)).^2); energy2=sum((abs(coefs2)).^2); energy3=energy1+energy2; for i=1:count2 recoefs1(i)=coefs1(i)/energy3; end for i=1:count3 recoefs2(i)=coefs2(i)/energy3; end %低频系数进行语音信号清浊音的判别 zhen=160; count4=fix(count2/zhen); for i=1:count4 n=160*(i-1)+1:160+160*(i-1); s=sound(n); w=hamming(160); sw=s.*w; a=aryule(sw,10); sw=filter(a,1,sw); sw=sw/sum(sw); r=xcorr(sw,'biased'); corr=max(r); %为清音(unvoice)时,输出为1;为浊音(voice)时,输出为0 if corr>=0.8 output1(i)=0; elseif corr<=0.1

小波变换的原理及matlab仿真程序

基于小波变换的信号降噪研究 2 小波分析基本理论 设Ψ(t)∈L 2( R) ( L 2( R) 表示平方可积的实数空间,即能量有限的信号空间) , 其傅立叶变换为Ψ(t)。当Ψ(t)满足条件[4,7]: 2 () R t dw w C ψψ =<∞? (1) 时,我们称Ψ(t)为一个基本小波或母小波,将母小波函数Ψ(t)经伸缩和平移后,就可以得到一个小波序列: ,()( )a b t b t a ψ -= ,,0a b R a ∈≠ (2) 其中a 为伸缩因子,b 为平移因子。 对于任意的函数f(t)∈L 2( R)的连续小波变换为: ,(,),()( )f a b R t b W a b f f t dt a ψψ-=<>= ? (3) 其逆变换为: 211()(,)()f R R t b f t W a b dadb C a a ψ ψ+-= ?? (4) 小波变换的时频窗是可以由伸缩因子a 和平移因子b 来调节的,平移因子b,可以改变窗口在相平面时间轴上的位置,而伸缩因子b 的大小不仅能影响窗口在频率轴上的位置,还能改变窗口的形状。小波变换对不同的频率在时域上的取样步长是可调节的,在低频时,小波变换的时间分辨率较低,频率分辨率较高:在高频时,小波变换的时间分辨率较高,而频率分辨率较低。使用小波变换处理信号时,首先选取适当的小波函数对信号进行分解,其次对分解出的参数进行阈值处理,选取合适的阈值进行分析,最后利用处理后的参数进行逆小波变换,对信号进行重构。 3 小波降噪的原理和方法 3.1 小波降噪原理 从信号学的角度看 ,小波去噪是一个信号滤波的问题。尽管在很大程度上小波去噪可以看成是低通滤波 ,但由于在去噪后 ,还能成功地保留信号特征 ,所以在这一点上又优于传统的低通滤波器。由此可见 ,小波去噪实际上是特征提取和低通滤波的综合 ,其流程框图如 图所示[6] : 小波分析的重要应用之一就是用于信号消噪 ,一个含噪的一维信号模型可表示为如下

矩阵的LU分解(自编MATLAB)实验报告

LU 分解原理 定理:设A C n n ,如果 A 的顺序主子式 A 11 ≠0, |a 11 a 12 a 21 a 22|≠0,…,|a 11a 12a 21a 22…a 12…a 22??a n?11a n?12? ?a n?1n?1 |≠0 则存在唯一的主对角线上元素全为 1 的下三角矩阵L 与唯一的上三角矩阵 U ,使得 A =LU . 证明:对矩阵A 的阶数使用数学归纳法. 显然,当 n=1 时,A 11=1 ?A 11 就是唯一的分解式。现假定对 n-1 阶矩阵,定理的结论成立。对 A 进行分块 A =( A A ?A A A A A A A AA ) 其中A A ,A A ∈A A ?A .由于 n-1 阶矩阵 A A ?A 的 k 阶顺序主子式就是 A 的 k 阶主子式(k=1,2,…,n-2),故它们都不为零.从而由归纳法假设,A A ?A 有唯一的 LU 分解 A A ?A =A A ?A A A ?A 其中A A ?A 的主对角线上的元素都1.由于 |A A ?A |=| A 11A 12A 21A 22 …A 12…A 22 ? ? A A ?11A A ?12 ? ? A A ?1A ?1 |=|A A ?A A A ?A |≠0 所以A A ?A 及A A ?A 是n-1阶可逆矩阵 先假设已有 A =LU ,其中 L =( A A ?A 0A A 1 ), U= ( A A ?A A A A A AA ) A ,A ∈A A ?A 是待定向量。作乘积 AA = (A A ?A A A ?A A A ?A A A A A A ?A A AA +A A A ) =(A A ?A A A A A A A AA )=A 则A,A 必须满足 A A ?A A =A A ,A A A A ?A =A A A ,A AA +A A A =A AA

小波变换图像去噪MATLAB实现

基于小波图像去噪的MATLAB 实现 一、 论文背景 数字图像处理(Digital Image Processing ,DIP)是指用计算机辅助技术对图像信号进行处理的过程。数字图像处理最早出现于 20世纪50年代,随着过去几十年来计算机、网络技术和通信的快速发展,为信号处理这个学科领域的发展奠定了基础,使得DIP 技术成为信息技术中最重要的学科分支之一。在现实生活中,DIP 应用十分广泛,医疗、艺术、军事、航天等图像处理影响着人类生活和工作的各个方面。 然而,在图像的采集、获取、编码和传输的过程中,都存在不同程度被各种噪声所“污染”的现象。如果图像被污染得比较严重,噪声会变成可见的颗粒形状,导致图像质量的严重下降。根据研究表明,当一张图像信噪比(SNR)低于14.2dB 时,图像分割的误检率就高于0.5%,而参数估计的误差高于0.6%。通过一些卓有成效的噪声处理技术后,尽可能地去除图像噪声,我们在从图像中获取信息时就更容易,有利于进一步的对图像进行如特征提取、信号检测和图像压缩等处理。小波变换处理应用于图像去噪外,在其他图像处理领域都有着十分广泛的应用。本论文以小波变换作为分析工具处理图像噪声,研究数字图像的滤波去噪问题,以提高图像质量。 二、 课题原理 1.小波基本原理 在数学上,小波定义为对给定函数局部化的新领域,小波可由一个定义在有限区域的函数()x ψ来构造,()x ψ称为母小波,(mother wavelet )或者叫做基本小波。一组小波基函数,()}{,x b a ψ,可以通过缩放和平移基本小波 来生成: ())(1 ,a b x a x b a -ψ=ψ (1) 其中,a 为进行缩放的缩放参数,反映特定基函数的宽度,b 为进行平移的平移参数,指定沿x 轴平移的位置。当a=2j 和b=ia 的情况下,一维小波基函数序列定义为: ()() 1222,-ψ=ψ--x x j j j i (2) 其中,i 为平移参数,j 为缩放因子,函数f (x )以小波()x ψ为基的连续小波变换定义为函数f (x )和()x b a ,ψ的内积:

matlab小波去噪实现的函数原理

函数wdencmp 功能:小波去噪,得到去噪后的图像 [XC,CXC,LXC,PERF0,PERFL2] = WDENCMP('gbl',X,'wname',N,THR,SORH,KEEPAPP) 其中XC为去噪后的图像信号 在wdencmp中通过xc = waverec2(cxc,lxc,w) ,重构函数得到信号xc Waverec2如何工作的呢? X = W A VEREC2(C,S,'wname') reconstructs the matrix X based on the multi-level wavelet decomposition structure [C,S] 利用经过阈值处理过得系数C和它对应的长度S按照分解时选择的小波来重构;Waverec2涉及到的函数x = appcoef2(c,s,varargin{:},0) Appcoef2函数得到x的方法:x= idwt(a,d,Lo_R,Hi_R,l(imax-p)),综合滤波器重构 Idwt中包含了上采用和卷积函数upsconv1 x = upsconv1(a,Lo_R,lx,dwtEXTM,shift) + upsconv1(d,Hi_R,lx,dwtEXTM,shift); 里面分别调用了采样函数和卷积函数 完成!! 函数wavedec2 功能:返回N层小波分解系数,使用指定滤波器 [C,S] = WA VEDEC2(X,N,'wname') returns the wavelet decomposition of the matrix X at level N,using the wavelet named in string 'wname' ,输出C小波系数,S是对应的系数长度;Wavedec2中通过dwt获得低频系数和小波系数 for i=1:n [x,h,v,d] = dwt2(x,Lo_D,Hi_D); % decomposition c = [h(:)' v(:)' d(:)' c]; % store details s = [size(x);s]; % store size end % Last approximation. c = [x(:)' c]; s = [size(x) ; s]; Dwt2函数如何实现此功能?包含卷积conv2和下采样convdown函数 根据二维mallat变换 输入信号先与滤波器卷积conv2,再下采样得到系数[x,h,v,d] ;

矩阵分解的MATALAB实现

5.3.3 矩阵分解的MATALAB实现 矩阵分解(decomposition, factorization)是多半将矩阵拆解为数个三角形矩阵(triangular matrix),依使用目的的不同,可分为三种矩阵分解法:1)三角分解法(Triangular Factorization),2)QR分解法(QR Factorization),3)奇异值分解法(Singular Value Decompostion)。 (1) 三角分解法 三角分解法是将原正方(square) 矩阵分解成一个上三角形矩阵或是排列(permuted) 的上三角形矩阵和一个下三角形矩阵,这样的分解法又称为LU 分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程,求反矩阵,和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵。 我们举以下二个矩阵为例: 利用三角分解法可将A和B二矩阵分别拆解为上下三角形矩阵 注意B分解的矩阵得到的第一个矩阵[LB]是排列的下三角形矩阵,如果第二、三列互换,则此变成完全的下三角形矩阵。 以MATLAB函数计算上述的LU分解法,其语法为[L,U]=lu(A),其中L代表下三角形矩阵U代表上三角形矩阵。我们来看一个例子。 >> A = [1 2 -1; -2 -5 3; -1 -3 0]; B=[1 3 2; -2 -6 1; 2 5 7]; >> [L1,U1] = lu(A); [L2,U2] = lu(B); >> L1; U1 L1 = % 注意这个矩阵L1和之前的[LA]不相同 -0.5 1 0 1 0 0 0.5 1 1

矩阵的Cholesky分解的Matlab实现

矩阵的Cholesky分解的Matlab实现 下次再补上改进的%Cholesky分解法, %Cholesky分解法 function [X]=m_chol(A,b) [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1); for i=1:N A(i,i)=sqrt(A(i,i)-A(i,1:i-1)*A(i,1:i-1)'); if A(i,i)==0 'A is singular. no unique solution' break end for j=i+1:N A(j,i)=(A(j,i)-A(j,1:i-1)*A(i,1:i-1)')/A(i,i); end end A b %前代法 for j=1:N Y(j)=(b(j)-A(j,1:j-1)*Y(1:j-1))/A(j,j); end Y % A=A' for k=N:-1:1 X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); end Cholesky分解 如果矩阵X是对称正定的,则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R,则下三角矩阵为其转置,即X=R'R。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解,其调用格式为: R=chol(X):产生一个上三角阵R,使R'R=X。若X为非对称正定,则输出一个出错信息。

[R,p]=chol(X):这个命令格式将不输出出错信息。当X为对称正定的,则p=0,R与上述格式得到的结果相同;否则p为一个正整数。如果X 为满秩矩阵,则R为一个阶数为q=p-1的上三角阵,且满足 R'R=X(1:q,1:q)。 实现Cholesky分解后,线性方程组Ax=b变成R‘Rx=b,所以x=R\(R’\b)。例7-4 用Cholesky分解求解例7-1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; R=chol(A) ??? Error using ==> chol Matrix must be positive definite 命令执行时,出现错误信息,说明A为非正定矩阵。 转自: http://203.208.37.104/search?q=cache:EfIWKymoWB8J:www.math.o https://www.doczj.com/doc/d53364852.html,/forums/index.php%3Fact%3DAttach%26type%3Dpost%26id%3 D213920+matlab+%E7%9F%A9%E9%98%B5LU%E5%88%86%E8 %A7%A3&hl=zh-CN&ct=clnk&cd=15&gl=cn&client=firefox-a&st_usg =ALhdy2-W-UGiapmEd7-JkiCNACw5NK2Gew

小波去噪matlab学习指令

MATLAB中实现阈值获取的函数有ddencmp、thselect、wbmpen和wwdcbm,下面对它们的用法进行简单的说明。 ddencmp的调用格式有以下三种: (1)[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,IN2,X) (2)[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,'wp',X) (3)[THR,SORH,KEEPAPP,CRIT]=ddencmp(IN1,'wv',X) 函数ddencmp用于获取信号在消噪或压缩过程中的默认阈值。输入参数X为一维或二维信号;IN1取值为'den'或'cmp','den'表示进行去噪,'cmp'表示进行压缩;IN2取值为'wv'或'wp',wv表示选择小波,wp表示选择小波包。返回值THR是返回的阈值;SORH是软阈值或硬阈值选择参数;KEEPAPP表示保存低频信号;CRIT是熵名(只在选择小波包时使用)。 函数thselect的调用格式如下: THR=thselect(X,TPTR); 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=wbmpen(C,L,SIGMA,ALPHA)返回去噪的全局阈值THR。THR通过给定的一种小波系数选择规则计算得到,小波系数选择规则使用Birge-Massart的处罚算法。{C,L]是进行去噪的信号或图像的小波分解结构;SIGMA是零均值的高斯白噪声的标准偏差;ALPHA是用于处罚的调整参数,它必须是一个大于1的实数,一般去ALPHA=2。 设t*使crit(t)=-sum(c(k)^2,k<=t) + 2 * SIGMA^2 * t*(ALPHA+log(n/t))的最小值,其中c(k)是按绝对值从大到小排列的小波包系数,n是系数的个数,则THR=|c(t*)|。 wbmpen(C,L,SIGMA,ALPHA,ARG)计算阈值并画出三条曲线。 2 * SIGMA^2 * t*(ALPHA+log(n/t)) sum(c(k)^2, k<=t) crit(t) 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层的系数个数。一般

MATLAB中矩阵LU分解

一、 题目 编写实现对N 阶非奇矩阵A 进行LU 分解的程序。 二、 算法组织 若n 阶方阵的各阶顺序主子行列式不为零则存在唯一的单位上三角矩阵L 和上三角矩阵L 式的A=LU 。其基本思想是GAUSS 消去法。参照《计算方法》第38页L 、U 各项计算公式编写公式。 1. 输入带分解矩阵A 2. For i=1,2,……n 2.1 将L 对角线元素赋值L (i ,i )=1; 3. For j=1,2,……n 3.1 将U 第一行元素赋值U (1,j )=A (1,j ); 4. For k=2,……n 4.1 将L 第一列元素赋值L(k,1)=A(k,1)/U(1,1); 5. For i=2,……n 5.1 For j=i ,……n 5.1.110 k kj ki ij kj i A L A U -=- ?∑ 5.2 For k=i+1,……n 5.2.1 1 0k jk ki ij kk jk i A L A U L -=?? -? ??? ∑ 三、 程序实现 clear all clc A=input('请输入一个方阵 ');%输入一个n 阶方阵 [n,n]=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n %将L 的主对角线元素赋值1 L(i,i)=1; end for j=1:n %求矩阵U 的第一行元素 U(1,j)=A(1,j); end for k=2:n %求矩阵L 的第一列元素 L(k,1)=A(k,1)/U(1,1); end for i=2:n %求L 、U 矩阵元素 for j=i:n s=0; for t=1:i-1 s=s+L(i,t)*U(t,j); end U(i,j)=A(i,j)-s; end for k=i+1:n r=0; for t=1:i-1 r=r+L(k,t)*U(t,i); end L(k,i)=(A(k,i)-r)/U(i,i); end end %输出矩阵L 、U L U

相关主题
文本预览
相关文档 最新文档