基于Matlab实现的地震波场边界处理软件
基于Matlab实现的地震波场边界处理软件
姓名：姚嘉德学号：2015301130007
院系：资源与环境科学学院
摘要：用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。但我们在实验室进行波场数值模拟时有限差分网格是限制在人工边界里面，即引入了人工边界条件。本文采用Clayton_Engquist_Majda二阶吸收边界条件，通过MATLAB编程实现了这一算法。依靠MATLAB具有更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境，方便同行们交流。利用Matlab本身所具有可视化功能以及像素识别功能，可以将生成的动画电影进行识别，用于地震局实时分析有着深远意义。
关键词：有限差分法，地震波场，吸收边界条件，MATLAB矢量帧,像素识别
Abstract：Modeling seismic wave field with the Finite Difference Method (FDM) is an effective method to study theseismic wave propagation in the earth medium. When we model seismic wave field in the laboratory, the finitedifference grids are restricted in the artificial boundary. So it should introduce the artificial boundary conditions. This paper adopts Clayton_Engquist_Majda second absorbing boundary conditions and realizes the arithmetic with MATLAB.
The MATLAB codes are direct and accord with our thinking custom. So it can provide the friendlyand succinct programming environment and is easy to communicate with https://www.wendangku.net/doc/521826527.html,ing the functions of Matlab that make visualization come true and identify the pixel,we can identify the earthquake wave field.
Key words: finite difference method, seismic wave field, numerical modeling, absorbing boundary conditions,MATLAB
一、引言
用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。但我们在实验室进行波场数值模拟时，只能在有限的空间进行，所以有限差分网格是限制在人工边界里面，即引入了人为的边界条件。这种人为边界条件的引入将对有限区域内的波场值的计算带来严重影响，所以必须进行特殊的边界处理。边界条件处理的好坏直接影响地震正演模拟的最终效果。本文中我们采用Clayton_Engquist_Majda二阶吸收边界条件[2]。
被称作是第四代计算机语言的MATLAB语言，利用其丰富的函数资源把编程工作者从繁琐的程序代码中解放出来。MATLAB用更加直观的、符合大众思维习惯的代码，为用户提供了友好、简洁的程序开发环境。本文介绍运用MATLAB实现带有吸收边界条件的地震波场数值模拟方法和步骤，便于同行们交流，亦可用于本科地震理论的教学中，让学生们在程序演示中理解地震波的传播规律。
二.、Clayton_Engquist_Majda二阶吸收边界条件
我们给定二维标量声波波动方程（含震源）：
（1）
式中：是声波波场，是声波速度，是震源。对（1）
式进行时间和空间2阶精度有限差分离散（见图1），整理后可得
（2）
式中，，为别为空间、时间离散步长，，
，为震源函数。
震源函数：
（3）
Clayton_Engquist_Majda 二阶吸收边界条件的微分表达式可参见文献[2]，其左、右、
上、下 边界的差分格式分别为：
三、基本算法步骤
从图1可以看出，k +1时刻的波场值由k 时刻和k -1时刻的波场值决定。所以在MATLAB
里实现的基本算法步骤如下：
(1) 初始时刻的全波场值均为零，P (i , j , d t )=0（在MATLAB 中初始从d t 开始，不能从0
开始）；
(2) 时刻2d t 时，在炮点S (m , n )给出一个脉冲震源Src (t )（见式（3）），其它点波场P =0；
(3) 时刻t 大于或等于3d t 时，P (i , j , k+1)根据式（2）计算，其它点波场P =0；
(4) 在波传播到四周边界时，左、右、上和下边界的波场值分别由式（4-1）、（4-2）、（4-3）
和（4-4）计算出来。
四、数值模拟
由于是计算机模拟，为了能说明问题且便于计算，我们设地质模型边界为1，具体详细参数如下见表1：
第一步：速度文件的载入及相关整理（替换）
clc; clear; %清除工作空间及显示屏幕
load vm_0.mat; % 载入速度文件，里面包含v(j, i)
Nx=101; Nz=101; Nt=800;hx=0.01;hz=0.01;dt=0.002; % 模拟参数见表1
for i=1:Nx
for j=1:Nz
r(j,i)=v(j,i)*dt/hx;
r2(j,i)=(v(j,i)*dt/hx)^2;
s(j,i)=2-4*(v(j,i)*dt/hx)^2;
end
end % 简缩“常量”
u=zeros(Nz,Nx,Nt); % 分布空间，预值充零
第二步：赋初值
for k=1:2
for j=1:Nz
for i=1:Nx
u(j,i,k)=0;
end
end
end
第三步：边界条件处理及7点差分计算波场延拓
for k=2:Nt-1
for j=1:Nz
for i=2:Nx-1
if j==1
u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j+1,i,k)-r2(j,i)*u(j+2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-…
2*r(j,i)*u(j+1,i,k-1); % 上边界吸收
elseif j==Nz
u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j-1,i,k)-r2(j,i)*u(j-2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-…
2*r(j,i)*u(j-1,i,k-1);
% 下边界吸收
elseif j==30&i==51&(k-1)*dt <= 6*pi/100 %炮点S(30,51)，子波持续时间条件
u(j,i,k+1)=c(j,i)^2*dt^2*sin(50*(k-1)*dt)*exp(-188*((k-1)*dt-3*pi/100)^2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+…
u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);
else
u(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);
end
end
end ；for i=1:Nx 4.1 地质模型的构造
本文选取了两个较为简单的地质模型，分别是均匀介质模型（见图2）和层状均匀介质模型（见图3）。 4.2 程序及相关说明
根据上述我们建立的地质模型，生成相应的速度文件，再联合表 1中的模拟参数和吸收边界条
件，就可以编程实现波场模拟。平时表示波场的习惯u(x,z,t)对应在matlab 程序中，为了便于成图则
被表示为u(z,x,t)，即u(i,j,k)变为u(j,i,k)。详细过程如下：
for j=2:Nz-1
if i==1
u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j,i+1,k)-r2(j,i)*u(j,i+2,k)+(2*r(j,i)-1)*u(j,i,k-1)-…
2*r(j,i)*u(j,i+1,k-1); % 左边界吸收
elseif i==Nx
u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j,i-1,k)-r2(j,i)*u(j,i-2,k)+(2*r(j,i)-1)*u(j,i,k-1)-…
2*r(j,i)*u(j,i-1,k-1); % 右边界吸收
elseif j==30&i==51&(k-1)*dt <= 6*pi/100 % 同上
u(j,i,k+1)=c(j,i)^2*dt^2*sin(50*(k-1)*dt)*exp(-188*((k-1)*dt-3*pi/100)^2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+…
u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);
else
u(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1);
end
end
end
end
第四步：四个角点的处理
for k=1:Nt
u(1,1,k)=1/2*(u(1,2,k)+u(2,1,k));
u(1,Nx,k)=1/2*(u(1,Nx-1,k)+u(2,Nx,k));
u(Nz,1,k)=1/2*(u(Nz-1,1,k)+u(Nz,2,k));
u(Nz,Nx,k)=1/2*(u(Nz,Nx-1,k)+u(Nz-1,Nx,k));
end
第五步：保存结果及成图
u % 在主窗口显示数值结果，按照时间切片
save U u ; % 保存数值结果
% 波场快照图显示
for k=1:Nt % 表示所有时刻，也可以是等间隔时间如k=1:10:Nt
figure(k)
imagesc(u(:,:,k)) ; colormap(gray); % 同理可用surf(u(:,:,k))
xlabel('x'); ylabel('z’');% zlabel('Amplitude') 在surf(u(:,:,k))中用到
title('Wave Field Snapshot');
axis square
end
% 二维电影动画放映，除了imagesc,还有contour、surf等绘图命令
clf; shg;
M=moviein(Nt); %预设画面矩阵
for k=1:Nt
imagesc(u(:,:,k))
colormap(gray)
xlabel('x');ylabel('z');
title('Wave Field Snapshot');
axis([0 Nx 0 Nz]);
M(:,k)=getframe; % 捕获画面
end
% movie(M,5,10) %%电影回放，每秒10帧，重复播放5次
xin_xi_ti_shi=('*****************该程序已经运行结束********************')
4.3 运行结果显示
运行以上程序，可以得到所有时刻的波场快照图，在这里只选择部分图件进行展示。
图4是同一速度均匀介质模型的波场快照图。a、b、c和d是分别200ms、400ms、600ms和800ms 时刻的波场快照，震源位于模型中心，波速为V1。从四个小图上可以看出在波逐步向外传播的过程中，当波传播到边界时逐渐被吸收从而没有人为的反射波。
图5是二层均匀介质模型的波场快照图。a、b、c和d是分别200ms、400ms、600ms和800ms时刻的波场快照，震源位于模型上中部，波速分别为V1和V2。从四个小图上可以看出在波逐步向外传
播的过程中，当波传播到地层边界时有很强的反射，而当波传播到人工边界时逐渐被吸收没有人为
的反射波。
五、成像总结
本文侧重于通过数值模拟介绍运用MATLAB实现带吸收边界条件的波场模拟的简要方法和基本
步骤，因而文中数值模型和其计算差分节点设置都仅考虑十分简单的情形。在实际应用中，除计算
节点数增多外，人们不断通过引入差分格式，差分精度等一系列措施改进正演模拟，从而提高计算
效率，减少计算误差。以上数值模拟程序在MATLAB7.0中运行通过，可供同行们交流。
对于二层均匀介质模型，运行的结果以P(z, x)形式保存成dat文件，然后用地震成图软
六、Matlab的动画实现原理
MATLAB中，创建电影动画的过程分为以下四步：
step1：调用moviein函数对内存进行初始化(该步骤在Matlab5.3以上均可省略)，创建一个足够大的矩阵，使之能够容纳基于当前坐标轴大小的一系列指定的图形（此处称为帧）。
step2：调用getframe函数生成每个帧。该函数返回一个列矢量，利用这个矢量，就可以创建一个电影动画矩阵。
getframe函数可以捕捉动画帧，并保存到矩阵中。一般将该函数放到for循环中得到一系列的动画帧。
该函数格式有：
(1)F=gefframe，从当前图形框中得到动画帧
(2)F=gefframe(h)，从图形句柄h中得到动画帧
(3)F=getframe(h，rect)，从图形句柄h的指定区域rec中得到动画帧
step3：调用movie函数按照指定的速度和次数运行该电影动画。
当创建了一系列的动画帧后，可以利用movie函数播放这些动画帧。
该函数的主要格式有：
(1)movie(M)，将矩阵M中的动画帧播放一次
(2)movie(M,n)，将矩阵M中的动画帧播放n次
(3)movie(M,n,fps)，将矩阵M中的动画帧以每秒fps帧的速度播放n次
step4：调用movie2avi函数可以将矩阵中的一系列动画帧转换成视频文件avi文件。这样，即使脱离了matlab环境都可以播放动画。
参考文献
[1] 孙若昧.地震波传播有限差分模拟的人工边界问题.地球物理学报，1996，11（3）.
[2] 邵秀民,刘臻.带吸收边界条件的声波方程显式差分格式的稳定性分析.计算数学，2001，23（2）.
[3] R.Clayton and B. Enquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bull. Seis. Soc. Am.,
67:6(1977), pp.1529-1540.
[4] R. P. Bording and L. R. Lines, Seismic modeling and imaging with the complete wave equation, USA: SEG, 1997.
[5] 张志涌,徐彦琴.MATLAB教程——基于6.x版本.北京：北京航空航天大学出版社，2001，4.
推荐-Broyden方法求解非线性方程组的Matlab实现 精品
Broyden方法求解非线性方程组的Matlab实现 注：matlab代码来自网络，仅供学习参考。 1.把以下代码复制在一个.m文件上 function [sol, it_hist, ierr] = brsola(x,f,tol, parms) % Broyden's Method solver, globally convergent % solver for f(x) = 0, Armijo rule, one vector storage % % This code es with no guarantee or warranty of any kind. % % function [sol, it_hist, ierr] = brsola(x,f,tol,parms) % % inputs: % initial iterate = x % function = f % tol = [atol, rtol] relative/absolute % error tolerances for the nonlinear iteration % parms = [maxit, maxdim] % maxit = maxmium number of nonlinear iterations % default = 40 % maxdim = maximum number of Broyden iterations % before restart, so maxdim-1 vectors are % stored % default = 40 % % output: % sol = solution % it_hist(maxit,3) = scaled l2 norms of nonlinear residuals % for the iteration, number function evaluations, % and number of steplength reductions % ierr = 0 upon successful termination % ierr = 1 if after maxit iterations % the termination criterion is not satsified. % ierr = 2 failure in the line search. The iteration % is terminated if too many steplength reductions % are taken. % % % internal parameter: % debug = turns on/off iteration statistics display as % the iteration progresses
MATLAB自适应滤波去噪
《MATLAB自适应滤波去噪》 课程设计报告 1.课程设计目的 此次课程设计目的是为了让我们学会使用MATLAB进行计算机仿真，使用自适应滤波法设计一个语音去噪声电路。培养我们的电路设计思路及其算法，明白理论与实践相结合的重要性，培养了我们的实际操作能力以及锻炼我们对实际问题的分析与解决的能力。 2.课程设计内容 2.1 LMS自适应算法原理 自适应过程一般采用典型LMS自适应算法，但当滤波器的输入信号为有色随机过程时，特别是当输入信号为高度相关时，这种算法收敛速度要下降许多，这主要是因为输入信号的自相关矩阵特征值的分散程度加剧将导致算法收敛性能的恶化和稳态误差的增大。此时若采用变换域算法可以增加算法收敛速度。变换域算法的基本思想是：先对输入信号进行一次正交变换以去除或衰减其相关性，然后将变换后的信号加到自适应滤波器以实现滤波处理，从而改善相关矩阵的条件数。因为离散傅立叶变换DFT 本身具有近似正交性，加之有FFT快速算法，故频域分块LMS FBLMS 算法被广泛应用。 FBLMS算法本质上是以频域来实现时域分块LMS算法的，即将时域数据分组构成N个点的数据块，且在每块上滤波权系数保持不变。其原理框图如图2所示。FBLMS 算法在频域内可以用数字信号处理中的重叠保留法来实现，其计算量比时域法大为减少，也可以用重叠相加法来计算，但这种算法比重叠保留法需要较大的计算量。块数据的任何重叠比例都是可行的，但以50％的重叠计算效率为最高。对FBLMS算法和典
型LMS算法的运算量做了比较，并从理论上讨论了两个算法中乘法部分的运算量。本文从实际工程出发，详细分析了两个算法中乘法和加法的总运算量，其结果为： 复杂度之比＝FBLMS实数乘加次数／LMS实数乘加次数＝（25Nlog2N＋2N－4）／[2N(2N－1)] 采用ADSP的C语言来实现FBLMS算法的程序如下： for(i＝0;i＜＝30;i＋＋) {for(j＝0;j＜＝n－1;j＋＋) {in[j]＝input[i×N＋j;] rfft(in,tin,nf,wfft,wst,n); rfft(w,tw,wf,wfft,wst,n); cvecvmlt(inf,wf,inw,n); ifft(inw,t,O,wfft,wst,n); for(j＝0,j＜＝N－1;j＋＋) {y[i×N＋j]＝O[N＋j]．re; e[i×N＋j]＝re fere[i×N＋j]－y[i×N＋j]; temp[N＋j]＝e[i×N＋j;}
MATLAB代码 解线性方程组的迭代法
解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;
if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵
自适应滤波LMS与RLS地matlab实现
MATLAB 仿真实现LMS 和RLS 算法的二阶AR 模型 及仿真结果分析 一、题目概述：二阶AR 模型如图1a 所示，可以如下差分方程表示： )()()2()1()()(21n d n v n x a n x a n v n x +=----= （1） 图1a 其中，v(n)是均值为0、方差为0.965的高斯白噪声序列。 ， 为描述性参数， .95,0,195.021=-=a a 设x(-1)=x(-2)=0，权值 ， =0.04①推导最优 滤波权值（理论分析一下）。②按此参数设置，由计算机仿真模拟权值收敛曲线并画出，改变步长在此模拟权值变化规律。③对仿真结果进行说明。④应用RLS 算法再次模拟最优滤波权值。 解答思路： （1）高斯白噪声用normrnd 函数产生均值为0、方差为0.965的正态分布随机1*N 矩阵来实现。随后的产生的信号用题目中的二阶AR 模型根据公式（1）产生，激励源是之前产生的高斯白噪声。 （2）信号长度N 取为2000点，用以观察滤波器权值变化从而估计滤波器系数，得到其收敛值。 （3）仿真时分别仿真了单次LMS 算法和RLS 算法下的收敛性能以及100次取平均后的LMS 和RLS 算法的收敛性能，以便更好的比较观察二者的特性。 （4）在用不同的分别取3个不同的 值仿真LMS 算法时， 值分别取为
0.001,0.003,0.006；用3个不同的λ值仿真RLS 算法时λ值分别取为1,0.98,0.94，从而分析不同步长因子、不同遗忘因子对相应算法收敛效果的影响。 二、 算法简介 1．自适应算法的基本原理 自适应算法的基本信号关系如下图所示： 图 1b 自适应滤波器框图 输入信号x(n)通过参数可调的数字滤波器后产生输出信号y(n)，将其与参考信号d(n)进行比较，形成误差信号e(n)。e(n)通过某种自适应算法对滤波器参数进行调整，最终是e(n)的均方值最小。当误差信号e(n)的均方误差达到最小的时候，可以证明信号y(n)是信号d(n)的最佳估计。 2. LMS 算法简介 LMS 算法采用平方误差最小的原则代替最小均方误差最小的原则，信号基本关系如下： 1 ()()() ()()()(1)()2()()(0,1,2,....1) N i i i y n w n x n i e n d n y n w n w n e n x n i i N μ-=-=-+=+-=-∑ （2） 写成矩阵型式为： ()()()()()() (1)()2()() T y n W n X n e n d n y n W n W n e n X n μ==-+=+ （3） 式（3）中，W(n) 为n 时刻自适应滤波器的权值， 011()[(),(),....()]T N W n w n w n w n -=，N 为自适应滤波器的阶数，本设计中取为2000；X( n) 为n 时刻自适应滤波器的参考输入矢量，由最近N 个信号采样值构成，()[(),(1),....(1)]T X n x n x n x n N =--+；d ( n) 是期望的输出值；e ( n) 为自适应滤波器的输出误差调节信号(简称失调信号) ；μ是控制自适应速度与稳定性的增益常数，又叫收敛因子或步长因子。
MatLab求解线性方程组
MatLab解线性方程组一文通 当齐次线性方程AX=0,rank(A)=r 1.均匀线阵方向图 （1）matlab 程序 clc; clear all; close all; imag=sqrt(-1); element_num=32;%阵元数为8 d_lamda=1/2;%阵元间距d与波长lamda的关系 theta=linspace(-pi/2,pi/2,200); theta0=0;%来波方向 w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]'); for j=1:length(theta) a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]'); p(j)=w'*a; end patternmag=abs(p); patternmagnorm=patternmag/max(max(patternmag)); patterndB=20*log10(patternmag); patterndBnorm=20*log10(patternmagnorm); figure(1) plot(theta*180/pi,patternmag); grid on; xlabel('theta/radian') ylabel('amplitude/dB') title([num2str(element_num) '阵元均匀线阵方向图','来波方向为' num2str(theta0*180/pi) '度']); hold on; figure(2) plot(theta,patterndBnorm,'r'); grid on; xlabel('theta/radian') ylabel('amplitude/dB') title([num2str(element_num) '阵元均匀线阵方向图','来波方向为' num2str(theta0*180/pi) '度']); axis([-1.5 1.5 -50 0]); 在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法. 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 成都理工大学工程技术学院本科毕业论文 自适应载波同步及其Matlab仿真 作者姓名： 专业名称： 指导老师： 年月日 摘要 自适应滤波算法的研究是现在社会自适应信号处理中最为活跃的研究课题之一。找寻收敛速度快，计算简单，数值稳定性好的自适应滤波算法是研究人员不断努力追求的目标。本设计在论述自适应滤波基本原理的基础上，说明了几种当前几种典型的自适应滤波算法和应用。并对这几种典型自适应滤波算法的性能特点进行简单的比较，给出了算法性能的综合评价。 载波同步是无线通信接收机的主要功能之一，其对通信系统质量的提高至关重要。随着新算法涌现和芯片处理速度的提高，不同的解决方案不断的提出。自适应载波同步是一种依据自适应算法的同步方法，内容新颖。本课题在介绍自适应算法和载波同步问题的基础上，详细讨论了平方差分环路法和锁相环路法，具体包括代价函数、代价函数的导数、迭代公式和原理图等，并在论文的第三部分给出了这两种方法的Matlab仿真。仿真结果验证了这两种方法在跟踪载波相位方面是满足要求的，且收敛速度较快。 关键词：自适应滤波载波同步平方差分环路锁相环路法 Abstact The research of adaptive filtering algorithm is one of the most activity tasks, the goal that researchers want to pursue is to find an adaptive filtering algorithm that converge fast and compute simplely. Based on the basis adaptive filtering principle, this paper introduces several typical adaptive algorithms and applications, then compares those algorithm's characters and gives the orithm performance evaluation. Carrier synchronization is one of the main functions of Wireless communications receiver，it is essential for the improvement in the quality of the communication system. With the emergence of new algorithms and the speed improvement of chip processing, different solutions is proposed continuously. Adaptive carrier synchronization is a synchronization method based on adaptive algorithms, and its content is innovative. Based on the introducing of adaptive algorithm and carrier synchronization, this issue has a detailed discussion of the square difference method and the PLL loop method, including its cost function, cost function derivative, iterative formula and schematic, etc. And the third part of the paper gives two methods of Matlab simulation.Simulation results show the two methods with tracking the carrier phase is to meet the requirements, and convergence speedly. Keywords:adaptive filter, carrier synchronization, differential circle square , phase-locked loop method 《MATLAB语言》课成论文 利用MATLAB求线性方程组 姓名：郭亚兰 学号：12010245331 专业：通信工程 班级：2010级通信工程一班 指导老师：汤全武 学院：物电学院 完成日期：2011年12月17日 利用MATLAB求解线性方程组 （郭亚兰 12010245331 2010 级通信一班） 【摘要】在高等数学及线性代数中涉及许多的数值问题，未知数的求解，微积分，不定积分，线性方程组的求解等对其手工求解都是比较复杂，而MATLAB语言正是处理线性方程组的求解的很好工具。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是讨论矩阵理论、与矩阵结合的有限维向量空间及其线性变换理论的一门学科。 【关键字】线性代数MATLAB语言秩矩阵解 一、基本概念 1、N级行列式A：A等于所有取自不同性不同列的n个元素的积的代数和。 2、矩阵B：矩阵的概念是很直观的，可以说是一张表。 3、线性无关：一向量组（a1,a2,…,an）不线性相关，既没有不全为零的数 k1,k2,………kn使得：k1*a1+k2*a2+………+kn*an=0 4、秩：向量组的极在线性无关组所含向量的个数成为这个向量组的秩。 5、矩阵B的秩：行秩，指矩阵的行向量组的秩；列秩类似。记：R(B) 求解线性方程组 solve，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff（fun，var，n）：对表达式fun中的变量var求n阶导数。 例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式 非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名；fun其中 x0位求解方程的初始向量或矩阵； option为设置命令参数 建立文件fun.m： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注： ...为续行符 m文件必须以function为文件头，调用符为@；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。如： X=A\B表示求矩阵方程AX＝B的解； 的解。XA=B表示矩阵方程B/A＝X． 对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。 如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； m %用matlab实现自适应图像阈值分割最大类方差法代码clear; warning off; SE = strel('diamond',4); BW1 = imread('cameraman.tif'); BW2 = imerode(BW1,SE); BW3 = imdilate(BW2,SE); BW4 = BW1-BW3; %rgb转灰度 if isrgb(BW4)==1 I_gray=rgb2gray(BW4); else I_gray=BW4; end figure,imshow(I_gray); I_double=double(I_gray);%转化为双精度 [wid,len]=size(I_gray); colorlevel=256; %灰度级 hist=zeros(colorlevel,1);%直方图 %threshold=128; %初始阈值 %计算直方图 for i=1:wid for j=1:len m=I_gray(i,j)+1; hist(m)=hist(m)+1; end end hist=hist/(wid*len);%直方图归一化 miuT=0; for m=1:colorlevel miuT=miuT+(m-1)*hist(m); end xigmaB2=0; for mindex=1:colorlevel threshold=mindex-1; omega1=0; omega2=0; for m=1:threshold-1 omega1=omega1+hist(m); end omega2=1-omega1; miu1=0; miu2=0; for m=1:colorlevel if m 第18卷第7期 系统仿真学报?V ol. 18 No. 7 2006年7月Journal of System Simulation July, 2006 基于MATLAB/Simulink机器人鲁棒自适应控制系统仿真研究 高道祥，薛定宇 （东北大学教育部暨辽宁省流程工业综合自动化重点实验室，沈阳 110004） 摘要：介绍了一种在MATLAB/Simulink环境下进行机器人鲁棒自适应控制系统仿真的方法，利 用Matlab软件强大的数值运算功能，将系统模型用Matlab语言编写成M-Function（或S-Function） 文件，通过User-Defined-Function模块嵌入到Simulink仿真环境中，可以充分发挥Simulink模块 实时的动态仿真功能，简化仿真模型的设计，修改和调整。基于M-Function建立机器人系统模型 的方法可以推广到其他复杂控制系统的建模，SimMechanics在建立多自由度连杆机器人受控对象 仿真模型时，简单可靠。 关键词：机器人；Matlab/Simulink；SimMechanics；仿真；鲁棒自适应控制 中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X(2006) 07-2022-04 Simulation Research of Robust Adaptive Control System for Robotic Manipulators Based on MATLAB/Simulink GAO Dao-xiang, XUE Ding-yü (Key Laboratory of Process Industry Automation, Ministry of Education, Northeast University, Shenyang 110004, China) Abstract: A simulation method of robust adaptive control was proposed for the robotic manipulator system. The method took the advantage of the powerful computing function of Matlab to programme M-function (or S-Function) for the system model by Matlab language and embedded it to the Simulink by User-Defined-Function module. The real time dynamic simulating function of Simulink would be exerted adequately and the design, modification and adjust of the system model could be greatly simplified. The method of constructing manipulator control system model based on M-Function could be generalized to the other complicated control system and SimMechanics would make the n-links manipulator model conveniently and credibly. Key words: robotic manipulator; Matlab/Simulink; SimMechanics; simulation; robust adaptive control 引言 一个新的控制算法在付诸使用之前，无论从经济原因还是技术角度，都需要经过仿真阶段来测试控制系统的性能和缺陷。尤其对复杂系统控制的研究，虽然仿真并不能说明控制算法是绝对合理与可靠的，但随着仿真技术的发展，仿真的确是系统设计必不可少的中间步骤。 Matlab/Simulink以其强大的数学运算能力，方便实用的绘图功能及语言的高度集成性成为系统仿真和自动控制领域首选的计算机辅助设计工具。Simulink可以将可视化的模块很方便地组成系统模型的仿真框图，对于使用普通Simulink模块不易搭建的复杂控制系统，用Matlab语言编写M-Function或S-Function文件，通过User-Defined-Function 模块嵌入到Simulink仿真环境中，大大扩充了Simulink的功能。对于机器人这类的复杂控制系统，通过Simulink可以很方便的建立其仿真模型。 机器人控制系统仿真模型中不易采用普通Simulink模块搭建的部分是控制器模型和受控对象—机器人模型，可以采用Matlab语言编写M-Function实现控制器和机器人模型。 收稿日期：2005-05-09 修回日期：2005-08-02 作者简介：高道祥（1972-），男，山东蓬莱人，博士生，研究方向为机器人鲁棒自适应控制。薛定宇（1963-），男，辽宁沈阳人，教授，博导，研究方向为控制系统CAD，机器人控制。另外，由于SimMechanics提供了机构的仿真模块集，对于n自由度的连杆机器人，也可以采用SimMechanics模块进行组建。 1 n连杆机器人的仿真模型 如果不考虑摩擦力等外界干扰的作用，机器人的动力学方程可由下式描述， τ= + +) ( ) , ( ) (q G q q q C q q M (1) 式中，q q q , ,是1 × n向量，表示各个关节的位置，速度，加速度。) (q M是n n×阶对称正定的惯量矩阵。q q q C ) , (是1 × n向量，表示离心力和哥氏力项。) (q G是1 × n向量，表示重力项。τ表示外界输入的控制力矩。由式(1)可以看出，机器人的动力学模型是一个高度复杂，强耦合的非线性时变方程，尤其模型的复杂程度随着自由度的增加呈指数上升趋势。因此，在用Matlab/Simulink进行机器人控制系统的仿真研究时，需要寻求一种简单可靠却行之有效的方法建立机器人控制系统仿真模型。 采用M-Function定制的Simulink模块与普通模块一样具有输入和输出向量，控制器与机器人受控对象的仿真模型函数可用如下关系式描述， ) , , , , (0q q q q q f d d d = τ(2) ) , , (1τ q q f q =(3) 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 第7章 求解非线性方程 7.1 多项式运算在MATLAB 中的实现 一、多项式的表达 n 次多项式表达为：n a +??++=x a x a x a p(x )1-n 1-n 1n 0，是n+1项之和 在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示 [a0, a1,……an-1,an] 二、多项式的加减运算 设 有 两 个 多 项 式 n a +??++=x a x a x a p1(x )1-n 1-n 1n 0和 m b +??++=x b x b x b p2(x )1-m 1-m 1m 0。它们的加减运算实际上就是它们的对应系 数的加减运算。当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。 例2 计算()()1635223-+++-x x x x a=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b 例3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐 f+g1, f-g1 三、多项式的乘法运算 conv(p1,p2) 例4 在上例中，求f(x)*g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g) 四、多项式的除法运算 [Q, r]=deconv(p1, p2) 表示p1除以p2，给出商式Q(x)，余式r(x)。Q,和r 仍为多项式系数向量 例4 在上例中，求f(x)/g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数 p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P·Q 的导函数 自适应滤波器 MATLAB仿真 摘要 : 本文介绍了自适应滤波器的工作原理，以及推导了著名的LMS( Least mean squares )算法。以一个例子演示了自适应滤波器的滤波效果。实验结果表明，该滤波器滤波效果较好。 关键词:自适应滤波器 MATLAB7.0 LMS 算法 Simulate of adaptive filter based on MATLAB7.0 Abstract: This article described the working principle of adaptive filter and deduced the well-known LMS algorithm. Take an example to demonstrate the adaptive filters filtering effects. The results show that the filter has an effective way to filter single. Key words: LMS algorithm Adaptive Filter Matlab7.0 1引言 由 Widrow B 等提出的自适应滤波理论，是在维纳滤波、卡尔曼滤波等线性滤波基础上发展起来的一种最佳滤波方法。由于它具有更强的适应性和更优的滤波性能，从而广泛应用于通信、系统辨识、回波消除、自适应谱线增强、自适应信道均衡、语音线性预测和自适应天线阵等诸多领域[1]。自适应滤波器最大的优点在于不需要知道信号和噪声的统计特性的先验知识就可以实现信号的最佳滤波处理。本文通过一个具体例子和结果论证了自适应滤波器的滤波效果。 2自适应滤波原理及 LMS算法 2.1 自适应滤波原理 图 1 自适应滤波原理图 在自适应滤波器中，参数可调的数字滤波器一般为 FIR 数字滤波器， IIR 数字滤波器或格型数字滤波器。自适应滤波分 2 个过程。第一，输入信号想 x(n) 通过参数可调的数字滤波器后得输出信号 y(n) ，y(n) 与参考信号 d(n) 进行比较得误差信号 e(n) ；第二，通过一种自适应算法和 x(n) 和 e(n) 的值来调节参数可调的数字滤波器的参数，即加权系 实验1.1 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref （A ） %将矩阵A 化为阶梯形的最简式 null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基 础解系 【例１】 求下列齐次线性方程组的一个基础解系，并写出通解： 我们可以通过两种方法来解： 解法1： >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x 取x2，x4为自由未知量，扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系，记 所以齐次方程组的通解为 解法2： clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用, 若去掉r ，是什么结果？ 执行后可得结果： B= 1 0 1 0 0 1 0 1 ?? ?=-=-0 04321x x x x ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x +=自适应波束形成与Matlab程序代码注解
MATLAB解线性方程组的直接方法
自适应载波同步及其Matlab仿真
利用MATLAB求线性方程组
Matlab求解线性方程组非线性方程组
用matlab实现自适应图像阈值分割最大类方差法代码
基于MATLAB_Simulink机器人鲁棒自适应控制系统仿真研究
线性方程组求解matlab实现
MATLAB应用 求解非线性方程
自适应滤波器MATLAB仿真
实验一用matlab求解线性方程组
- Matlab求解线性方程组非线性方程组
- 伪弧长法求解非线性方程matlab程序
- Matlab求解线性方程组、非线性方程
- matlab程序设计实践 牛顿法解非线性方程
- 推荐-Broyden方法求解非线性方程组的Matlab实现 精品
- matlab非线性方程求解
- matlab程序设计实践-牛顿法解非线性方程
- 用matlab求解非线性方程组的几种方法之程序
- 实验2 利用matlab解(非)线性、微分方程(组)
- 基于matlab的非线性方程组求解的方法
- MATLAB解方程
- 第二章非线性方程(组)的数值解法的matlab程序
- 牛顿法求解非线性方程组matlab源程序
- Matlab求解线性方程组、非线性方程组
- ——matlab——二分法求解非线性方程的近似解
- MATLAB应用 求解非线性方程
- 用matlab对非线性方程求解
- 非线性方程组求解及matlab实现分解
- MATLAB求解非线性方程
- 非线性方程组求解及matlab实现