基于Matlab的序列比对分析3-25
实验四基于Matlab的序列比对分析
实验目的
1．了解MATLAB7.x生物信息工具箱中的序列比对方法；
2．熟悉从数据库获取序列信息, 查找序列的开放阅读框, 将核普酸序列转换为氨基酸序列, 绘制比较两氨基酸序列的散点图, 用Needleman-wunsch算法和Smith-Waterman算法进行比对, 以及计算两序列的同一性的方法；
3．熟悉与序列比对相关的生物信息学函数。
所需软件
MATLAB 7.0或MATLAB 7.0以上的版本
实验内容
序列比对是生物信息学的重要基础。进行序列比对的目的之一是判断两个序列之间是否具有足够的相似性，从而判定二者之间是否具有同源性。序列比对的基本算法主要有两个，一个是用于全局比对的Needleman-Wunsch算法，另一个是主要用于局部比对的Smith-Waterman算法，而后者又是在前者的基础上发展起来的。在MATLAB生物信息工具箱中，序列比对主要用这两种算法。
确定两个序列的相似性是生物信息学的基础工作，通过序列比对（又称序列联配），可以确定两个序列是否具有同源性。
1.查找序列信息
Tay-Sachs症是一种由于缺乏?-氨基已糖苷酶A（Hex A）而导致的常染色体隐性遗传疾病。这种酶能分解大脑和神经细胞中的神经节苷脂（GM2）。基因HEXA编码该酶的?亚基，而第三个基因GM2A编码活化剂蛋白质GM2。
1.1查找目的基因Tay-Sachs
在NCBI（https://www.360docs.net/doc/377465322.html,）上查找信息，在Search列表中选择[Nucleotide],在for框中输入[Tay-Sachs], 点击Go。
1.2读入序列数据
查找结果返回编码酶HexA的α和β亚基的基因和编码活化剂酶的相关页面。NCBI中人类基因HEXA的登录号是NM_000520。用fastaread或genbankread函数可将基因信息被以结构列表的形式导入MATLAB工作区。
方式1：
HumanHEXA = fastaread('NM_000520.txt');
humanHEXA=getfield(HumanHEXA,'Sequence')；
方式2：
HumanHEXA = genbankread('NM_000520.gb')；
humanHEXA=getfield(HumanHEXA,'Sequence')
1.3读入另一序列的信息mouseHEXA
许多基因的序列和功能通过同源基因在进化过程中被保留下来。同源基因就是有共同祖先或是相似序列的基因。查找公共数据库的目的之一就是找出相似的基因。如果用户能在数据库中定位一个未知的基因，那么这个未知基因和已知基因的功能和特征很可能是相同的。
用fastaread或genbankread函数可将鼠类HEXA基因信息被以结构列表的形式导入MATLAB工作区（NCBI中鼠类基因HEXA的序列号是AK080777）。
方式1：
MouseHEXA = fastaread('AK080777.txt');
mouseHEXA=getfield(MouseHEXA,'Sequence')
方式2：
MouseHEXA = genbankread('AK080777.gb')；
mouseHEXA=getfield(MouseHEXA,'Sequence')
2．确定蛋白质编码序列
一个核苷酸序列在蛋白质编码段的前后都包含了调控序列。通过分析这个序
列，可以确定在编码最终蛋白质中亚氨基酸的核苷酸。
2.1 查找人类HEXA的ORF
使用seqshoworfs函数输出人类HEXA的所有阅读框中ORF中起始和终止密码子的位置。
humanORFs = seqshoworfs(humanHEXA)
结果显示了三个阅读框的ORF, 分别以蓝色、红色和绿色标记, 其中最长的ORF在第1个阅读框。
阅读框部分省略
阅读框部分省略
阅读框部分省略
2.2确定鼠类HEXA 的ORF
使用seqshoworfs 函数输出人类HEXA 的所有阅读框中ORF 中起始和终止密码子的位置。
mouseORFs = seqshoworfs(mouseHEXA)
结果得到三个阅读框的ORF, 分别以蓝色、红色和绿色标记, 其中最长的ORF 在第一个阅读框。
Frame 1
阅读框部分省略
阅读框部分省略
阅读框部分省略
3. 比较氨基酸序列
在确定核苷酸序列中的ORF之后，就可以将核苷酸序列的蛋白质编码段转换为相应的氨基酸序列。并使用比对功能来确定两序列的相似性。
3.1 将ORF转换为氨基酸序列
mouseProtein = nt2aa(mouseHEXA);
由于人类的ORF在第一个阅读框, 所以需要指出其位置
humanProtein = nt2aa(humanHEXA,'Frame',1);
3.2 绘制散点图
比较人类和鼠类的氨基酸序列。
seqdotplot(humanProtein,mouseProtein,4,1)
ylabel('Human hexosaminidase A');
xlabel('Mouse hexosaminidase A');
散点图是确定两序列相似性最简单的方法之一。图中对角线平直连续, 表示这两个序列相似性较好。
3.3 比对这两个氨基酸序列
下面nwalign函数有目的地比对两序列。采用的是Needleman-wunsch算法, 可返回全局比对的计算统计量。
[globalscore, globalAlignment] = nwalign(humanProtein,mouseProtein) showalignment(globalAlignment);
Identities = 486/753 (65%), Positives = 570/753 (76%)
3.4 截短序列
寻找终点：
humanStops = find(humanProtein == '*')
mouseStops = find(mouseProtein == '*')
下面将序列截短至只含第一个甲硫氨酸至第一个停止符,进行局部比对。截短序列至只包含蛋白质的氨基酸序列和停止符。
humanSeq = humanProtein(70:humanStops(2));
humanSeqFormatted = seqdisp(humanSeq)
mouseSeq = mouseProtein(11:mouseStops(1));
mouseSeqFormatted = seqdisp(mouseSeq)
3.5 比对被截短的氨基酸序列
[globalscore, globalalignment] = nwalign(humanSeq,mouseSeq);
showalignment(globalalignment);
Identities = 450/540 (83%), Positives = 507/540 (94%)
3.6 局部比对两氨基酸序列
下面swalign函数有目的地比对两序列。采用的是Smith-Waterman算法, 可返回局部比对的计算统计量。
[localscore, localAlignment] = swalign(humanProtein,mouseProtein);
showalignment(localAlignment);
Identities = 454/547 (83%), Positives = 514/547 (94%)
作业
1．进入NCBI任意搜索两条基因组序列，按照序列比对的步骤进行序列比对，进行如下操作并列出结果：
查找序列的开放阅读框, 将核普酸序列转换为氨基酸序列, 绘制比较两氨基酸序列的散点图, 用Needleman-wunsch算法和Smith-Waterman算法进行比对, 以及计算两序列的同一性。
MATLAB中产生高斯白噪声
MATLAB中产生高斯白噪声，涉及到awgn和wgn函数 MATLAB中产生高斯白噪声非常方便，可以直接应用两个函数，一个是WGN，另一个是AWGN。WGN用于产生高斯白噪声，AWGN则用于在某一信号中加入高斯白噪声。 1. WGN：产生高斯白噪声 y = wgn(m,n,p) 产生一个m行n列的高斯白噪声的矩阵，p以dBW为单位指定输出噪声的强度。 y = wgn(m,n,p,imp) 以欧姆(Ohm)为单位指定负载阻抗。 y = wgn(m,n,p,imp,state) 重置RANDN的状态。 在数值变量后还可附加一些标志性参数： y = wgn(…,POWERTYPE) 指定p的单位。POWERTYPE可以是'dBW', 'dBm'或 'linear'。线性强度(linear power)以瓦特(Watt)为单位。 y = wgn(…,OUTPUTTYPE) 指定输出类型。OUTPUTTYPE可以是'real'或 'complex'。 2. AWGN：在某一信号中加入高斯白噪声 y = awgn(x,SNR) 在信号x中加入高斯白噪声。信噪比SNR以dB为单位。x的强度假定为0dBW。如果x是复数，就加入复噪声。 y = awgn(x,SNR,SIGPOWER) 如果SIGPOWER是数值，则其代表以dBW为单位的信号强度；如果SIGPOWER为'measured'，则函数将在加入噪声之前测定信号强度。y = awgn(x,SNR,SIGPOWER,STATE) 重置RANDN的状态。 y = awgn(…,POWERTYPE)指定SNR和SIGPOWER的单位。POWERTYPE可以是'dB'或'linear'。如果POWERTYPE是'dB'，那么SNR以dB为单位，而SIGPOWER以dBW为单位。如果POWERTYPE是'linear'，那么SNR作为比值来度量，而SIGPOWER 以瓦特为单位。 注释 1. 分贝(decibel,dB)：分贝（dB）是表示相对功率或幅度电平的标准单位，换句话说，就是我们用来表示两个能量之间的差别的一种表示单位，它不是一个绝对单位。例如，电子系统中将电压、电流、功率等物理量的强弱通称为电平，电平的单位通常就以分贝表示，即事先取一个电压或电流作为参考值（0dB），用待表示的量与参考值之比取对数，再乘以20作为电平的分贝数（功率的电平值改乘10）。 2. 分贝瓦(dBW, dB Watt)：指以1W的输出功率为基准时，用分贝来测量的功率放大器的功率值。 3. dBm (dB-milliWatt)：即与1milliWatt（毫瓦）作比较得出的数字。 0 dBm = 1 mW 10 dBm = 10 mW 20 dBm = 100 mW 也可直接用randn函数产生高斯分布序列，例如： 程序代码 y=randn(1,2500); y=y/std(y);
matlab实现：常见的离散时间信号
1． 单位抽样序列，或称为离散时间冲激，单位冲激： ? ??=01)(n δ 00≠=n n 如果)(n δ在时间轴上延迟了k 个单位，得到)(k n -δ即： ???=-01)(k n δ 0≠=n k n 2．单位阶跃序列 ? ??01)(n u 00<≥n n 在MATLAB 中可以利用ones( )函数实现。 );,1(N ones x = 3．正弦序列 )(cos )(0φω+=n A n x 这里， ,,0ωA 和φ都是实数，它们分别称为本正弦信号)(n x 的振幅，角频率和初始相位。 πω200=f 为频率。 4．复正弦序列 n j e n x ω=)( 5．实指数序列 n A n x α=)( 6. 随机序列 长度为N 的随机序列 基本数学函数参考教材P69页以及随后的使用说明。 注意使用行向量，特别是冒号运算符。 举例，长度为N 的实指数序列在MATLAB 中实现： n a x N n .^1 :0=-= 1. 单位采样 长度为N 的单位采样序列u(n)可以通过下面的MATLAB 命令获得：
u=[1 )1,1(-N zeros ]； 延迟M 个采样点的长度为N 的单位采样序列ud(n)(M goptions vsize=7cm hsize=10cm; data b; format time monyy5.; input monyy7. asr; dif=dif(asr) ; keep time asr dif; cards; Jan1999 50 Feb1999 54.5 Mar1999 51 Apr1999 49 May1999 50 Jun1999 52 Jul1999 49 Aug1999 49 Sep1999 55 Oct1999 58 Nov1999 60 Dec1999 67.6 Jan2000 62 Feb2000 58.4 Mar2000 55 Apr2000 52.7 May2000 54.4 Jun2000 55.9 Jul2000 53.6 Aug2000 53.4 Sep2000 58.7 Oct2000 62.8 Nov2000 64.2 Dec2000 73.9 Jan2001 66.9 Feb2001 61.7 Mar2001 58.5 Apr2001 56.3 May2001 60.1 Jun2001 60.3 Jul2001 58 Aug2001 58.5 Sep2001 64.3 Oct2001 68.5 Nov2001 70.6 Dec2001 79.2 Jan2002 72.4 Feb2002 67.3 Mar2002 62.9 Apr2002 60.7 May2002 65.9 Jun2002 65.8 Jul2002 62.9 Aug2002 63.6 Sep2002 70.5 Oct2002 76 Nov2002 79 Dec2002 85.1 Jan2003 79.9 Feb2003 73.5 Mar2003 69.5 Apr2003 64.8 May2003 67.6 Jun2003 73.4 Jul2003 70.2 Aug2003 71.6 Sep2003 79.3 Oct2003 85.5 Nov2003 88.5 Dec2003 98.4 Jan2004 90.8 Feb2004 81.8 Mar2004 78.8 Apr2004 75 May2004 81 Jun2004 83.9 Jul2004 80.1 Aug2004 81.1 Sep2004 89.7 Oct2004 98.7 Nov2004 101.7 Dec2004 116.3 Jan2005 103.7 Feb2005 94.2 Mar2005 89.1 Apr2005 86.2 May2005 91.9 Jun2005 98.6 Jul2005 92.2 Aug2005 96.1 Sep2005 103.5 5 常用信号的MATLAB表示5.1单位冲激函数、单位冲激序列 示例7： t = -5:0.01:5; y = (t==0); subplot(121); plot(t, y, 'r'); n = -5:5; x = (n==0); subplot(122); stem(n, x); 图5 运行结果如图5所示。 程序说明： (1)由n = -5:5得到一个1×11数组n；而在x = (n==0)中，n==0是一个向量运算，即向量n中的每一个元素与0比较是否相等，其比较结果0或1放在x中。这样得到的向量x也是1×11数组，且正好就是单位冲激序列。 (2)在MATLAB中，任何向量x的下标是从1开始的，不能取零或负值，而x(n)中的时间变量n则不此受限制。因此向量x的下标与时间变量n是两个概念，如本例中向量x(n)的下标是从1到11，而时间变量n是从-5到5。所以必须用一个与向量x等长的定位时间变量n，以及向量x，才能完整地表示序列x(n)。在信号的表示和运算中，这一点请务必注意；只有当序列x(n)的时间变量正好是从1开始时，才能省去时间变量n，因为此时向量的下标与时间变量相同。 (3)单位冲激函数的实现方法实际上与单位冲激序列是完全相同的，都是用序列表示。只不过表示连续时间信号的序列中两相邻元素所对应的时间间隔更小，如本例中t的间隔为0.01，而表示离散时间信号的序列中两相邻元素所对应的时间间隔一般为1。 由于单位冲激序列在信号与系统中经常使用，我们专门编制一个函数文件delta.m，在后面的实验部分直接调用该函数即可产生需要的波形。 % delta.m function [x, n] = delta(n1,n2,k) % 产生冲激序列δ(n-k),其中n1<=n<=n 2, n1<=k<=n2 2012.3 26 基于MATLAB 的太阳黑子时间序列 分析与仿真 周园 肖洪祥 董俊飞 桂林理工大学信息科学与工程学院 广西 541004 摘要：本文研究了时间序列的分析方法，具体分析了基于最大Lyapunov 指数的方法在太阳黑子时间序列分析中的应用。介绍利用MATLAB 对太阳黑子时间序列进行分析与仿真的方法，并给出相关的流程、程序和相应的仿真结果。最终证明太阳黑子时间序列是一个混沌时间序列。 关键词：混沌时间序列；最大Lyapunov 指数；太阳黑子数；仿真 0 引言 在非线性系统中，初始条件的微小变化，往往会导致结果以指数级的大小发生分离，这时我们称这个系统存在混沌。时间序列是非线性动力系统的一种模型。如果时间序列对初始条件敏感，采用传统线性时间序列分析方法将很难予以分析，因此传统时间序列预测模型对混沌时间序列的拟合和预测准确度都很差。经过混沌学的发展，可以使用序列本身的规律对其进行预测。Lyapunov 指数法即是其中之一。通过最大Lyapunov 指数的数值，可以判断一个时间序列是否是混沌时间序列，亦即该非线性系统中是否存在着混沌。本文对太阳黑子序列进行分析，证明其是一个混沌时间序列。 1 基于Lyapunov 指数的时间序列分析方法 对时间序列进行分析，首先必须进行相空间重构。根据有限的数据重构吸引子以研究系统动力行为的方法即是相空间重构。主要思想为：系统中每个分量的演化皆是由与之联系的其他分量所决定的，相关分量的信息隐含在任意其他分量的变化过程中，即是运用系统的任何一个观察量可以重构出整个系统的模型。 设时间序列为{}t x ，其中1,2,...,t N =。重构相空间m R 的元素组为： (1)(,,)(,,...,),T=1,2,3,...,T T T T m X m N X X X p τττ++-= (1) 其中，N 为重构相空间维数；τ为延迟时间间隔数，且为正整数；(1)p N m τ=--为时间序列嵌入相空间的向量数，N 为时间序列的数据点数。 由Tokens 定理，在理论条件下可任选τ。但在现实条件下时间序列都是有限长且有噪声的。因而在重构相空间时， τ的选取至关重要。目前所采用的方法大多是通过经验来选 择τ， 从而使得T X 和T X τ+相互独立并不完全相关。 Lyapunov 指数是描述奇异吸引子性质的数据量。在m 维离散系统中存在m 个Lyapunov 指数，即Lyapunov 指数族。正的Lyapunov 指数意为在此维度方向，系统以指数级速度分离。1983年，G.Grebogi 证明了若最大Lyapunov 指数 max 0λ>，则系统一定存在着混沌。因此要判断一个时间序 列是否为混沌时间序列，必须求出其最大Lyapunov 指数。为了保证领域点沿着不同的轨道运动，最近邻域点间必须有分离间隔。此处取分离间隔为/w T t =?，其中T 为用FFT 计算出的序列平均周期；t ?为序列的采样周期。 2 计算机仿真步骤 仿真步骤如图1所示。输入太阳黑子年平均序列，通过 FFT 算法计算得到其平均周期T 。计算分离间隔作为时间窗 ωτ。由公式得到嵌入维数m 。运用所得的参数使用Wolf 法 算出最大Lyapunov 指数。进而判断该序列是否是混沌序列。 白噪声的测试MATLAB程序 学术篇 2009-11-13 22:18:03 阅读232 评论0 字号：大中小订阅 clear; clc; %生成各种分布的随机数 x1=unifrnd(-1,1,1,1024);%生成长度为1024的均匀分布 x2=normrnd(0,1,1,1024);%生成长度为1024的正态分布 x3=exprnd(1,1,1024);%生成长度为1024的指数分布均值为零 x4=raylrnd(1,1,1024);%生成长度为1024的瑞利分布 x5=chi2rnd(1,1,1024);%生成长度为1024的kaifang分布%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求各种分布的均值 m1=mean(x1),m2=mean(x2),m3=mean(x3),m4=mean(x4),m5=mean(x5) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求各种分布的方差 v1=var(x1),v2=var(x2),v3=var(x3),v4=var(x4),v5=var(x5) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求各种分布的自相关函数 figure(1);title('自相关函数图'); cor1=xcorr(x1);cor2=xcorr(x2);cor3=xcorr(x3);cor4=xcorr(x4);cor5=xcorr(x5); subplot(3,2,1),plot(1:2047,cor1);title('均匀分布自相关函数图'); subplot(3,2,2),plot(1:2047,cor2);title('正态分布'); subplot(3,2,3),plot(1:2047,cor3);title('指数分布'); subplot(3,2,4),plot(1:2047,cor4);title('瑞利分布'); subplot(3,2,5),plot(1:2047,cor5);title('K方分布'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求各种分布的概率密度函数 y1=unifpdf(x1,-1,1); y2=normpdf(x2,0,1); y3=exppdf(x3,1); y4=raylpdf(x4,1); y5=chi2pdf(x5,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %各种分布的频数直方图 figure(2); subplot(3,2,1),hist(x1);title('均匀分布频数直方图'); subplot(3,2,2),hist(x2,[-4:0.1:4]);title('正态分布'); subplot(3,2,3),hist(x3,[0:.1:20]);title('指数分布'); subplot(3,2,4),hist(x4,[0:0.1:4]);title('瑞利分布'); subplot(3,2,5),hist(x5,[0:0.1:10]);title('K方分布'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %各种分布的概率密度估计 figure(3); 实验1 常见离散信号的MATLAB 产生和图形显示 一、实验目的：加深对常用离散信号的理解 二、实验原理： 1.单位抽样序列：???=01)(n δ 00 ≠=n n 在MATLAB 中可以利用zeros()函数实现。 ; 1)1();,1(==x N zeros x 如果)(n δ在时间轴上延迟了k 个单位，得到)(k n ?δ即：? ??=?01)(k n δ 0≠=n k n 2．单位阶越序列：???0 1 )(n u 00<≥n n 在MATLAB 中可以利用ones()函数实现。(1,)x ones N = 3．正弦序列：)/2sin()(?π+=Fs fn A n x 在MATLAB 中 ) /***2sin(*1 :0fai Fs n f pi A x N n +=?= 4．复正弦序列：n j e n x ?=)( 在MATLAB 中 ) **exp(1 :0n w j x N n =?= 5．指数序列：n a n x =)( 在MATLAB 中 n a x N n .^1:0=?= 三、实验内容： 1、编制程序产生上述5种信号（长度可输入确定），并绘出其图形。 2、讨论正弦序列、复指数序列的性质： （1）绘出信号()zn x n e =，当1126z j π=?+、1126z j π=+、112 z =、62πj z +=、6π j z =时 信号的实部和虚部图；当6 π j z =时信号的周期为多少？ （2）绘出信号() 1.5sin(2*0.1)x n n π=的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期。 3、使用帮助功能学习square(方波),sawtooth(锯齿波)和sinc 函数，并绘图。 四、实验要求： 1、预先阅读MATLAB 基础； 2、讨论复指数序列的性质。 §7.利用SPSS 和Matlab 进行时间序列预测 1.移动平均和滑动平均计算 例1：下表给出了某地区1990～2004年粮食产量数据（表1）。试分别用Matlab 和SPSS 软件，对该地区的粮食产量进行移动平均和和滑动平均计算。 表1 某地区1990～2004年粮食产量及其平滑结果 移动平均 滑动平均 年份 自然序号 粮食产量y （单位：104 t ）三点移动 五点移动 三点滑动 五点滑动 1990 1 3149.44 1991 2 3303.66 3154.47 1992 3 3010.30 3010.30 3141.19 3242.44 1993 4 3109.61 3154.47 3253.04 3263.32 1994 5 3639.21 3141.19 3334.21 3295.88 1995 6 3253.80 3253.04 3242.44 3453.17 3461.80 1996 7 3466.50 3334.21 3263.32 3520.07 3618.81 1997 8 3839.90 3453.17 3295.88 3733.69 3692.89 1998 9 3894.66 3520.07 3461.80 3914.72 3892.78 1999 10 4009.61 3733.69 3618.81 4052.51 4019.78 2000 11 4253.25 3914.72 3692.89 4121.45 4075.78 2001 12 4101.50 4052.51 3892.78 4158.21 4148.58 2002 13 4119.88 4121.45 4019.78 4160.01 4227.01 2003 14 4258.65 4158.21 4075.78 4260.11 2004 15 4401.79 4160.01 4148.58 利用spss 进行移动平均计算主要有以下步骤： （1） 在菜单中依次选择transform->Create time series…,在弹出的对话框中 的单击function 下面的下拉条，选择Prior moving average 方法，span 框中输入数值为3（表示进行三点滑动平均）。 第二章的白噪声产生程序 例2.2 用乘同余法产生(见光盘FLch2bzsheg2.m) ①编程如下： A=6; x0=1; M=255; f=2; N=100；%初始化； x0=1; M=255; for k=1: N %乘同余法递推100次； x2=A*x0; %分别用x2和x0表示x i+1和x i-1； x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（x i）中； v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中； ）减去0.5再乘以存储器f中的系数，存放在v(:,k)=(v1-0.5 )*f; %将v1中的数（ i 矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环 次数变化； x0=x1; % x i-1= x i； v0=v1; end %递推100次结束； v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且 可直接显示在MATLAB的window中； k1=k; %grapher %以下是绘图程序； k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声') ②程序运行结果如图2.6所示。 图2.6 采用MA TLAB产生的(-1,+1)均匀分布的白噪声序列 ③产生的(-1，1)均匀分布的白噪声序列 在程序运行结束后，产生的(-1，1)均匀分布的白噪声序列，直接从MATLAB的window 界面中copy出来如下（v2中每行存6个随机数）： v2 = -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 -0.0156 0.9219 0.5703 0.4531 -0.2500 -0.4844 0.1016 -0.3672 0.8047 -0.1328 0.2188 0.3359 -0.9531 -0.7188 0.6875 -0.8359 *另外，书中图2.3白噪声的产生如下： 显然，只要在例2.2程序的初始化部分中给N=300，f=6，运行程序就可以得到如图2.3所示的（-3，3）的白噪声过程. ①编程如下： A=6; x0=1; M=255; f=6; N=300；%初始化； x0=1; M=255; for k=1: N %乘同余法递推100次； x2=A*x0; %分别用x2和x0表示x i+1和x i-1； x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（x i）中； v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中； ）减去0.5再乘以存储器f中的系数，存放在v(:,k)=(v1-0.5 )*f; %将v1中的数（ i 矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环 次数变化； x0=x1; % x i-1= x i； v0=v1; end %递推100次结束； v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且 可直接显示在MATLAB的window中； k1=k; %grapher %以下是绘图程序； k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声') 时间序列MATLAB程序 时间序列 移动平均法 clc,clear y=[533.8 574.6 606.9 649.8 705.1 772.0 816.4 892.7 963.9 1015.1 1102.7]; m=length(y); n=[4,5];%n为移动平均的项数 for i=1:length(n) %由于n的取值不同，下面使用了细胞数组 for j=1:m-n(i)+1 yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i); end y12(i)=yhat{i}(end);%提出第12月份的预测值 s(i)=sqrt(mean((y(n(i)+1:end)-yhat{i}(1:end-1)).^2));%求预测的标准误差end y12, s %分别显示两种方法的预测值和预测的标准误差 指数平滑 一次指数平滑 程序： clc,clear yt=load('dianqi.txt'); %读取dianqi.txt的数据， n=length(yt); %求yt的长度 alpha=[0.2 0.5 0.8]; %输入a的值 m=length(alpha) yhat(1,[1:m])=(yt(1)+yt(2))/2; %求第一个预测值索引 for i=2:n yhat(i,:)=alpha*yt(i-1)+(1-alpha).*yhat(i-1,:); end yhat %求预测值 err=sqrt(mean((repmat(yt,1,m)-yhat).^2)) %求预测的标准误差 xlswrite('dianqi.xls',yhat) %把预测数据写到Excel文件，准备在word表格中使用yhat1988=alpha*yt(n)+(1-alpha).*yhat(n,:) %求1988的预测值 二次指数平滑 现代通信原理作业一 姓名：张英伟学号：8036 班级：13级理工部3班 利用matlab完成： ●产生正弦波信号、均匀白噪声以及高斯白噪声并分别将两种噪声叠加到正弦 波信号上，绘出波形。 ●分别求取均匀白噪声序列和高斯白噪声序列的自相关及功率谱密度，绘出波 形。 一、白噪声区别及产生方法 1、定义： 均匀白噪声：噪声的幅度分布服从均匀分布，功率谱密度在整个频域内均匀分布的噪声。 高斯白噪声：噪声的幅度分布服从正态分布，功率谱密度在整个频域内均匀分布的噪声。 2、matlab仿真函数： rand函数默认产生是区间在[0,1]的随机数，这里需要利用公式： z2=a+(b-(a))*rand(m,n)............（公式1） randn函数默认产生均值是0、方差是1的随机序列，所以可以用其来产生均值为0、方差为1的正态分布白噪声，即N(0，12)。利用公式： z1=a+b*randn(1,n).................（公式2） 可以产生均值为a，方差为b2 高斯白噪声，即N(a，b2)。 二、自相关函数与功率谱密度之间的关系 1、功率谱密度：每单位频率波携带的功率，这被称为信号的功率谱密度。 2、自相关函数：描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度。 3、维纳-辛钦定理： 由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。幸运的是维纳-辛钦定理提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。 4、平稳随机过程：是在固定时间和位置的概率分布与所有时间和位置的概率分布相同的随机过程。（就是指得仅一个随机过程，中途没有变成另外一个统计特性的随机过程） 小波分析—时间序列的多时间尺度分析 一、问题引入 1.时间序列（Time Series ） 时间序列是指将某种现象某一个统计指标在不同时间上的各个数值，按时间先后顺序排列而形成的序列。在时间序列研究中，时域和频域是常用的两种基本形式。其中： 时域分析具有时间定位能力，但无法得到关于时间序列变化的更多信息； 频域分析（如Fourier 变换）虽具有准确的频率定位功能，但仅适合平稳时间序列分析。 然而，许多现象（如河川径流、地震波、暴雨、洪水等）随时间的变化往往受到多种因素的综合影响，大都属于非平稳序列，它们不但具有趋势性、周期性等特征，还存在随机性、突变性以及“多时间尺度”结构，具有多层次演变规律。对于这类非平稳时间序列的研究，通常需要某一频段对应的时间信息，或某一时段的频域信息。显然，时域分析和频域分析对此均无能为力。 2.多时间尺度 河流因受季节气候和流域地下地质因素的综合作用的影响,就会呈现出时间尺度从日、月到年,甚至到千万年的多时间尺度径流变化特征。推而广之，这个尺度分析，可以运用到对人文历史的认识，以及我们个人生活及人生的思考。 3.小波分析 产生：基于以往对于时间序列分析的各种缺点，融合多时间尺度的理念，小波分析在上世纪80年代应运而生，为更好的研究时间序列问题提供了可能，它能清晰的揭示出隐藏在时间序列中的多种变化周期，充分反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计。 优点： 相对于Fourier 分析：Fourier 分析只考虑时域和频域之间的一对一的映射，它以单个变量（时间或频率）的函数标示信号；小波分析则利用联合时间-尺度函数分析非平稳信号。 相对于时域分析：时域分析在时域平面上标示非平稳信号，小波分析描述非平稳信号虽然也在二维平面上，但不是在时域平面上，而是在所谓的时间尺度平面上，在小波分析中，人们可以在不同尺度上来观测信号这种对信号分析的多尺度观点是小波分析的基本特征。 应用范围： 目前，小波分析理论已在信号处理、图像压缩、模式识别、数值分析和大气科学等众多的非线性科学领域内得到了广泛的应用。在时间序列研究中，小波分析主要用于时间序列的消噪和滤波，突变点的监测和周期成分的识别以及多时间尺度的分析等。 二、小波分析基本原理 1. 小波函数 小波分析的基本思想是用一簇小波函数系来表示或逼近某一信号或函数。因此，小波函数是小波分析的关键，它是指具有震荡性、能够迅速衰减到零的一类函数，即小波函数)R (L )t (2 ∈ψ(有限能量空间)且满足： ?+∞ ∞-=0dt )t (ψ （1） 式中，)t (ψ为基小波函数，它可通过尺度的伸缩和时间轴上的平移构成一簇函数系： )a b t (a )t (2/1b ,a -=-ψψ 其中，0a R,b a,≠∈ （2） 实验名称：常见连续信号的MATLAB 表示 报告人： 姓名班级学号 一、实验目的 1、熟悉常见连续时间信号的意义、特性及波形； 2、学会使用MATLAB 表示连续时间信号的方法； 3、学会使用MATLAB 绘制连续时间信号的波形。 二、实验内容及运行结果 1、运行以上5个例题的程序，保存运行结果。 2、已知信号()t f 的波形如下图所示，试用MATLAB 绘出满足下列要求的信号波形。 <1）()t f -； <2）()2-t f ； <3）()at f <其中a 的值分别为 21= a 和2=a ）； <4）? ?? ??+12 1t f 。 第一题 例题1 程序如下： >> t1=-10:0.5:10。 >> f1=sin(t1>./t1。 >> figure(1> >> plot(t1,f1> >> xlabel('取样间隔p=0.5'>。 >> title('f(t>=Sa(t>=sin(t>/t'>。>> t2=-10:0.1:10。 >> f2=sin(t2>./t2。 >> figure(2> >> plot(t2,f2> >> xlabel('取样间隔p=0.1'>。 >> title('f(t>=Sa(t>=sin(t>/t'>。运行结果如下： f(t)=Sa(t)=sin(t)/t 取样间隔p＝0.5 f(t)=Sa(t)=sin(t)/t 取样间隔p=0.1例题2 程序如下： >> syms t >> f=sin(t>/t。 >> ezplot(f,[-10,10]>运行结果如下： sin(t)/t t 例题3： 程序如下： >> t=-1:0.01:4。 >> t0=0。 >> ut=stepfun(t,t0>。>> plot(t,ut> >> axis([-1,4,-0.5,1.5]>运行结果如下: 时间序列 移动平均法 clc,clear y=[533.8 574.6 606.9 649.8 705.1 772.0 816.4 892.7 963.9 1015.1 1102.7]; m=length(y); n=[4,5];%n为移动平均的项数 for i=1:length(n) %由于n的取值不同，下面使用了细胞数组 for j=1:m-n(i)+1 yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i); end y12(i)=yhat{i}(end);%提出第12月份的预测值 s(i)=sqrt(mean((y(n(i)+1:end)-yhat{i}(1:end-1)).^2));%求预测的标准误差end y12, s %分别显示两种方法的预测值和预测的标准误差 指数平滑 一次指数平滑 程序： clc,clear yt=load('dianqi.txt'); %读取dianqi.txt的数据， n=length(yt); %求yt的长度 alpha=[0.2 0.5 0.8]; %输入a的值 m=length(alpha) yhat(1,[1:m])=(yt(1)+yt(2))/2; %求第一个预测值索引 for i=2:n yhat(i,:)=alpha*yt(i-1)+(1-alpha).*yhat(i-1,:); end yhat %求预测值 err=sqrt(mean((repmat(yt,1,m)-yhat).^2)) %求预测的标准误差 xlswrite('dianqi.xls',yhat) %把预测数据写到Excel文件，准备在word表格中使用yhat1988=alpha*yt(n)+(1-alpha).*yhat(n,:) %求1988的预测值 二次指数平滑 M序列是工程中常用的输入信号，它的性质类似于白噪声，而白噪声是理论上最好的输入信号，可见M序列的价值。下面介绍M序列的matlab产生方法。 idinput函数 产生系统辨识常用的典型信号。 格式 u = idinput(N,type,band,levels) [u,freqs] = idinput(N,'sine',band,levels,sinedata) N 产生的序列的长度，如果N=[N nu]，则nu为输入的通道数，如果N=[P nu M]，则nu 指定通道数，P为周期，M*P为信号长度。默认情况下，nu=1，M=1，即一个通道，一个周期。 Type 指定产生信号的类型，可选类型如下 Band 指定信号的频率成分。对于’rgs’、’rbs’、’sine’，band = [wlow, whigh]指定通带的范围，如果是白噪声信号，则band=[0, 1]，这也是默认值。指定非默认值时，相当于有色噪声。 对于’prbs’，band=[0, B]，B表示信号在一个间隔1/B（时钟周期）内为恒值，默认为[0, 1]。 Levels 指定输入的水平。Levels=[minu, maxu]，在type=’rbs’、’prbs’、’sine’时，表示信号u的值总是在minu和maxu之间。对于type=’rgs’，minu指定信号的均值减标准差，maxu指定信号的均值加标准差，对于0均值、标准差为1的高斯白噪声信号，则levels=[-1, 1]，这也是默认值。 说明 对于PRBS信号，如果M>1，则序列的长度和PRBS周期会做调整，使PRBS的周期为对应一定阶数的最大值（即2^n-1，n为阶数）；如果M=1，PRBS的周期是大于N的相应阶数的值。在多输入的情形时，信号被最大平移，即P/nu为此信号能被估计的模型阶次的上界。 上面的意思可如下理解：对于M=1时， ms = idinput(12, 'prbs', [0 1], [0 1]); figure stairs(ms) title('M序列') ylim([-0.5 1.5]) 基于MATLAB的直接序列扩频通信系统仿真 1.实验原理：直接序列扩频（DSSS）是直接利用具有高码率的扩频码系列采用各种调 制方式在发端与扩展信号的频谱，而在收端，用相同的扩频码序去进行解扩，把扩展宽的扩频信号还原成原始的信息。它是一种数字调制方法，具体说，就是将信源与一定的PN码（伪噪声码）进行摸二加。例如说在发射端将"1"用11000100110，而将"0"用00110010110去代替，这个过程就实现了扩频，而在接收机处只要把收到的序列是11000100110就恢复成"1"是00110010110就恢复成"0"，这就是解扩。这样信源速率就被提高了11倍，同时也使处理增益达到10DB以上，从而有效地提高了整机倍噪比。 1.1 直扩系统模型 直接序列扩频系统是将要发送的信息用伪随机码(PN码)扩展到一个很宽的频带上去，在接收端用与发送端相同的伪随机码对接收到的扩频信号进行相关处理，恢复出发送的信号。对干扰信号而言，与伪随机码不相关，在接收端被扩展，使落入信号通频带的干扰信号功率大大降低，从而提高了相关的输出信噪比，达到了抗干扰的目的。直扩系统一般采用频率调制或相位调制的方式来进行数据调制，在码分多址通信中，其调制多采用BPSK、DPSK、QPSK、MPSK等方式，本实验中采取BPSK方式。 直扩系统的组成如图1所示，与信源输出的信号a(t)是码元持续时间为Ta的信息流，伪随机码产生器产生伪随机码c(t)，每个伪随机码的码元宽度为Tc (Tc< 实验一 常见离散信号的MATLAB 产生和图形显示 授课课时：2学时 一、实验目的： （1）熟悉MATLAB 应用环境，常用窗口的功能和使用方法。 （2）掌握MATLAB 在时域内产生常用离散时间信号的方法。 （3）掌握离散信号的基本运算。 （4）掌握简单的绘图命令。 二、实验原理： （一）信号的表示和产生 ① 单位抽样序列 ? ??=01)(n δ 00 ≠=n n 如果)(n δ在时间轴上延迟了k 个单位，得到)(k n -δ即： ???=-0 1)(k n δ ≠=n k n 参考程序： 例1-1：)2010(()(<<-=n n n x ）δ clear all n1=-10;n2=20;n0=0;%在起点为n1,终点为n2的范围内，于n0处产生冲激。 n=n1:n2;%生成离散信号的时间序列 x=[n==n0];%生成离散信号x(n) stem(n,x);%绘制脉冲杆图 xlabel(' n');ylabel('x(n)');%横坐标和纵坐标的标注说明。 title('Unit Sample Sequence');%图形上方标注图名 axis([-10 20 0 1.2]);%确定横坐标和纵坐标的取值范围 ② 单位阶跃序列 ???=0 1 )(n u 00<≥n n 例1-2：)202((u )(<<-=n n n x ） clear all n1=-2;n2=20;n0=0; n=n1:n2;%生成离散信号的时间序列 x=[n>=n0];%生成离散信号x(n) stem(n,x,'filled'); xlabel('n');ylabel('x(n)'); title('Unit step Sequence'); axis([-2 20 0 1.2]); ③ 正弦序列 )sin()(?+=wn A n x 例1-3：一正弦信号的频率为1HZ ，振幅值幅度A 为1V ，在窗口显示2个周期的信号波形，并对该信号的一个周期进行32点采样获得离散信号并显示该连续信号和离散信号的波形。 参考程序： clear f=1; A=1;nt=2; 现代通信原理作业一 姓名：张英伟学号：133320085208036 班级：13级理工部3班 利用matlab完成： ●产生正弦波信号、均匀白噪声以及高斯白噪声并分别将两种噪声叠加到正弦 波信号上，绘出波形。 ●分别求取均匀白噪声序列和高斯白噪声序列的自相关及功率谱密度，绘出波 形。 一、白噪声区别及产生方法 1、定义： 均匀白噪声：噪声的幅度分布服从均匀分布，功率谱密度在整个频域内均匀分布的噪声。 高斯白噪声：噪声的幅度分布服从正态分布，功率谱密度在整个频域内均匀分布的噪声。 2、matlab仿真函数： rand函数默认产生是区间在[0,1]的随机数，这里需要利用公式： z2=a+(b-(a))*rand(m,n)............（公式1） randn函数默认产生均值是0、方差是1的随机序列，所以可以用其来产生均值为0、方差为1的正态分布白噪声，即N(0，12)。利用公式： z1=a+b*randn(1,n).................（公式2） 可以产生均值为a，方差为b2 高斯白噪声，即N(a，b2)。 二、自相关函数与功率谱密度之间的关系 1、功率谱密度：每单位频率波携带的功率，这被称为信号的功率谱密度。 2、自相关函数：描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度。 3、维纳-辛钦定理： 由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。幸运的是维纳-辛钦定理提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。 4、平稳随机过程：是在固定时间和位置的概率分布与所有时间和位置的概率分布相同的随机过程。（就是指得仅一个随机过程，中途没有变成另外一个统计特性的随机过程）时间序列ARIMA模型的SAS程序编写
常用信号的MATLAB表示
基于MATLAB的太阳黑子时间序列与仿真
白噪声的测试MATLAB程序
实验1 常见离散信号的MATLAB产生和图形显示
利用SPSS和Matlab进行时间序列预测
白噪声产生程序
时间序列MATLAB程序
matlab 正弦波 高斯白噪声 均匀白噪声 功率谱密度 自相关函数
matlab时间序列的多时间尺度小波分析
常见连续信号的MATLAB表示
时间序列MATLAB程序
M序列的matlab产生方法
基于matlab的直接序列扩频通信系统仿真
实验一--常见离散信号的MATLAB产生和图形显示
matlab 正弦波 高斯白噪声 均匀白噪声 功率谱密度 自相关函数