用MATLAB实现数字信号处理.docx
《用MATLAB实现数字信号处理.docx》由会员分享，可在线阅读，更多相关《用MATLAB实现数字信号处理.docx（12页珍藏版）》请在冰豆网上搜索。
用MATLAB实现数字信号处理
摘要
MATLAB（MatrixLaboratory-----矩阵实验室）是由美国Mathworks公司于1984年正式推出的一种以矩阵运算为基础的交互式程序语言。
与其他的语言相比她有着“简洁和智能化，适应科技专业人员的思维方式和书写习惯”的优点，所以为科技人员所乐于接受。
更因为它可适用于各种平台，并且随计算机硬、软件的更新而及时升级。
所以MATLAB逐渐成为当前我国大学教学和科学研究中一种非常重要的工具软件。
首先通过的MATLAB和DSP（DigitalSignalProcessing）------数字信号处理的简单概括，在第二章里利用MATLAB中几个简单而又强大的命令（比如：
plot,函数的赋值，M函数等）来对数字信号处理中信号的取样和还原进行模拟。
本文先从理论上说明：
原函数与冲激函数相乘得到的函数fs（t）其频谱是原函数频谱的无限个频移项组成，其幅值是原函数的1/Ts.而且，只有当取样频率fs大于原函数fm的两倍时，那么相后的频谱就不会混叠了，也就是说，取样函数（fs=f（t）*s（t））包含了原函数的全部信息。
那么我们就可以在接收端用一个低通滤波器（其频率响应的幅度是Ts）,就可以把原信号还原出来了。
反之，如果采样的时候，采样频率小于2*fm,那么采样后得到的函数的频谱会出现混叠。
这样，采样函数就不能包括原函数的全部信息，使得接收断的低通滤波器所截取的信号不能真实的反映原函数。
关键词：
MATLAB,DSP,取样频率，信号最高频率
Abstract
MATLAB（MatrixLaboratory）wasascientistcomputerlanguagewhichwasinventedbytheAmericancompanyMathworksin1984.ThislanguageisbasedonMatrixprocess.Comparedwithothercomputersciencelanguages,ithastheadvantagesinitsclearanceintheexpressionofitslanguageandfitsthethinkingmethodofmanyscientists.So,asthedaysgoesby,itbecomesanacceptableandenjoyablelanguageforthescientists.Forthemostimportantaspects,itcouldbeinstalledinmanydifferentoperatesystemsandbeupdateditsversionaccordingtothedevelopmentofcomputersoftwareandhardware,thus,MATLABbecomesanimportantteachingtoolsinmanyChinauniversity.
Inthefirstchapter,thisarticleprovidesusanoveralldescriptionoftheMATLABandDSP（DigitalSignalProcessing）.Then,inthesecondchapter,thisarticleusesomesimplebutpowerfulcommandstosimulatetheexperimentofthesamplingandreconstruction.Thereaderwillfindinthesecondchapter,thearticlefirstprovidesthetheorywhythesignalcouldbereconstructed:
whenthesamplingfrequency（fs）islargerorequalthanthetwotimesofthemaximumfrequencyofthesignal（fm）.Then,thesamplingsignalcouldincludethewholeinformationoftheoriginalsignalandwhenweusealowpassfilterintheend,wecouldreconstructiontheoriginalsignal.Otherwise,ifthefs<=2fm,thesamplingsignalcouldnotprovideusthewholeinformationoftheoriginalsignal,thusbasedonthiswecouldnotreconstructtheoriginalsignal.
Keyword:
MATLAB,DSP,samplingfrequency，themaximumfrequencyofthesignal
第一章绪论
第一节MATLAB语言简介
1.1.1MATLAB语言概述
MATLAB（MatrixLaboratory-----矩阵实验室）是一种科学计算软件，适用于工程应用各领域的分析设计与复杂计算，它使用方便，输入简捷运算高效且内容丰富，很容易由用户自行扩展。
因此，当前以成为美国和其他发达国家大学教学和科学研究中最常用且必不可少的工具。
MATLAB是由美国Mathworks公司于1984年正式推出的，到1988年有了3.x（DOS）版本；1992年推出4.x（Windows）版本；1997年推出5.1（Windows）版本。
以后又升级到5.3（也称R11）版本。
2000年下半年，Mathworks公司推出了他们的最新产品MATLAB6.0（R12）试用版，并于2001年初推出了正式版。
随着版本的升级，内容不断扩充，功能更加强大。
另一方面对使用环境也提出了更高的要求。
近几年来，Mathworks公司在将MATLAB语言运用于系统仿真和实时运行方面，取得了很多成绩，更扩大了它的应用前景。
MATLAB是“矩阵实验室”（MATrixLABoratoy）的缩写，它是一种以矩阵运算为基础的交互式程序语言，着重针对科学运算，工程计算和绘图的需求。
与其他计算机语言相比，其特点是简洁和智能化，适应科技专业人员的思维方式和书写习惯，使得编程和调试效率大大提高。
它用解释方式工作，键入程序立即得出结果，人机交互性能好，为科技人员所乐于接受。
特别是它可适用于各种平台，并且随计算机硬、软件的更新而及时升级。
因此，MATLAB语言在国外的大学工学院中，特别是数值计算用的最频繁的电子信息类学科中，已成为每个学生都应掌握的工具了。
它大大提高了课程教学、解题作业、分析研究的效率。
1.1.2MATLAB语言的特点
MATLAB语言有以下特点：
1．起点高
（1）每个变量代表一个矩阵。
当前的科学计算中，几乎无处不用矩阵运算，这使它的优势得到充分的体现。
在MATLAB中，每个变量代表一个矩阵，它可以有n*m个元素。
（2）每个元素都看作复数，这个特点在其他语言中是不多见的。
（3）所有的运算都对矩阵和复数有效，包括加、减、乘、除、函数运算等。
2．人机界面适合科技人员。
（1）MATLAB的语言规则和笔算式相似。
（2）矩阵的行列无需定义。
要输入一个矩阵，用其他语言时必须先定义矩阵的阶数，而MATLAB则不必有阶数定义语句。
（3）键入算式立即得结果，无需编译。
MATLAB是以解释方式工作的，即它对每条语句解释后立即执行，若有错误也立即做出反应，便于编程者立即改正。
3.强大而简易的做图功能。
（1）能根据输入数据自动确定坐标绘图。
（2）能规定多种坐标系。
（3）能绘制三维坐标中的曲线和曲面。
（4）可设置不同颜色、线型、视角等。
4.智能化程度高。
（1）绘图时自动选择最佳坐标，大大方便了用户。
（2）做数值积分时自动按精度选择步长。
（3）自动检测和显示程序错误的能力强，易于调试。
5.功能丰富，可扩展性能强。
扩展部分称为工具箱。
它实际上是用MATLAB的基本语句编成的各种子程序集，用于解决某一方面的专门问题，或实现某一类的新算法。
现在已经有控制系统、信号处理、图象处理、系统辩识、模糊集合、神经元网络及小波分析等工具箱。
第二节数字信号处理简介
1.2.1数字信号处理概述
数字信号处理，或者说对信号的数字处理，是60年代前后发展起来的一门新兴学科，进入70年代以来，随着电子计算机、大规模集成电路（LSI）和超大规模集成电路（VLSI），以及微处理器技术的迅速发展、数字信号处理无论在理论上还是在工程应用中，都是目前发展最快的学科之一，而且日趋完善和成熟。
目前对数字信号处理没有一个严格的定义，我们不妨把用包括软件在内的数字技术来处理模拟信号叫做数字信号处理。
或者简单的说，数字信号处理就是用数字的方法对信号进行变换和处理。
从学科的内容来看，数字滤波器和快速傅立叶变换是数字信号处理的两个最重要的学科分支，也是数字信号处理的核心内容。
事实上，对信号进行调制、滤波、变换、增强、压缩、检测、估计、识别等，都是数字信号处理的研究范围。
用数字的方法实现具有各种功能的信号处理处理系统，可以使系统做到更加经济和小型化，进而有可能实现过去用模拟技术难以实现的各种功能。
图1给出了了模拟信号处理数字化处理系统的简化框图。
此系统先将模拟信号变换为数字信号，经数字信号处理后，在变换成模拟信号输出，其中抗混叠滤波器的作用是将输入信号x（t）中高于折叠频率（其值等于采样频率的一半）的分量滤除，以防信号频谱的混叠。
随后，信号经采样和A/D变换后，变成数字信号x（n）,数字信号处理器对x（n）进行处理，得到输出数字信号，经D/A变换器变成模拟信号。
此信号经低通滤波器，滤除不需要的高频分量，然后输出平滑的模拟信号y（t）.
都是信号处理、ryy_____________________________________________________________________________________________________________
1.2.2信号是如何被处理的？
我们在实际中所遇到的信号大部分是模拟信号。
它们的幅度随时间连续地变化。
通常用由无源或者有源器件所组成的电网络对他们作处理。
这中方法称为模拟信号处理。
例如无线电和电视接收机。
模拟信号：
x（t）y（t）
它们也可以用包括加法器、乘法器和逻辑元件或用微处理器组成的数字硬件来处理。
但是必须先把模拟信号变成适合于数字硬件来处理的形式。
这样的信号形式成为数字信号。
它在一定的时间间隔内只取有限数目的数值之一，因而它能用二进制或位来表示。
数字信号的处理称为DSP,它可以用框图来表示。
图中各个方框的意义介绍如下：
PrF:
这是一个前置滤波器或反泄漏滤波器，它限制了模拟信号以防止泄漏。
ADC:
模-数转换器，它是由模拟信号产生出一个二进制流。
数字信号处理器：
这是DSP的心脏，它可以代表一个通用计算机或专用处理器，或数字硬件等等。
DAC:
这是ADC的逆运算，称为数-模转换器，它由一个二进制流产生一个阶梯波形，作为形成模拟信号的第一步。
PoF:
这是一个后置滤波器，它把阶梯波形平滑成为预期的模拟信号。
1.2.3DSP胜于ASP的优点
ASP的主要缺点是在做复杂信号处理时只有有限的能力。
这造成了处理的不灵活性和系统时间的复杂性。
这也造成了产品的昂贵。
另一方面，采用DSP的方法，它可以把一个便宜的个人计算机变成一个强大的信号处理器。
DSP的一些重要的优点如下：
（1）用DSP方法的系统的开发可以用通用计算机上的软件来进行。
因此DSP比较容易开发和测试，而软件是可移植的。
（2）DSP运算是单纯的基于加法和乘法，这导致特别稳定的处理性能——例如不受温度影响。
（3）DSP运算可以方便地作实时的修改，通常只要改变程序，或者对寄存器重新加载。
（4）由于超大规模集成电路（VLSI）降低了寄存器、门、微处理器的价格，DSP比较便宜。
DSP的主要缺点是运算速度低，特别是对高频的信号。
主要是由于它的上述有限，DSP现在已经是很多技术的首选方法，比如消费电子、通信、无线电话和医学图象等。
第二章用MATLAB实现信号处理
*注：
本文中的所有程序均在matlab6.5下调试运行
信号经过采样和量化运算后变成离散信号。
这些离散信号由数字信号处理器处理后，再用一个重构运算后变成模拟信号。
用傅立叶分析可以从频域的观点来描述采样运算，分析其效果，然后重新把信号还原出来。
我们假定量化的等级足够的多因而离散信号量化误差可以忽略不计。
第一节信号的取样：
取样定理论述了在一定的条件下,一个连续时间的信号完全可以用该信号在等时间间隔上的瞬时值表示,这些样本值完全包含了该连续时间信号的全部信息,利用这些样本值可以恢复原信号.可以说,取样定理在连续时间信号和离散时间信号之间架起了一座桥梁.可以说取样定理为连续时间信号与离散时间信号的相互转换提供了理论依据.
所谓的”取样”就是利用取样脉冲序列s（t）从连续时间信号f（t）中”抽取”一系列离散样本值的过程.这样得到的离散信号称为取样信号.如果图2.1所示:
图2.1取样的模型
式中取样脉冲序列s（t）也称开关函数。
如果其各脉冲间隔的时间相同，均为Ts,就称为均匀取样。
Ts称为取样周期，fs=1/Ts称为取样频率或取样率。
Ws=2*pi*fs=2*pi/Ts.
如果令f（t）和s（t）的傅立叶变换为F（jw）和S（jw）,则由频域卷积定理，得取样信号fs（t）的频谱函数:
Fs（jw）=（1/2*pi）*F（jw）*S（jw）（式2.1）
2.1.1冲激取样：
如果取样序列s（t）是周期为Ts的冲激函数序列x（t）,则称为冲激取样。
且冲激序列x（t）的频谱函数也是周期序列，即：
F[s（t）]=F（x（t））=F[∑x（t-n*Ts）]=ws∑x（w-n*ws）其中n从-∞到+∞（式2.2）
把（式2.2）代入（式2.1）中，可以得到：
Fs（jw）=（1/2*pi）*F（jw）*ws*∑x（w-nws）=1/Ts*∑F（j（w-nws））其中n从-∞到+∞（式2.3）
由（式2.3）可以知道，取样信号fs（t）的频谱由原信号频谱F（jw）的无限个频移项组成，其频谱的角频率为nws,其幅值为原频谱的1/Ts。
2.1.2矩形脉冲取样：
如果取样脉冲序列s（t）是幅度为1，脉宽为L（LS（jw）=P（jw）=F（p（t））=（2*pi*L）/Ts*∑Sa（n*ws*L/2）x（w-nws）（式2.4）
上式代入式（式2.1）可以得到取样信号fs（t）的频谱函数为：
Fs（jw）=1/（2*pi）*F（jw）*（2*pi*L/Ts）*∑Sa（nwsL/2）*x（w-nws）=L/T*∑Sa（nws*L/2）*F[j（w-nws）]（式2.5）
比较式2.5和式2.1可以发现，经过冲激取样或矩形脉冲取样后，其取样信号fs（t）的频谱很相似。
因此当ws>2wm时取样信号的频谱Fs（jw）也不会出现混叠，从而能从取样信号fs（t）中恢复信号f（t）。
可以看出，如果ws>2wm,那么各相邻频谱不会发生互相重叠，这时候就设法（比如用低通滤波器），从取样信号的频谱Fs（jw）中得到原信号的频谱，即从取样信号fs（t）中恢复原信号f（t）,如果ws<2wm，那么频移后的各相邻频谱将互相重叠，如图所示。
这样就无法将它们分开，因而也就不能再恢复原信号。
频谱重叠的这种现象叫做混叠现象。
可见，为了不发生混叠现象，ws>2wm。
第二节信号的重构
从采样定理和上述的例子可以清楚的看到，如果对有限带宽信号xa（t）以高于2fc的频率进行采样，就可以从其采样序列x（n）重构原模拟信号。
重构可以考虑为两个步骤：
1.先把样本集转换为一个加权脉冲串列。
2.然后把这个脉冲串列通过一个带宽为[-Fs/2,Fs/2]的低通滤波器进行滤波。
现在以冲激取样为例，研究如何从取样信号fs（t）恢复信号f（t）。
设有冲激信号fs（t）,其取样角频率ws>2wm（wm为原信号的最高频率）。
Fs（jw）为fs（t）的频谱函数，为了从Fs（jw）中无失真地恢复F（jw）,选择一个理想的低通滤波器，其频率响应的幅度为Ts,截止角频率为wc（wmTs,|w|H（jw）={0,|w|>wc
因此：
F（jw）=Fs（jw）*H（jw）（式2.6）
即恢复了原信号的频谱函数F（jw）。
根据卷积定理，（式2.6）的相应于时域为：
f（t）=fs（t）*h（t）（式2.7）
由于冲激取样信号
fs（t）=f（t）·s（t）=f（t）·∑x（t-n*Ts）=∑f（n*Ts）x（t-nTs）（式2.8）其中n从-∞到+∞
为简便，选wc=ws/2,则Ts=2*pi/ws=pi/wc，得：
h（t）=Sa（ws*t/2）（式2.9）
把（式2.8）和（式2.9）代入（式2.7）中可以得到：
f（t）=∑f（nTs）x（t-nTs）*Sa（wst/2）=∑f（nTs）Sa[ws/2*（t-nTs）]=∑f（nTs）Sa（ws*t/2-n*pi）
上式表明，连续信号f（t）可以展开成正交取样函数（Sa函数）的无穷级数，该级数的系数等于取样值f（nTs）。
也就是说，若在取样信号fs（t）的每个样点出，画一个最大峰值为f（nTs）的波形函数，那么其合成波形就是原信号f（t）.因此，只要以知各取样值f（nTs）,就能唯一的确定出原信号f（t）。
2.2.1时域取样定理
一个频谱在区间（-wm,wm）以外为零的频带有限信号f（t）,可以由其样点值f（nTs）确定。
需要注意的是，为了能从取样信号fs（t）中恢复原信号f（t）,需满足两个条件：
（1）f（t）必须是带限信号，其频谱函数|w|>wm各处为零；
（2）取样频率不能过低，必须fs>2fm,或者说取样间隔不能太长，必须Ts<1/（2fm）,否则将发生混叠。
第三节MATLAB对信号重建和还原的实现
2.3.1离散傅立叶变换
因为信号的采样中需要用到离散傅立叶变换，故这里单独拿出来，并且将它编成一个子程序。
严格的说MATLAB不能直接处理无限长连续的信号，但是根据离散傅立叶变换的两个性质：
周期性和对称性，我们时间间隔dt或者频率间隔dw取的很小，然后来估计表达式在[0,pi]上，或者一个时间段里的特性。
假定x（n）是有限长的，我们在[0，pi]间按等间隔频点来估计X（exp（jw））,则可以用矩阵—向量相乘的运算来实现。
假设
x（n）在[n1,nN]之间，且wk=（pi/M）*k,k=0,1,……M
那么X（exp（jwk））=∑exp（-j*（pi/M）*k*nl）x（nl）,其中l范围为[1，N]
在MATLAB中，我们把序列和下标排成行向量，得到：
X=x*exp（-j*w*n’）
编写一个离散傅立叶变换子程序dtft.m:
functionX=dtft（x,w）
X=x*exp（-j*[1:
length（x）’*w]）;
下面我们以函数x（t）=A*exp（-a*t）*sin（b*t）（其中a=b=50*sqrt
（2）*pi）为例来模拟信号的取样和还原，并对上面的理论进行验证。
2.3.2采样
同样，因为MATLAB不能直接处理无限长连续的信号，我们把dt取的很小，dt<xc（m）=xa（m*dt）所以
Xa（jΩ）=∑xc（m）e（-j*Ωmdt）dt，在这个程序中，分别设定四个采样频率：
fs=10kHz,1kHz,400Hz,200Hz.
程序如下：
clearall;
closeall;
fs=10000;fs1=1000;fs2=400;fs3=200;
t=0:
1/fs:
0.1;
A=444.128;a=50*sqrt
（2）*pi;b=a;
xa=A*exp（-a.*t）.*sin（b.*t）;
k=0:
511;f=fs*k/512;
Xa=dtft（xa,2*pi*k/512）;
T1=1/fs1;t1=0:
T1:
0.1;
x1=A*exp（-a.*t1）.*sin（b.*t1）;
X1=dtft（x1,2*pi*k/512）;
T2=1/fs2;t2=0:
T2:
0.1;
x2=A*exp（-a.*t2）.*sin（b.*t2）;
X2=dtft（x2,2*pi*k/512）;
T3=1/fs3;t3=0:
T3:
0.1;
x3=A*exp（-a.*t3）.*sin（b.*t3）;
X3=dtft（x3,2*pi*k/512）;
figure
（1）;
subplot（4,2,1）;plot（t,xa）
subplot（4,2,3）;stem（t1,x1,'.'）;
subplot（4,2,5）;stem（t2,x2,'.'）;
subplot（4,2,7）;stem（t3,x3,'.'）;
subplot（4,2,2）;plot（f/pi,abs（Xa）/max（abs（Xa）））;
subplot（4,2,4）;plot（f/pi,abs（X1）/max（abs（X1）））;
subplot（4,2,6）;plot（f/pi,abs（X2）/max（abs（X2）））;
subplot（4,2,8）;plot（f/pi,abs（X3）/max（abs（X3）））;
结果为
上图中，自上到下的四个图分别为模拟信号，采样频率分别为1kHz,400Hz,200Hz的图以及他们的幅度谱。
可以看出，fs=1kHz的失真较小，fs=400Hz的失真较大，而fs=200Hz失真比较严重。
2.3.3重建
根据前面的公式xa（t）=∑x（n）sinc[Fs（t-nTs）];设定一个ti的值求xa（ti）的问题，可归结为一个行向量x（n）和一个同长的由n’构成的列向量g（ti-n’T）相乘。
这与前面的dtft的思想比较类似。
“t-n*T”这个的实现可以把t设成行向量，nT’为列向量。
我们的目的是把它构成一个行数与n一样，列数与t一样的矩阵。
因此可以把t右乘ones（length（n）,1）;把nT左乘向量ones（1,length（t））就可以了。
可以令TNT=ones（length（n）,1）t-n’*T*ones（1,length（t））
程序如下：
clearall;
close
升级会员 