现代数字信号处理实验报告
数字信号处理实验报告
实验一 信号、系统及系统响应一、实验目的1、熟悉理想采样的性质，了解信号采样前后的频谱变化，加深对时域采样定理的理解。
2、熟悉离散信号和系统的时域特性。
3、熟悉线性卷积的计算编程方法：利用卷积的方法，观察、分析系统响应的时域特性。
4、掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对离散信号、系统及其系统响应进行频域分析。
二、 实验原理1.理想采样序列：对信号x a (t)=A e −αt sin(Ω0t )u(t)进行理想采样，可以得到一个理想的采样信号序列x a (t)=A e −αt sin(Ω0nT ),0≤n ≤50，其中A 为幅度因子，α是衰减因子，Ω0是频率，T 是采样周期。
2.对一个连续时间信号x a (t)进行理想采样可以表示为该信号与一个周期冲激脉冲的乘积，即x ̂a (t)= x a (t)M(t),其中x ̂a (t)是连续信号x a (t)的理想采样；M(t)是周期冲激M(t)=∑δ+∞−∞(t-nT)=1T ∑e jm Ωs t +∞−∞,其中T 为采样周期，Ωs =2π/T 是采样角频率。
信号理想采样的傅里叶变换为X ̂a (j Ω)=1T ∑X a +∞−∞[j(Ω−k Ωs )],由此式可知：信号理想采样后的频谱是原信号频谱的周期延拓，其延拓周期为Ωs =2π/T 。
根据时域采样定理，如果原信号是带限信号，且采样频率高于原信号最高频率分量的2倍，则采样以后不会发生频率混叠现象。
三、简明步骤产生理想采样信号序列x a (n),使A=444.128，α=50√2π，Ω0=50√2π。
（1） 首先选用采样频率为1000HZ ，T=1/1000，观察所得理想采样信号的幅频特性，在折叠频率以内和给定的理想幅频特性无明显差异，并做记录；（2） 改变采样频率为300HZ ，T=1/300，观察所得到的频谱特性曲线的变化，并做记录；（3） 进一步减小采样频率为200HZ ，T=1/200，观察频谱混淆现象是否明显存在，说明原因，并记录这时候的幅频特性曲线。
现代数字信号处理实验报告
现代数字信号处理实验报告1、估计随机信号的样本自相关序列。
先以白噪声()x n 为例。
(a) 产生零均值单位方差高斯白噪声的1000个样点。
(b) 用公式：9991ˆ()()()1000x n r k x n x n k ==-∑估计()x n 的前100个自相关序列值。
与真实的自相关序列()()x r k k δ=相比较，讨论你的估计的精确性。
(c) 将样本数据分成10段，每段100个样点，将所有子段的样本自相关的平均值作为()x n 自相关的估值，即：999001ˆ()(100)(100) , 0,1,...,991000x m n r k x n m x n k m k ===+-+=∑∑与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列()()x r k k δ=吗？(d) 再将1000点的白噪声()x n 通过滤波器11()10.9H z z-=-产生1000点的y (n )，试重复(b)的工作，估计y (n )的前100个自相关序列值，并与真实的自相关序列()y r k 相比较，讨论你的估计的精确性。
仿真结果：(a)图1.1 零均值单位方差高斯白噪声的1000个样本点分析图1.1：这1000个样本点是均值近似为0，方差为1的高斯白噪声。
(b)图1.2 ()x n的前100个自相关序列值分析上图可知：当k=0时取得峰值，且峰值大小比较接近于1，而当k ≠0时估计的自相关值在0附近有小幅度的波动，这与真实自相关序列r x (k)=δ(k)比较接近，k ≠0时估计值非常接近0，说明了估计的结果是比较精确的。
(c)图1.3基于Bartlett 法的前100个自相关序列值与(b)的结果相比，同样在k=0时达到峰值，k ≠0时0值附近上下波动；估计值的方差比较小，随着k 的增大波动幅度逐渐变小，在k 较大时它更接近真实自相关序列()()x r k k δ=。
即采用分段方法得到的自相关序列的估计值更加接近r x (k)=δ(k)。
数字信号处理实验报告(西电)
数字信号处理实验报告班级：****姓名：郭**学号：*****联系方式：*****西安电子科技大学电子工程学院绪论数字信号处理起源于十八世纪的数学，随着信息科学和计算机技术的迅速发展，数字信号处理的理论与应用得到迅速发展，形成一门极其重要的学科。
当今数字信号处理的理论和方法已经得到长足的发展，成为数字化时代的重要支撑，其在各个学科和技术领域中的应用具有悠久的历史，已经渗透到我们生活和工作的各个方面。
数字信号处理相对于模拟信号处理具有许多优点，比如灵活性好，数字信号处理系统的性能取决于系统参数，这些参数很容易修改，并且数字系统可以分时复用，用一套数字系统可以分是处理多路信号；高精度和高稳定性，数字系统的运算字符有足够高的精度，同时数字系统不会随使用环境的变化而变化，尤其使用了超大规模集成的DSP 芯片，简化了设备，更提高了系统稳定性和可靠性；便于开发和升级，由于软件可以方便传送，复制和升级，系统的性能可以得到不断地改善；功能强，数字信号处理不仅能够完成一维信号的处理，还可以试下安多维信号的处理；便于大规模集成，数字部件具有高度的规范性，对电路参数要求不严格，容易大规模集成和生产。
数字信号处理用途广泛，对其进行一系列学习与研究也是非常必要的。
本次通过对几个典型的数字信号实例分析来进一步学习和验证数字信号理论基础。
实验一主要是产生常见的信号序列和对数字信号进行简单处理，如三点滑动平均算法、调幅广播（AM ）调制高频正弦信号和线性卷积。
实验二则是通过编程算法来了解DFT 的运算原理以及了解快速傅里叶变换FFT 的方法。
实验三是应用IRR 和FIR 滤波器对实际音频信号进行处理。
实验一●实验目的加深对序列基本知识的掌握理解●实验原理与方法1.几种常见的典型序列：0()1,00,0(){()()(),()sin()j n n n n u n x n Aex n a u n a x n A n σωωϕ+≥<====+单位阶跃序列：复指数序列：实指数序列：为实数 正弦序列：2.序列运算的应用：数字信号处理中经常需要将被加性噪声污染的信号中移除噪声，假定信号 s(n)被噪声d(n)所污染，得到了一个含噪声的信号()()()x n s n d n =+。
数字信号处理实验报告_五个实验
实验一 信号、系统及系统响应一、 实验目的1、熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解；2、熟悉时域离散系统的时域特性；3、利用卷积方法观察分析系统的时域特性；4、掌握序列傅立叶变换的计算机实现方法，利用序列的傅立叶变换对连续信号、离散信号及系统响应进行频域分析。
二、 实验原理及方法采样是连续信号数字处理的第一个关键环节。
对采样过程的研究不仅可以了解采样前后信号时域和频域特性发生变化以及信号信息不丢失的条件，而且可以加深对傅立叶变换、Z 变换和序列傅立叶变换之间关系式的理解。
对一个连续信号)(t x a 进行理想采样的过程可用下式表示：)()()(^t p t t xx aa=其中)(^t x a 为)(t x a 的理想采样，p(t)为周期脉冲，即∑∞-∞=-=m nT t t p )()(δ)(^t x a的傅立叶变换为)]([1)(^s m a m j X T j a XΩ-Ω=Ω∑∞-∞=上式表明^)(Ωj Xa为)(Ωj Xa的周期延拓。
其延拓周期为采样角频率（T /2π=Ω）。
只有满足采样定理时，才不会发生频率混叠失真。
在实验时可以用序列的傅立叶变换来计算^)(Ωj X a 。
公式如下：Tw jw ae X j X Ω==Ω|)()(^离散信号和系统在时域均可用序列来表示。
为了在实验中观察分析各种序列的频域特性，通常对)(jw e X 在[0，2π]上进行M 点采样来观察分析。
对长度为N 的有限长序列x(n)，有：n jw N n jw k ke m x eX--=∑=)()(1其中，k Mk πω2=，k=0,1,……M-1 时域离散线性非移变系统的输入/输出关系为 ∑∞-∞=-==m m n h m x n h n x n y )()()(*)()(上述卷积运算也可在频域实现)()()(ωωωj j j e H e X eY =三、 实验程序s=yesinput(Please Select The Step Of Experiment:\n 一.(1时域采样序列分析 s=str2num(s); close all;Xb=impseq(0,0,1); Ha=stepseq(1,1,10);Hb=impseq(0,0,3)+2.5*impseq(1,0,3)+2.2*impseq(2,0,3)+impseq(3,0,3); i=0;while(s);%时域采样序列分析 if(s==1) l=1; k=0;while(1)if(k==0)A=yesinput('please input the Amplitude:\n',...444.128,[100,1000]); a=yesinput('please input the Attenuation Coefficient:\n',...222.144,[100,600]); w=yesinput('please input the Angle Frequence(rad/s):\n',...222.144,[100,600]); end k=k+1;fs=yesinput('please input the sample frequence:\n',...1000,[100,1200]); Xa=FF(A,a,w,fs); i=i+1;string+['fs=',num2str(fs)]; figure(i)DFT(Xa,50,string); 1=yesinput 1=str2num(1); end%系统和响应分析else if(s==2)kk=str2num(kk);while(kk)if(kk==1)m=conv(Xb,Hb);N=5;i=i+1;figure(i)string=('hb(n)');Hs=DFT(Hb,4,string);i=i+1;figure(i)string('xb(n)');DFT(Xb,2,string);string=('y(n)=xb(n)*hb(n)');else if (kk==2)m=conv(Ha,Ha);N=19;string=('y(n)=ha(n)*(ha(n)');else if (kk==3)Xc=stepseq(1,1,5);m=conv(Xc,Ha);N=14;string=('y(n)=xc(n)*ha(n)');endendendi=i+1;figure(i)DFT(m,N,string);kk=yesinputkk=str2num(kk);end卷积定理的验证else if(s==3)A=1;a=0.5;w=2,0734;fs=1;Xal=FF(A,a,w,fs);i=i+1;figure(i)string=('The xal(n)(A=1,a=0.4,T=1)'); [Xa,w]DFT(Xal,50,string);i=i+1;figure(i)string =('hb(n)');Hs=DFT(Hb,4,string);Ys=Xs.*Hs;y=conv(Xal,Hb);N=53;i=i+1;figure(i)string=('y(n)=xa(n)*hb(n)');[yy,w]=DFT(y,N,string);i=i+1;figure(i)subplot(2,2,1)plot(w/pi,abs(yy));axis([-2 2 0 2]);xlabel('w/pi');ylabel('|Ys(jw)|');title(FT[x(n)*h(n)]');subplot(2,2,3)plot(w/pi,abs(Ys));axis([-2 2 0 2]);xlabel('w/pi');ylabel('|Ys(jw)|');title('FT[xs(n)].FT[h(n)]');endendend子函数:离散傅立叶变换及X(n),FT[x(n)]的绘图函数function[c,l]=DFT(x,N,str)n=0:N-1;k=-200:200;w=(pi/100)*k;l=w;c=x*Xc=stepseq(1,1,5);子函数：产生信号function c=FF(A,a,w,fs)n=o:50-1;c=A*exp((-a)*n/fs).*sin(w*n/fs).*stepseq(0,0,49); 子函数：产生脉冲信号function [x,n]=impseq(n0,n1,n2)n=[n1:n2];x=[(n-n0)==0];子函数：产生矩形框信号function [x,n]=stepseq(n0,n1,n2) n=[n1:n2];x=[(n-n0>=0)];四、 实验内容及步骤1、认真复习采样理论，离散信号与系统，线性卷积，序列的傅立叶变换及性质等有关内容，阅读本实验原理与方法。
数字信号处理综合实验报告
综合实验1. 实验目的能综合利用信号处理的理论和Matlab 工具实现对信号进行分析和处理(1)熟练对信号进行时域和频域分析；(2)熟练进行滤波器设计和实现；(3)掌握对信号的滤波处理和分析。
2．实验原理设计并实现滤波器对信号进行分析和处理是信号处理课程学习的主要内容。
通过对信号进行频谱分析，能发现信号的频率特性，以及组成信号的频率分量。
对信号进行滤波处理，能改善信号的质量，或者为数据处理(如传输，分类等)提供预处理，等。
本次实验是对特定信号进行分析并进行滤波处理，需要综合应用之前的实验内容，主要有以下几个方面。
(1)离散时间信号与系统的时域分析Matlab 为离散时间信号与系统的分析提供了丰富且功能强大的计算函数和绘图分析函数，便于离散时间信号和系统的时域表示和分析。
(2)信号的频域分析信号处理课程主要学习了离散信号和系统的频域分析方法与实现，以及滤波器的设计与实现。
离散信号与系统的频域分析包括DTFT DFT Z变换等，FFT则是DFT的快速实现。
用Matlab分析信号的频谱可以用freqz函数或者FFT函数。
(3)滤波器设计滤波器的设计首先要确定滤波器的类型，即低通、高通、带通还是带阻。
滤波器的边缘频率可以通过对信号的频谱分析得到，滤波器的幅度指标主要有阻带最小衰减As 和通带最大衰减Ap。
一般来说，As越大，对截止通过的频率分量的衰减越大；Ap越小，对需要保留的频率分量的衰减越小。
因此，As 越大，Ap 越小，滤波器的性能越好，但随之而来，滤波器的阶数越大，实现的代价(包括计算时间和空间)越大。
由此，滤波器的设计需要对滤波器性能和实现代价进行均衡考虑。
另外根据冲激响应的长度可以分为IIR 和FIR 两种类型。
两种类型的滤波器各有特点。
用FIR 滤波器可以设计出具有严格线性相位的滤波器，但在满足同样指标的条件下，FIR 滤波器的阶数高于IIR 滤波器。
Matlab 为各种类型的滤波器的设计提供了丰富的函数，可以借助这些函数方便地设计出符合要求地滤波器。
数字信号处理实验报告
数字信号处理实验报告引言数字信号处理（Digital Signal Processing，DSP）是一门研究数字信号的获取、分析、处理和控制的学科。
在现代科技发展中，数字信号处理在通信、图像处理、音频处理等领域起着重要的作用。
本次实验旨在通过实际操作，深入了解数字信号处理的基本原理和实践技巧。
实验一：离散时间信号的生成与显示在实验开始之前，我们首先需要了解信号的生成与显示方法。
通过数字信号处理器（Digital Signal Processor，DSP）可以轻松生成和显示各种类型的离散时间信号。
实验设置如下：1. 设置采样频率为8kHz。
2. 生成一个正弦信号：频率为1kHz，振幅为1。
3. 生成一个方波信号：频率为1kHz，振幅为1。
4. 将生成的信号通过DAC（Digital-to-Analog Converter）输出到示波器上进行显示。
实验结果如下图所示：（插入示波器显示的正弦信号和方波信号的图片）实验分析：通过示波器的显示结果可以看出，正弦信号在时域上呈现周期性的波形，而方波信号则具有稳定的上下跳变。
这体现了正弦信号和方波信号在时域上的不同特征。
实验二：信号的采样和重构在数字信号处理中，信号的采样是将连续时间信号转化为离散时间信号的过程，信号的重构则是将离散时间信号还原为连续时间信号的过程。
在实际应用中，信号的采样和重构对信号处理的准确性至关重要。
实验设置如下：1. 生成一个正弦信号：频率为1kHz，振幅为1。
2. 设置采样频率为8kHz。
3. 对正弦信号进行采样，得到离散时间信号。
4. 对离散时间信号进行重构，得到连续时间信号。
5. 将重构的信号通过DAC输出到示波器上进行显示。
实验结果如下图所示：（插入示波器显示的连续时间信号和重构信号的图片）实验分析：通过示波器的显示结果可以看出，重构的信号与原信号非常接近，并且能够还原出原信号的形状和特征。
这说明信号的采样和重构方法对于信号处理的准确性有着重要影响。
最新数字信号处理实验报告
最新数字信号处理实验报告一、实验目的本次实验旨在加深对数字信号处理（DSP）理论的理解，并通过实践操作掌握数字信号处理的基本方法和技术。
通过实验，学习如何使用相关软件工具进行信号的采集、分析、处理和重构，提高解决实际问题的能力。
二、实验内容1. 信号采集与分析- 使用数字示波器采集模拟信号，并将其转换为数字信号。
- 利用傅里叶变换（FFT）分析信号的频谱特性。
- 观察并记录信号的时域和频域特性。
2. 滤波器设计与实现- 设计低通、高通、带通和带阻滤波器。
- 通过编程实现上述滤波器，并测试其性能。
- 分析滤波器对信号的影响，并调整参数以优化性能。
3. 信号重构实验- 应用所学滤波器对采集的信号进行去噪处理。
- 使用逆傅里叶变换（IFFT）重构经过滤波处理的信号。
- 比较重构信号与原始信号的差异，评估处理效果。
三、实验设备与材料- 计算机及DSP相关软件（如MATLAB、LabVIEW等）- 数字示波器- 模拟信号发生器- 数据采集卡四、实验步骤1. 信号采集- 连接并设置好数字示波器和模拟信号发生器。
- 生成一系列不同频率和幅度的模拟信号。
- 通过数据采集卡将模拟信号转换为数字信号。
2. 滤波器设计- 在DSP软件中设计所需的滤波器，并编写相应的程序代码。
- 调整滤波器参数，如截止频率、增益等，以达到预期的滤波效果。
3. 信号处理与重构- 应用设计的滤波器对采集的数字信号进行处理。
- 利用IFFT对处理后的信号进行重构。
- 通过对比原始信号和重构信号，评估滤波器的性能。
五、实验结果与分析- 展示信号在时域和频域的分析结果。
- 描述滤波器设计参数及其对信号处理的影响。
- 分析重构信号的质量，包括信噪比、失真度等指标。
六、实验结论- 总结实验中所学习到的数字信号处理的基本概念和方法。
- 讨论实验中遇到的问题及其解决方案。
- 提出对实验方法和过程的改进建议。
七、参考文献- 列出实验过程中参考的书籍、文章和其他资源。
数字信号处理实验报告_完整版
实验1 利用DFT 分析信号频谱一、实验目的1.加深对DFT 原理的理解。
2.应用DFT 分析信号的频谱。
3.深刻理解利用DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。
二、实验设备与环境 计算机、MATLAB 软件环境 三、实验基础理论1.DFT 与DTFT 的关系有限长序列 的离散时间傅里叶变换 在频率区间 的N 个等间隔分布的点 上的N 个取样值可以由下式表示：212/0()|()()01N jkn j Nk N k X e x n eX k k N πωωπ--====≤≤-∑由上式可知，序列 的N 点DFT ,实际上就是 序列的DTFT 在N 个等间隔频率点 上样本 。
2.利用DFT 求DTFT方法1：由恢复出的方法如下：由图2.1所示流程可知：101()()()N j j nkn j nN n n k X e x n eX k W e N ωωω∞∞----=-∞=-∞=⎡⎤==⎢⎥⎣⎦∑∑∑ 由上式可以得到：IDFTDTFT( )12()()()Nj k kX e X k Nωπφω==-∑ 其中为内插函数12sin(/2)()sin(/2)N j N x eN ωωφω--= 方法2：实际在MATLAB 计算中，上述插值运算不见得是最好的办法。
由于DFT 是DTFT 的取样值，其相邻两个频率样本点的间距为2π/N ，所以如果我们增加数据的长度N ，使得到的DFT 谱线就更加精细，其包络就越接近DTFT 的结果，这样就可以利用DFT 计算DTFT 。
如果没有更多的数据，可以通过补零来增加数据长度。
3.利用DFT 分析连续信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续信号离散化，这里需要进行两个操作：一是采样，二是截断。
对于连续时间非周期信号,按采样间隔T 进行采样，阶段长度M ，那么：1()()()M j tj nT a a a n X j x t edt T x nT e ∞--Ω-Ω=-∞Ω==∑⎰对进行N 点频域采样，得到2120()|()()M jkn Na a M kn NTX j T x nT eTX k ππ--Ω==Ω==∑因此，可以将利用DFT 分析连续非周期信号频谱的步骤归纳如下： （1）确定时域采样间隔T ，得到离散序列（2）确定截取长度M ，得到M 点离散序列,这里为窗函数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
现代数字信号处理实验报告1、估计随机信号的样本自相关序列。
先以白噪声()x n 为例。
(a) 产生零均值单位方差高斯白噪声的1000个样点。
(b) 用公式：9991ˆ()()()1000x n r k x n x n k ==-∑估计()x n 的前100个自相关序列值。
与真实的自相关序列()()x r k k δ=相比较，讨论你的估计的精确性。
(c) 将样本数据分成10段，每段100个样点，将所有子段的样本自相关的平均值作为()x n 自相关的估值，即：999001ˆ()(100)(100) , 0,1,...,991000x m n r k x n m x n k m k ===+-+=∑∑与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列()()x r k k δ=吗？(d) 再将1000点的白噪声()x n 通过滤波器11()10.9H z z -=-产生1000点的y (n )，试重复(b)的工作，估计y (n )的前100个自相关序列值，并与真实的自相关序列()y r k 相比较，讨论你的估计的精确性。
仿真结果：(a)图1.1 零均值单位方差高斯白噪声的1000个样本点分析图1.1：这1000个样本点是均值近似为0，方差为1的高斯白噪声。
(b)图1.2 ()x n的前100个自相关序列值分析上图可知：当k=0时取得峰值，且峰值大小比较接近于1，而当k≠0时估计的自相关值在0附近有小幅度的波动，这与真实自相关序列r(k)=δ(k)x比较接近，k≠0时估计值非常接近0，说明了估计的结果是比较精确的。
(c)图1.3基于Bartlett 法的前100个自相关序列值与(b)的结果相比，同样在k=0时达到峰值，k ≠0时0值附近上下波动；估计值的方差比较小，随着k 的增大波动幅度逐渐变小，在k 较大时它更接近真实自相关序列()()x r k k δ=。
即采用分段方法得到的自相关序列的估计值更加接近r x (k)=δ(k)。
分析仿真图也可以看出：将样本数据分段，将所有子段的样本自相关的平均值作为()x n 自相关的估值时，可以有效的降低自相关估计的方差，而分段样本估计的优点在于，估计自相关序列与实际自相关序列的方差减小，且当分段数越大，估计值越趋向于无偏估计。
(d)图1.4 y(n)的前100个自相关序列值与真实值的对比从图中可以看出在k=0时估计与真实的自相关序列之间有较小的误差，随着k的增大，估计得到的值有较大的波动，存在一定误差。
源程序clcclear%%产生1000个高斯白噪声的样本点x=randn(1,1000);K=1000;figure(1);k=0:K-1;stem(k,x,'.'); %绘制1000个高斯白噪声title('零均值单位方差高斯宝噪声，1000个样本点');xlabel('k');ylabel('x[k]');mean_x=mean(x)var_x=var(x)%%for k=0:99for n=k+1:1000y_ess(n)=x(n)*x(n-k);endr_ess(k+1)=sum(y_ess)/1000;endfigure(2);k=[0:99];stem(k,r_ess,'.');title('根据样本点估计出的前100自相关序列值');xlabel('k');ylabel('r_ess[k]');hold on;realvalue=[1,zeros(1,99)];stem(k,realvalue,'r','.');legend('根据样本点估计出的前100自相关序列值','真实的自相关序列');error1=r_ess-realvalue;mean_error_b=mean(error1)var_error_b=var(error1)%%for k=0:99for m=0:9for n=k+1:100y_ess2(m+1,n)=x(n+100*m)*x(n-k+100*m);endendr_ess2(k+1)=sum(sum(y_ess2))/1000;endfigure(3);k=0:99;stem(k,r_ess2,'b.');hold on;realvalue2=[1,zeros(1,99)];stem(k,realvalue2,'r.','.');title('Bartlett法估计功率谱方法得出的前100个自相关序列值');xlabel('k');ylabel('r_ess2[k]');legend('Bartlett法估计功率谱方法得出的前100个自相关序列值','真实的自相关序列');error2=r_ess2-realvalue2;mean_error_c=mean(error2)var_error_c=var(error2)%%y=zeros(1,1000);B=[1]; A=[1,-0.9];y=filter(B,A,x);r_ess3=zeros(1,100); for k=0:99for n=(k+1):1000r_ess3(k+1)=r_ess3(k+1)+y(n)*y(n-k); endr_ess3(k+1)=r_ess3(k+1)/1000; endfigure(4);stem(r_ess3,'.');title('y[n]前100个自相关序列估计值'); xlabel('k'),ylabel('r_ess3(k)'); hold on;p=[1,zeros(1,99)]; h=filter(B,A,p); for i=1:100h1(i)=h(101-i); endrh=conv(h,h1); rh=rh(100:199);realvalue3=conv(p,rh);realvalue3=realvalue3(1:100); stem(realvalue3,'r.','.');legend('y[n]前100个自相关序列估计值','y[n]的真实自相关序列');2、计算机练习2：AR 过程的线性建模与功率谱估计。
考虑AR 过程：()(1)(1)(2)(2)(3)(3)(4)(4)(0)()x n a x n a x n a x n a x n b v n =-+-+-+-+()v n 是单位方差白噪声。
(a) 取b (0)=1, a (1)=0.-7348, a (2)=-1.8820, a (3)=-0.7057, a (4)=-0.8851，产生x (n )的N =256个样点。
(b) 计算其自相关序列的估计ˆ()x rk ，并与真实的自相关序列值相比较。
(c) 将ˆ()x rk 的DTFT 作为x (n )的功率谱估计，即： 1211ˆˆ()()|()|N j jk j xx k N P e rk e X e Nωωω--=-+==∑。
(d) 利用所估计的自相关值和Yule-Walker 法(自相关法)，估计(1), (2), (3), (4)a a a a 和(0)b 的值，并讨论估计的精度。
(e) 用(d)中所估计的()a k 和(0)b 来估计功率谱为：2241|(0)|ˆ()1()j xjk k b P e a k e ωω-==+∑。
(f) 将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。
(g) 重复上面的(d)~(f)，只是估计AR 参数分别采用如下方法：(1) 协方差法；(2) Burg 方法；(3) 修正协方差法。
试比较它们的功率谱估计精度。
仿真结果： （a ）图2.1 x (n )的N =256个样点（b ）图2.2自相关序列的估计值与真实的对比图2.2中估计的自相关序列与真实的自相关序列差异较大；在k>100后，真实的自相关序列就波动得很小，而估计的自相关序列则仍有较大的波动，但总体上来言两者都随着k的增大而逐渐衰减，当k>150时，真实自相关序列逐渐趋于0，而估计出的自相关序列却仍有较大的波动，这是因为估计的点数较少，使得估计精度不够。
（c）图2.3 估计的功率谱与真实功率谱对比（d）Yule-Walker法(自相关法)估计的参数值为：b(0)= 1.1537[a(1) a(2) a(3) a(4)]=[-0.7174 -1.7952 -0.6387 -0.8167]图2.4 Yule-Walker法估计的功率谱与真实功率谱对比Yule-Walker法(自相关法)估计的参数，其系数的符号正确，但数值大小相差较大，并且多次实验的相差值也很大，参数估计的精度远远不够。
因此从图2.4中也能看出，该方法得到功率谱与真实的谱相差很大（e）协方差法图2.5 协方差法估计的功率谱与真实功率谱对比采用协方差法估计的参数，其系数与真实的系数非常接近，该方法能够较精确的估计出功率谱。
（f）修正协方差图2.6 修正的协方差法估计的功率谱与真实功率谱对比采用修正的协方差法估计的参数，其系数虽然没有协方差法和burg法那么精确，但是估计误差也很小。
从图2.6中也能看出，该方法能够较精确的估计出功率谱。
（g）Burg算法图2.7 burg法估计的功率谱与真实功率谱对比采用burg估计的参数，其系数几乎等于真实的系数，分析图2.7中也能看出，该方法非常精确的估计出功率谱，几乎与真实的功率谱曲线重合。
源程序：clc;clear;N=256;NN=20000;v1=normrnd(0,1,50,NN);v=v1(:,1:N);r=zeros(1,N);r1=zeros(1,N);w=0:2*pi/100:2*pi;P=zeros(1,length(w));PP1=zeros(1,length(w));PP2=zeros(1,length(w));PP3=zeros(1,length(w));PP4=zeros(1,length(w));for s=1:50x1=filter([1],[1,0.7348,1.8820,0.7057,0.8851],v1(s,:)); x=x1(1:N);for k=1:Nrx(k)=0;for n=k:Nrx(k)=rx(k)+x(n)*x(n-(k-1));endrx(k)=rx(k)/(N);endr=r+rx;for k=1:Nrx1(k)=0;for n=k:NNrx1(k)=rx1(k)+x1(n)*x1(n-(k-1));endrx1(k)=rx1(k)/(NN);endr1=r1+rx1;for i=1:length(w)P(i)=P(i)+rx(1);for n=2:NP(i)=P(i)+rx(n)*exp(-j*(n-1)*w(i))+rx(n)*exp(j*(n-1)*w(i)); endendA=toeplitz(rx(1:4)',rx(1:4));B=-rx(2:5)';Ap=A\B;b0=rx(1);for i=1:4b0=b0+Ap(i)*rx(i+1);endb0=sqrt(b0);for i=1:length(w)P1(i)=1;for k=1:4P1(i)=P1(i)+Ap(k)*exp(-j*k*w(i));endP1(i)=b0^2/abs(P1(i))^2;endPP1=PP1+P1;A=[];for k=1:4c=[];for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n-k);endc=[c;rr];endA=[A c];endB=[];for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n);endB=[B;rr];endB=B*(-1);Ap=A\B;for i=1:length(w)P2(i)=1;for k=1:4P2(i)=P2(i)+Ap(k)*exp(-j*k*w(i));endP2(i)=x(1)^2/abs(P2(i))^2;endPP2=PP2+P2;A=[];for k=1:4c=[];for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n-k)+x(n-4+l)*x(n-4+k); endc=[c;rr];endA=[A c];endB=[];for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n)+x(n-4+l)*x(n-4);endB=[B;rr];endB=B*(-1);Ap=A\B;for i=1:length(w)P3(i)=1;for k=1:4P3(i)=P3(i)+Ap(k)*exp(-j*k*w(i));endP3(i)=x(1)^2/abs(P3(i))^2;endPP3=PP3+P3;p=4;ef=zeros(1+p,N);eb=zeros(1+p,N);ef(1,:)=x;eb(1,:)=x;km=[];for m=2:p+1mol=0;den=0;for n=m:Nmol=mol+(-2)*ef(m-1,n)*eb(m-1,n-1);den=den+(ef(m-1,n))^2+(eb(m-1,n-1))^2; endkm(m-1)=mol/den;for n=m:Nef(m,n)=ef(m-1,n)+km(m-1)*eb(m-1,n-1); eb(m,n)=eb(m-1,n-1)+km(m-1)*ef(m-1,n); endendaa=[1];for i=1:4aa=[aa;0];bb=aa(length(aa):-1:1);aa=aa+bb*km(i);endfor i=1:length(w)P4(i)=1;for k=2:5P4(i)=P4(i)+aa(k)*exp(-j*(k-1)*w(i));endP4(i)=1/abs(P4(i))^2;endPP4=PP4+P4;endfigure(1)stem(1:N,x,'.');title('x[n]的256个样本点');xlabel('n');ylabel('x[n]');figure(2)plot(0:N-1,r/50); hold on;plot(0:N-1,r1/50,'r');title('x[n]的256个样本点估计出的前256个自相关序列和真实值'); ylabel('Rx(k)');xlabel('k');legend('估计值','真实值');grid on;axis([0 256 -40 40]);hold off;figure(3)plot(w/pi,10*log10(P/50)); hold on;title('功率谱估计');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(PP1/50),'r');plot(w/pi,10*log10(PP2/50),'g');plot(w/pi,10*log10(PP3/50),'y');plot(w/pi,10*log10(PP4/50),'k');aap=[0.7348,1.8820,0.7057,0.8851];for i=1:length(w)P5(i)=1;for k=1:4P5(i)=P5(i)+aap(k)*exp(-j*k*w(i)); endP5(i)=1/abs(P5(i))^2;endplot(w/pi,10*log10(P5),':');legend('Rx(k)的DTFT','Yule-Walker');grid on;hold off;figure(4)plot(w/pi,10*log10(P/50)); hold on;title('功率谱估计比较');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('Rx(k)的DTFT','真实频谱');grid on;hold off;figure(5)plot(w/pi,10*log10(PP1/50)); hold on;title('Yule-Walker法功率谱估计比较'); ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('Yule-Walke法','真实频谱');grid on;hold off;figure(6)plot(w/pi,10*log10(PP2/50)); hold on;title('协方差法功率谱估计比较');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('协方差法','真实频谱');grid on;hold off;figure(7)plot(w/pi,10*log10(PP3/50)); hold on;title('修正协方差法功率谱估计比较');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('修正协方差法','真实频谱');grid on;hold off;figure(8)plot(w/pi,10*log10(PP4/50)); hold on;title('Burg 法功率谱估计比较');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('Burg 法','真实谱');grid on;hold off;3、计算机练习3：维纳噪声抑制(例6.6的扩展实验)。