基于MATLAB的信号合成与分解.doc

计算机与信息学院2011-2012-课程论文 MATLAB与信号处理系统课程设计课程性质：考察学号：2010156203姓名：蔡红亮专业：通信工程授课教师：孙水发完成日期：2012年6月25日 周期信号的分解与合成 作者：蔡红亮摘要：为了便于进行周期信号的分析与处理，常要把复杂的周期信号进行分解，即将周期信号分解为正余弦等此类基本信号的线性组合，通过对这些基本信号单元在时域和频域特性的分析来达到了解信号特性的 目的。本文主要阐述了傅立叶级数的推演过程，从而得出周期信号的分解与合成的基本原理。并利用Matlab仿真软件强大的数值分析和图形功能来对周期方波信号与周期三角波信号以及周期锯齿波信号的分解与合成进行演示，直观明了的观察周期信号分解与合成过程、周期信号的对称性与谐波成分的关系，以及对其中的误差程度和吉布斯现象做定量的分析，从而可以进行仿真结果与理论分析结论的对比，加深了对周期连续信号分解与合成的理解，描述了傅立叶级数分解合成信号的实现性，同时也展示了用MATLAB编写周期连续信号分解与合成的演示程序的优点关键词： 周期信号 分解 合成 吉布斯现象 正文：一、前言1、 对周期信号时频分析的重要意义 信号分析通过研究信号的描述、运算、特性以及信号发生某些变化时其特性相应的变化，来揭示信号自身的时域特性、频域特性等。信号分析的主要途径是研究信号的分解，即将信号分解为某些基本信号的线性组合，通过对这些基本信号单元在时域和频域特性的分析来达到了解信号特性的目的。信号的分解可以在时域、频域或变换域中进行，分别用到信号分析的时域方法（time domain analysis）、频域方法（frequency domain analysis）和变换域方法（transform domain analysis）。系统是若干相互依赖、相互作用的事物组合而成的、具有特定功能的整体。系统可以是物理系统，例如通信系统、自动控制系统、导航系统等；也可以是非物理系统，例如生产管理、司法等社会经济与管理方面的系统。系统分析的主要任务是分析系统对指定激励所产生的影响。其分析过程主要包括建立系统模型，根据模型建立系统的方程，求解出系统的响应，必要时对解得的结果给出物理解释。系统分析是系统综合与系统诊断的基础。2、周期信号合成与分解 信号的频谱分析是信号处理的基础，该方法被广泛应用于电力系统、工 业 控制、电信 等领域在信号与系统课程中，从周期信号的傅里叶级数入手，引入傅里叶分析，建立信号频谱的概念，给初学者打开了另一个认识“世界”的窗口。信号的分解与合成是信号与线性系统课程中最重要的内容之一。任何满足狄里赫利条件的周期信号都是由各种不同频率、幅度和初相的正弦波叠加而成的。对周期信号由它的傅立叶级数展开式可知，各次谐波为基波频率的整数倍。而非周期信号包括了从零到无穷大的所有频率成分，每一个频率成分的幅度均趋向无穷小，但其相对大小式不同的。3、 MATLAB仿真软件进行数字信号处理的优势信号的分解与合成是本课程的重点内容之一，信号分解的方式一般有 3种：(1)用数学方式推导；(2)用实验的方法实现；(3)片J软件仿真实现比较而言，数学推导步骤繁琐，难 以形象直观；硬件实验由于实验仪器本身的局限性对实验现象和实验结论的得出都会有一定的影响，且 实际实验教学中周期信号的分解一般只能观测直流分量和前几次谐波分量，同时对周期信号的频率亦 有一定的限制，在测量中由于波形和数据较复杂，学生眼睛易疲劳，难以进行误差分析；使用MATLAB 软件仿真可以克服这些缺点，形象直观的显示信号分解与合成的过程，定量分析其中的误差程度 由The Maths Works开发的MATLAB软件是集数值分析、矩阵运算、信号处理和图形显示于一 体的可视化软件，本文应用MATI AB强大的数值分析和图形功能分析信号分解与合成。其主要功能特点如下：（1）运算功能强大。MATLAB是以矩阵为基本编程编程元素的程序设计语言，它的数值运算要素不是单个数据，而是矩阵，每个变量代表一个矩阵，矩阵由m乘n个元素，每个元素都可看做复数，所有的运算包括加、减、乘、除、函数运算等都对矩阵和复数有效；另外，通过MATLAB的符号工具箱，可以解决在数学、应用科学和工程计算领域中常常遇到的符号计算问题，强大的运算功能使其成为世界顶尖的数学应用软件之一。（2）编程效率高。MATLAB的语言规则与笔算式相似，矩阵的行数无需定义，MATLAB的命令表达方式与标准的数学表达式非常相近，因此，易写易读并易于在科技人员之间交流。MATLAB是以解释方式工作的，即它对每条语句解释后立即执行，键入算式无需编译立即得出结果，若有错误也立即做出反应，便于编程者立即更正。这些都大大减少了编程和调试的工作量，提高了编程的效率。（3）强大而智能化的作图功能。MATLAB可以方便地用图形显示二维或三维数组，将工程计算的结果可视化，使数据间的内在联系明晰了。MATLAB能智能化地根据输入的数据自动确定最佳坐标，可规定多种坐标系（如极坐标系、对数坐标系），可设置不同颜色、线型、视角等。（4）可扩展性强。MATLAB由一套程序扩展系统和工具箱，具有良好的可扩展性。工具箱是MATLAB函数的子程序库，每个工具箱都有为某个学科领域的应用而制定的，MATLAB每年都会增加一些新的工具箱。2、 周期信号的分解与合成的基本原理:1、 正交函数的基本概念：（1）如有定义在（）区间两个函数和，若满足 式 1-1 则称和在区间（）内正交。（2） 如有个函数构成一个函数集，当这些函数在区间 （）内满足 式1-2（3）如果在正交函数集之外，不存在函数 满足等式 式1-3 则此函数集称为完备正交函数集。也就是说，如能找到一个函数，使得式1-2成立，即与函数集的每一个函数都正交，那么它本身就应属于此函数集。三角函数集在区间组成正交函数集，而且是完备的正交函数集。这是因为 即三角函数集满足正交特性式（2.1-1），因而是正交函数集。2、 连续周期周期信号的分解与合成：周期信号当满足狄里赫利（Dirichlet）条件时可以展开成傅里叶级数。傅里 叶级数分三角形式和指数形式两种。狄里赫利条件如下：（1）在一个周期内，是绝对可积，即。（2）在一个周期内，的最大值和最小值的数目是有限个。（3）在一个周期内，只有有限个间断点，而且在这些间断点上，函数值必须是有限个。周期为T的信号的三角形式的傅里叶级数表示方法如下：设有周期信号，它的周期为T，角频率，则的三角傅里叶级数表示（1）一般形式为： 式2-1其中（2） 表达式二： 式2-2 式（2-2）中的系数、称为傅里叶系数，为在函数中的分量（相对大小）；为在函数中的分量,可通过下面公式计算得到： 式2-3 （3）表达式三：用余弦函数表示 式2-4 其中 初相 式2-5式2-4表明，任何满足狄里赫利条件的周期函数都可以分解为直流和许多余弦（或正弦）分量。其中第一项是常数项，它是周期信号中所包含的直流分量；式中第二项称为基波或一次谐波，它的角频率和原周期信号相同，是基波振幅，是基波初相角；式中第三项称为二次谐波，它的频率是基波频率的二倍，是二次谐波的振幅，是其初相角。以此类推，还有四三次、四次、谐波。一般而言，称为次谐波，是次谐波的振幅，是其初相角。式2-4表明，周期信号可以分解为各次谐波分量。试验及结果部分：3、 通过MATLAB实现周期信号的分解与合成：1、方波信号1.1 方波信号的分解： 下图为一个周期的方波原始信号，其周期为2S，振幅为1： 根据课题要求，由式2-3、2-5，可求得周期方波信号的的余弦分解信号的振幅、频率、相位，由于选取的谐波次数K不同，求得的结果自然也各不相同，课题要求的不同情况下的余弦分解信号的频谱图如下。（1）谐波次数K=5时：（2）谐波次数K=10：（3）谐波次数K=15： 、1.2 方波周期信号的合成根据公式2-5可求得合成方波信号，由下图知，当选取的谐波次数越高，合成得到的波形越来越接近于理想方波波形。 （1） 谐波次数K=5(2)谐波次数K=10(2)谐波次数K=15 3、 三角波信号：3.1三角波信号的分解下图为三角波原始信号，其周期为2S，振幅为1： 根据课题要求，由式2-3、2-5，可求得周期三角波信号的的余弦分解信号的振幅、频率、相位，由于选取的谐波次数K不同，求得的结果自然也各不相同，课题要求的不同情况下的余弦分解信号的频谱图如下。（1）谐波次数K=5时： （2）谐波次数K=5时：（3）谐波次数K=15时： 3.2三角波信号的合成：根据公式2-5可求得合成方波信号，由下图知，当选取的谐波次数越高，合成得到的波形越来越接近于理想方波波形。(1)谐波次数K=5 (2)谐波次数K=10(3)谐波次数K=154、 锯齿波信号4.1锯齿波信号分解下图为一个周期的锯齿原始信号，其周期为2S，振幅为1：根据课题要求，由式2-3、2-5，可求得周期锯齿波信号的的余弦分解信号的振幅、频率、相位，由于选取的谐波次数K不同，求得的结果自然也各不相同，课题要求的不同情况下的余弦分解信号的频谱图如下。（1）谐波次数K=5时： （2） 谐波次数K=10时： （3） 谐波次数K=15时：4.2锯齿波信号的合成：根据公式2-5可求得合成方波信号，由下图知，当选取的谐波次数越高，合成得到的波形越来越接近于理想锯齿波波形。（1)谐波次数K=5（2)谐波次数K=10：（2)谐波次数K=15：5、观察分析吉布斯现象： 从合成信号的波形中，可以看到，波形平坦部分中起伏的次数是由部分和所包围的谐波次数决定的，另外注意到在方波阶跃变点附近，某些点的函数值大于1或小于-1，形成过冲现象，这种在间断点处的过冲现象，称为吉布斯现象。过冲值为合成波形的最高点超过原方波信号值之比。应用MATLAB软件编写程序计算上冲或下冲的过冲值。吉布斯现象是当用信号的谐波分量的和来表述具有间断点的波形时出现的，从上述三种波形的合成过程中可以看到，随着谐波次数n取值的增大，吉布斯现象越来越显著.4、 总结1、周期信号的频谱特点： （1）离散性 周期信号的频谱非连续谱，而是离散谱。 （2）谐波性 由于n只取整数，因而谱线在频谱轴上的位置是基频的整数倍。 （3） 收敛性 幅度谱中各谱线的高度尽管不一定岁随谐波次数的增高作单调的减小 中间有可能有起伏，但总的趋势是随n的增高而减小的，当n为时，高度趋于零。 2、 当选取傅里叶有限级数的相数越多，在所合成的的波形飞f(t)中出现的峰值愈靠近f(t)的不连续点，当所选取项数N很大时，该峰值趋于一个常数，并从不连续点开始以起伏振荡的形式逐渐 衰减下去，跳变峰面积减小，通过MATLAB是这种吉布斯现象得到清楚的表现。 3、 信号中频率较低的谐波幅度的幅值大，信号波形中所含有的频率较高的告辞谐波的幅度较小，信号波形中所含有的高次谐波越多，波形的边缘越陡峭，随着N值的增大，吉布斯现象越明显。 参考文献1 张德丰 MATLAB数字信号处理与应用. M 北京 清华大学出版社 2010.12 吴大正.信号与线性系统分析（第四版）M.北京：高等教育出版社，2005: 421478 3 潘双来，邢丽冬.信号与线性系统M.北京：清华大学出版社，2006：55644 陈怀琛,吴大正,高西全.MATLAB及在电子信息课程中的应用M.北京:电子工业出版社,2002:124-162.5游春霞.徐州师范大学学报J. 学术期刊, 2006年第24卷(第01期): 2132176刘保柱，苏彦华，张宏林. MATLAB 7.0从入门到精通(修订版) M.北京：人民邮电出版社，2010：2432463 胡晓冬，董辰辉.MATLAB从入门到精通M.北京：人民邮电出版社 ,2010：78242附录：方波的分解：clear all;close all;N=200; dt=2/N; %数据点数和采样间隔for n=1:N %得到方波信号if(n*dt=1) %半个周期幅值分别为 1、-1x(n)=1;else x(n)=-1;end endfigure(1);plot(1:N)*dt,x); axis(0 2 -1.5 1.5);hold on;plot(1:N)*dt,zeros(1,N),r);xlabel(时间/s);ylabel(振幅);title(原始信号);a=zeros(1,N);m=input(输入谐波的最高次数);for k=0:m-1a(k+1)=0;b(k+1)=0;for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b=zeros(1,N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);endAn(k+1)=sqrt(a(k+1).2+b(k+1).2);ph(k+1)=(atan(b(k+1)/a(k+1)+2*pi*k*i/N);end subplot(211); stem(0:m-1)/(N*dt),An); xlabel(频率(HZ);ylabel(振幅(An);title(Fourier变换); subplot(212); stem(0:m-1)/(N*dt),ph);title(Fourier Fourier变换);xlabel(频率(HZ);ylabel(相位（rad）); 方波信号的合成：clear all;close all;N=200; dt=2/N; %数据点数和采样间隔for n=1:N %得到方波信号if(n*dt=2) x(n)=1;else x(n)=-1;end endfigure(1);subplot(211); plot(1:N)*dt,x);hold on;plot(1:N)*dt,zeros(1,N),r);xlabel(时间/s);title(振幅);a=zeros(1,N);b=zeros(1,N);mm=floor(N/2)+1;for k=0:mm-1 a(k+1)=0;b(k+1)=0;for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);endc(k+1)=sqrt(a(k+1).2+b(k+1).2);endsubplot(212);plot(0:mm-1)/(N*dt),c);title(Fourier变换);xlabel(频率/HZ);ylabel(振幅);%绘制振幅谱axis(0 25 0 1);m=input(输入谐波最大阶数);%输入最大K值if(mfloor(N/2)+1)error(谐波次数必须小于频率对应的阶数);endif(mod(N,2)=1)a(mm)=a(mm)/2;end %此时b(nn)为零a(nn)减半for i=0:N-1 %合成信号xx(i+1)=a(1)/2;for k=1:mxx(i+1)=xx(i+1)+a(k+1)*cos(2*pi*k*i/N)+b(k+1)*sin(2*pi*k*i/N);endendfigure(2);plot(1:N)*dt,xx,(0:N-1)*dt,x);%绘制合成信号与方波信号图便于比较hold on;plot(1:N)*dt,zeros(1,N),r);xlabel(绘制横轴);ylabel(振幅An);title(信号合成图);三角波信号分解clear all;close all;clear all;close all;N=200; dt=2/N; %数据点数和采样间隔for n=1:N %得到三角波信号if(n*dt=1) x(n)=-n*dt+2;else x(n)=n*dt;end endfigure(1);%打开第一个窗口，绘制原始信号及Fourier分析系数plot(1:N)*dt,x); axis(0 2 0 1);xlabel(时间/s);ylabel(振幅);title(原始信号);a=zeros(1,N);b=zeros(1,N);m=input(输入谐波的最高次数);for k=0:m-1a(k+1)=0;b(k+1)=0;for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b=zeros(1,N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);endAn(k+1)=sqrt(a(k+1).2+b(k+1).2);ph(k+1)=(atan(b(k+1)/a(k+1)+2*pi*k*i/N);end subplot(211); stem(0:m-1)/(N*dt),An); xlabel(频率(HZ);ylabel(振幅(An);title(Fourier变换谐波次数K=10);subplot(212); stem(0:m-1)/(N*dt),ph);title(Fourier变换谐波次数K=,num2str(m);xlabel(频率(HZ);ylabel(相位（rad）); 三角波信号合成：clear all;close all;N=200;dt=2/N; %数据点数和采样间隔for n=1:N %得到三角波信号t=mod(n*dt,2);if(t1) x(n)=-n*dt+2;else x(n)=n*dt;end endfigure(1);%打开第一个窗口，绘制原始信号及Fourier分析系数subplot(211); plot(1:N)*dt,x);xlabel(时间/s);ylabel(振幅);axis(0 2 0 1.2)title(原始信号);a=zeros(1,N);b=zeros(1,N);mm=floor(N/2)+1;Fork=0:mm-1 a(k+1)=0;b(k+1)=0;for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N); %求解fourier 系数b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);endc(k+1)=sqrt(a(k+1).2+b(k+1).2);endsubplot(212);plot(0:mm-1)/(N*dt),c);title(Fourier变换);xlabel(频率/HZ);ylabel(振幅);%绘制振幅谱axis(0 25 0 1);m=input(输入谐波最大阶数);%输入最大k值if(mfloor(N/2)+1)%输入的最大K值如大于Nyquist 频率对应的点数，则显示出错信息error(谐波阶数必须小于频率对应的阶数);endif(mod(N,2)=1)a(mm)=a(mm)/2;end %此时b(nn)为零，a（nn）减半for i=0:N-1 %合成信号xx(i+1)=a(1)/2;for k=1:mxx(i+1)=xx(i+1)+a(k+1)*cos(2*pi*k*i/N)+b(k+1)*sin(2*pi*k*i/N);endendfigure(2);%打开第二个图形窗口，绘制合成信号图plot(1:N)*dt,xx,(0:N-1)*dt,x);%绘制合成信号和方波信号图便于比较axis(0 2 0 1.2);xlabel(绘出横轴 时间/s);ylabel(振幅An);title(信号合成图谐波次数K=,num2str(m); 锯齿波信号分解： clear all;close all;N=200; dt=2/N; %数据点数和采样间隔for n=1:N %得到锯齿角波信号if(n*dt2) x(n)=atan(1/2)*n*dt;else x(n)=0;end endfigure(1);%打开第一个窗口，绘制原始信号及Fourier分析系数plot(1:N)*dt,x); axis(