现代信号处理大作业

信号时频分析技术及matlab仿真电路与系统王冠军 201128013926153摘要:本文介绍了时频分析的一些基础理论,对短时傅里叶变换Wigner-Ville分布做了简单介绍，运用MATLAB语言实现了旨在构造一种时间和频率的密度函数,以揭示信号中所包含的频率分量及其演化特性的wigner-ville分布。

并对时频分析方法的优缺点进行了分析。

关键词:时频分析短时傅里叶变换wigner-ville分布1 引言基于Fourier变换的传统信号处理技术从信号频域表示及能量的频域分布的角度揭示了信号在频域的特征。

但Fourier变换是一种整体变换，只能为人们提供信号在时域或频域的全局特性而无法了解信号频谱随时间变化的情况。

因此，需要使用一种时间和频率的联合函数来表示信号，这种表示简称为信号，也就是信号的时频分析。

2 时频分析方法信号时频分析主要研究非平稳信号或时变信号的频谱含量是怎样随时间变化的。

时频分析是当今信号处理领域的一个主要研究热点，目前常用时频分析方法主要有短时傅里叶变换、Gabor展开、小波变换、Wigner-Ville分布。

本文主要介绍了短时傅里叶变换和Wigner-Ville分布两种分析方法。

2．1 短时傅立叶变换STFT从历史上看,信号的时频分析用的最多的是短时傅立叶变换,这种变换的基本思想是用一个窗函数乘时间信号,该窗函数的时宽足够窄,使取出的信号可以被看成是平稳的,然后进行的傅立叶变换可以反映该时宽中的频谱变化规律,如果让窗随时间轴移动,可以得到信号频谱随时间变化的规律。

对于时变信号,了解不同时刻附近的频域特征是至关重要的。

因此,人们采用时间—频率描述时变信号,将一维的时域信号映射到一个二维的时域平面,全面反映观测信号的时频联合特征。

短时傅立叶变换反映了这一思想,对于时变信号,采用某一滑动窗函数截取信号,并认为这些信号是准平稳的,然后,再分别对其进行傅立叶变换,构成时变信号的时变谱。

《现代信号处理》大作业姓名：王成志学号：1140349078一. L D 迭代算法的matlab 实现1.1 Levinson-Durbin 算法介绍功率谱估计大致可以分为经典谱估计和现代功率谱估计，经典谱估计方法存在着以下三点缺陷：（1）数据加窗或自相关加窗，都隐含着假定在窗外未观测到的数据或自相关系数为零，该假设不切实际。

（2）要性能好往往需要较长的数据，但实际数据长度有限（3）窗函数容易造成谱的模糊。

采用AR 模型的现代谱估计方法可以克服这些不足。

其中LD 递推算法可以在计算机上方便实现。

LD 递推算法具体计算步骤如下：（1） Yule-Walker 方程的矩阵形式(1)所示：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--+-----001)0()2()1()()1()1()0()1()()2()1()0(2,1,σk k k xx xx xx xx x xx xx xx xx xx xx xx a a r k r k r k r k r r r r k r r r r 系数矩阵xx Hxx R R =,为Hermitian 矩阵，对角线上元素相同，即为Topliez 矩阵。

（2） P-1阶Yule-Walker 方程为：21111(0)(1)(1),1(1)(0)(2)0,1(1)(2)(0)0x x x p p x x x p x x x R R R p a R R R p a p R p R p R σ-----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎢⎥⎣⎦⎣⎦⎣⎦ 其中，2211{()}p p E e l σ--=为误差功率。

写成联立方程：2111,0,0()0,1,,1p pp k xk m a R m k m p σ---=⎧=-=⎨=-⎩∑ 取共轭得：21**11,0,0()0,1,,1p pp kxk m aR m k m p σ---=⎧=-=⎨=-⎩∑变量替换，并利用*()()x x R l R l =得：21*11,10,1()0,0,,2p pp p kx k m p aR m k m p σ-----=⎧=--=⎨=-⎩∑ 表示成矩阵：*1*1210(0)(1)(1),10(1)(0)(2),2(1)(2)(0)1x x x p x x x p p x x x R R R p a p R R R p a p R p R p R σ-----⎡⎤⎡⎤⎡⎤-⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎣⎦⎣⎦⎣⎦ 求解得：*.1,1,,0,,p k p k p p p k a a K a k p ---=+=22*1p p p p K σσ-=+∆ 2210p p p K σ-=∆+,p p p K a =222*22111[][1]p p p p p p p K K K σσσσ---=+-=-（3） 当k=1时，即一阶递推为：⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-01)0()1()1()0(211,1σa R R R R x x x x求解可得：)1()0()0()1( ,11,1211,10,1x x x x R a R R R a a +=-==σ（4） 对于2≥p 时，递推为:10,≡p a , *,1,1,k p p p k p k p aK a a ---+=, ]1[2212p p p K -=-σσ 21,-∆-==p pp p p a K σ∑-=--+=∆11,1)()(p k x kp x p k p R ap R矩阵R x 已知，可得到各阶AR 模型系数为：)0())1(1( ,)0()1()1(2111xx xx xx r a r r a -=-=ρ11111)()()()(--=--∑-+-=∆-=k k l xx k xx k kk l k r l a k r k a ρρ1,,2,1)()()()(*11-=-+=--k i i k a k a i a i a k k k k12))(1(--=k k k k a ρρ1.2实验结果(1) 输入p=3，rr = [70,60,50,40] 时，求得AR 模型估计参数为:a =1.0000 -0.8571 0 0 1.0000 -0.5275 -0.3846 0 1.0000 -0.7572 -0.6996 0.5972 各阶求得的方差为：sigma = 18.5714 15.8242 10.18013阶时，a 3 (1)= -0.7572 a 3 (2)= -0.6996 a 3 (3)= -0.5972(2) 输入p=5，rr = [30,45,26,33,47,43]时，AR 模型估计参数为：a =1.0000 -1.5000 0 0 0 0 1.0000 0.2800 -1.1867 0 0 0 1.0000 0.8227 -1.3147 -0.4573 0 0 1.0000 1.9708 1.9858 -2.5226 -2.5105 0 1.0000 1.0869 1.0977 -1.8235 -1.8166 0.3521 各阶求得的方差为： sigma =37.5000 15.3067 12.1054 64.1881 56.23165阶时， a 5 (1)= 1.0869 a 5(2)= 1.0977 a 5(3)= -1.8235 a 5(4)= -1.8166 a 5(5)= 0.3521二． 一维平稳信号由两个高斯信号叠加而成12241122()()[exp(())exp(())]22z t t t j t t t j t αααωωπ=--++--+，其中12,t t >12ωω>，分别求出()z t 的WV 分布及其模糊函数，画出二者的波形图，指出并分析其信号项和交叉项。

研究生“现代信号处理”课程大型作业（以下四个题目任选三题做）1. 请用多层感知器（MLP ）神经网络误差反向传播（BP ）算法实现异或问题（输入为[00;01;10;11]X T =，要求可以判别输出为0或1），并画出学习曲线。

其中，非线性函数采用S 型Logistic 函数。

2. 试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。

滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。

3. 根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1） Levinson 算法2） Burg 算法3） ARMA 模型法4） MUSIC 算法4. 图1为均衡带限信号所引起失真的横向或格型自适应均衡器(其中横向FIR 系统长M =11), 系统输入是取值为±1的随机序列)(n x ，其均值为零；参考信号)7()(-=n x n d ；信道具有脉冲响应：12(2)[1cos()]1,2,3()20 n n h n W π-⎧+=⎪=⎨⎪⎩其它式中W 用来控制信道的幅度失真(W = 2～4, 如取W = 2.9,3.1,3.3,3.5等)，且信道受到均值为零、方差001.02=v σ(相当于信噪比为30dB)的高斯白噪声)(n v 的干扰。

试比较基于下列几种算法的自适应均衡器在不同信道失真、不同噪声干扰下的收敛情况(对应于每一种情况,在同一坐标下画出其学习曲线):1） 横向/格-梯型结构LMS 算法2） 横向/格-梯型结构RLS 算法并分析其结果。

图1 横向或格-梯型自适应均衡器参考文献[1] 姚天任, 孙洪. 现代数字信号处理[M]. 武汉: 华中理工大学出版社, 2001[2] 杨绿溪. 现代数字信号处理[M]. 北京: 科学出版社, 2007[3] S. K. Mitra. 孙洪等译. 数字信号处理——基于计算机的方法(第三版)[M]. 北京: 电子工业出版社, 2006[4] S.Haykin, 郑宝玉等译. 自适应滤波器原理(第四版)[M].北京: 电子工业出版社, 2003[5] J. G. Proakis, C. M. Rader, F. Y. Ling, etc. Algorithms for Statistical Signal Processing [M].Beijing: Tsinghua University Press, 2003一、请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为[00;01;10;11]，要求可以判别输出为0或1），并画出学习曲线。

现代信号处理作业1.(5″)证明下面定理：任何一个无偏估计子方差的下界叫作Cramer-Rao 下界 定理：令1(,,)N x x x =为一样本向量，(|)f x θ是x 的条件密度，若ˆθ是θ的一个无偏估计子，且(|)/f x θθ∂∂存在，则221ˆˆvar()()[ln (|)]E E f x θθθθθ=-≥∂∂式中ˆln (|)()()f x K θθθθθ∂=-∂。

其中()K θ是θ的某个不包含x 的正函数。

2.(10″)Wiener 滤波是信号处理中最常用和基础的波形估计工具之一，对其在自己研究领域的应用情况进行一个简单综述。

3.(5″)二阶滑动平均过程由2()()1(1)2(2),{()~(0,)}x n w n b w n b w n w n N σ=+-+-定义，式中2(0,)N σ表示正态分布，其均值为零、方差为2σ。

求x(n)的功率谱。

4.(20″)信号的函数表达式为：()sin(2100)sin(2300)()sin(2200)()()x t t t A t t dn t n t πππ=++++，其中，A(t)为一随时间变化的随机过程，dn(t)为经过390-410Hz 带通滤波器后的高斯白噪声，n(t)为高斯白噪声，采样频率为1kHz ，采样时间为2.048s 。

(1) 利用现代信号处理知识进行信号的谱估计； (2) 利用现代信号处理知识进行信号的频率提取； (3) 分别利用Wiener 滤波和Kalman 滤波进行去噪； (4) 利用Wigner-Ville 分布分析信号的时频特征。

5.(10″)附件中表sheet1 为某地2008年4月28日凌晨12点至2008年5月4日凌晨12点的电力系统负荷数据，采样时间间隔为1小时，利用ARMA 方法预测该地5月5日的电力系统负荷，并给出预测误差（5月5日的实际负荷数据如表sheet2）。

1、定理：令1(,,)N x x x =为一样本向量，(|)f x θ是x 的条件密度，若参数估计ˆθ是真实参数θ的一个无偏估计子，且(|)/f x θθ∂∂、22(|)/f x θθ∂∂存在，则ˆθ的均方误差所能达到的下界（称为Cramer-Rao 下界）等于Fisher 信息的导数，即：221ˆˆvar()()[ln (|)]E E f x θθθθθ=-≥∂∂ （1-1）不等式中等号成立的充分必要条件是：ˆln (|)()()f x K θθθθθ∂=-∂ （1-2） 其中()K θ是θ的某个正函数，与样本1(,,)N x x x =无关。

现代信号及其应用雷达信号处理技术发展综述姓名:******班级****** 学号******日期：2015/4/10摘要本文简单介绍了雷达发展、雷达工作原理、雷达分类、航迹、航迹处理、滤波器、滤波器工作原理以及滤波器在航迹处理中的应用。

ABSTRACTThis paper simply introduces the development of radar, the working principle of radar, the classification of radar,the track, the tracking processing, application of filter, the filter and filter in the working principle of the track processing.目录1前言 (4)2雷达 (4)2.1雷达意义简述 (4)2.2雷达分类 (4)2.3.1雷达原理 (4)2.3.2影响雷达监测目标的影响因素 (4)2.3.3探测范围 (4)2.3.4探测距离 (5)2.3.5发射距离 (5)2.3.6发射频率 (5)3航迹 (5)3.1航迹形成过程 (5)3.2航迹处理 (5)3.3.1滤波器 (5)3.3.2滤波器的工作原理 (6)3.3.3滤波器在航迹处理过程中的应用 (6)4总结 (7)雷达信号处理技术发展综述1 前言自第二次世界大战以来，雷达技术得到了迅猛的发展，雷达家族不断扩大。

雷达的用途早已超出了检测与定位的范畴。

特别是为适应现代战争和现代经济发展的需要研制的不同用途的雷达，更是种类繁多，五花八门，本文下面简单介绍雷达的发展史以及滤波器在航迹处理过程中的应用的部分内容。

2 雷达2.1 雷达意义简述雷达为英文Rader 一词的译音，来源于词组Radio Detection and Ranging ，意为无线电监测与定位，即雷达是利用无线电方法来发现目标并测定目标位置参数的一种电子系统。

1、已知X a (t) 2COS (2f o t)式中f o =1OOH Z,以采样频率f s =400Hz 对X a (t)进行采样，得到采样信号X a (t)和时域离散信号X(n)，试完成下面各题： (1)写出X a (t)的傅里叶变换表示式 X a (j )；(2) 写出X a (t)和x(n)的表达式；(3 )分别求出X a (t)的傅里叶变换和x(n)的傅里叶变换。

解:( 1)j tj tX a (j )X a (t)e j dt 2cos( o t)e j dt3、在时域对一有限长的模拟信号以 4KHZ 采样，然后对采到的N 个抽样做N 点DFT ，所得离散谱线的间距相当于模拟频率100HZ 。

某人想使频率能被看得清楚些，每50HZ 能有一根谱线，于是他用8KHZ 采样，对采到的2N 个样点做2N 点DFT 。

问：他的目的能达到吗？ 答：不能，因为他忽略了数字频率和模拟频率的区别。

提高采样频率f s ，N 固然大了，数字频率(单位圆)上的样点数确实增加了，但从模拟频率谱看，样点一点也没有变得密集，这是因为数字频率2总是对应模拟频率 f s 。

2 f sf s采样频率由f s 到2 f s 增加一倍，N 也增加一倍，但模拟频率的采样间隔s s100Hz2N N 2)，不能提 2N高模拟频率的分辨(e j 0t e j 0t上式中指数函数和傅里叶变换不存在, X a (j ) 2 [ ()(2)x a (t )X a (t) (t)e j t dt引入奇异函数)]函数，它的傅里叶变换可以表示成:nT)2cos( 0nT) (tnT)n2cos( 0nT),2、用微处理器对实数序列作谱分析,以下各参数：(1)x(n) 最小记录时间 (2) (3) (4) 解:( 1)Tpmin T max N min要求谱分辨率F最大取样时间 最少采样点数 在频带宽度不变的情况下将频率分辨率提高一倍的已知(2)F 50Hz1 F 1 T pmin10.02s 50 1(3)N min1 s minT PT 0.02s0.5 10 3s(4)辩率提高1倍(F 变成原来的12)T p0.04sN min~T 0.5 10 s频带宽度不变就意味着采样间隔 5OHZ ，信号最高频率1KHz,是确定N 值。

现代信号处理作业现代信号处理课程作业1.做⼀个⽹络检索,简述现代信号处理技术的主要特征和技术特点,并阐述信号处理在实际⼯程中的应⽤情况代信号处理技术的主要特征和技术特点:1)精度⾼:在模拟系统的电路中,元器件精度要达到10-3以上已经不容易了,⽽数字系统17位字长可以达到10-5的精度,这是很平常的?例如,基于离散傅⾥叶变换的数字式频谱分析仪,其幅值精度和频率分辨率均远远⾼于模拟频谱分析仪?2) 灵活性强:数字信号处理采⽤了专⽤或通⽤的数字系统,其性能取决于运算程序和乘法器的各系数,这些均存储在数字系统中,只要改变运算程序或系数,即可改变系统的特性参数,⽐改变模拟系统⽅便得多?3) 可以实现模拟系统很难达到的指标或特性:例如:有限长单位脉冲响应数字滤波器可以实现严格的线性相位;在数字信号处理中可以将信号存储起来,⽤延迟的⽅法实现⾮因果系统,从⽽提⾼了系统的性能指标;数据压缩⽅法可以⼤⼤地减少信息传输中的信道容量?4)可以实现多维信号处理:利⽤庞⼤的存储单元,可以存储⼆维的图像信号或多维的阵列信号,实现⼆维或多维的滤波及谱分析等?信号处理在实际⼯程中的应⽤情况:数字信号处理是利⽤计算机或专⽤计算机或专⽤处理设备,以数据形式对信号进⾏采集,变换,滤波,估值,增强,压缩,识别等处理,以得到符合⼈们需要的信号形式?数字信号处理是以众多科学为理论基础的,他所涉及的范围及其⼴泛?DSP 技术应⽤到我们的⽣活的每⼀个⾓落,从军⽤到民⽤,从航空航天到⽣产⽣活,都越来越多地使⽤DSP. DSP技术在航空⽅⾯,主要⽤于雷达和声纳信号处理;在通信⽅⾯,主要⽤于移动电话,IP电话,ADSL和HFC的信号传输;在控制⽅⾯,主要⽤于电机控制,光驱和硬盘驱动器;在测试/测量⽅⾯,主要⽤于虚拟仪器,⾃动测试系统,医疗诊断等;在电⼦娱乐⽅⾯,主要⽤于⾼清晰度电视,机顶盒,家庭影院,DVD 等应⽤;还有数字相机,⽹络相机等等都应⽤了SP技术?同时,SOC芯⽚系统,⽆线应⽤,嵌⼊式DSP都是未来DSP的发展⽅向和趋势?可以说,没有DSP就没有对互联⽹的访问,也不会有多媒体,也没有⽆线通信?因此DSP仍将是整个半导体⼯业的技术驱动⼒?现在,DSP应⽤领域不断拓宽,其涵盖⾯包括宽带Internet接⼊业务,下⼀代⽆线通信系统的发展,数字消费电⼦市场,汽车电⼦市场的发展等诸多多⽅⾯?现代数字信号处理器是执⾏⾼速数字信号系统的IC电路,它恰好适合多媒体信息化社会需求,迅速发展壮⼤?如今,世界电⼦器件市场上,各种各样的DSP器件已相当丰富?⼤⼤⼩⼩封装形式的DSP器件,已⼴泛⽤于各种产品的⽣产领域,⽽且DSP的应⽤领域仍在不断的扩⼤,发展速度异常?2?简述信号的频率分析技术及其应⽤,阐述实现精细频率分析的实现⽅法?考虑到数字信号分析中,虽然提⾼信号的采样频率可以改善信号分析的频率分辨率,但是提⾼信号的采样频率通常需要付出额外的硬件代价,往往受制于可实现性与成本问题⽽难以实现?因此,就需要使⽤频谱细化技术在尽可能低的采样频率下提⾼数字信号分析的频率分辨率的措施?频谱细化的基本思路是对信号频谱中的某⼀频段进⾏局部放⼤,也即在某⼀频率附近局部增加谱线密度,实现选带频段分析?频谱细化技术在⽣产实践和科学研究中获得了⽇益⼴泛的应⽤?例如,齿轮箱的故障诊断要求准确分辨齿轮各阶啮合振动的主频和边频等,其频谱图上的频率间隔很细,但频率分布⼜较宽,为了识别谱图的细微结构,就必须对信号进⾏细化分析;直升机?坦克?巡航导弹的声⾳具有显著的⾮平稳性,为了得到准确的时延量,信号的取样不能太长,⽽FFT计算的频谱存在栅栏效应?因此必须采⽤有效的⽅法对频谱进⾏细化,这样才能保证⾜够的相关计算精度;在⽆线电通信信号和其他的实际⼯程信号的分析中,为了获取更⾼的测量精度和实时检测能⼒,需要对信号频谱进⾏细化分析,以提供有⽤信息?因此对频谱细化技术的研究受到普遍重视,也是当前信号处理技术研究中的⼀个⼗分活跃的课题?常见的经典⽅法有:复调制细化法?Chirp-Z变换?FFT+FT细化法?DFT补零法等很多⽅法?复调制细化法:⼜称为选带频率细化选带频谱分析,是20世纪70年代发展起来的?其传统的分析步骤为:移频(复调制)低通滤波器重抽样--FFT及谱分析频率成分调整,因其物理概念⾮常明确,所以⼀直沿⽤⾄今?FFT+FT细化法:该⽅法的原理本质是将连续傅⾥叶变换经过将积分化成求和?时域离散化和时域截断为有限长三个步骤变换得到时间离散?频率连续的特殊傅⾥叶变换形式?FFT+FT连续细化分析傅⾥叶变换法先⽤FFT做全景谱,再对指定的⼀个频率区间进⾏细化计算:先确定频率分辨率,再确定计算频率序列,最后⽤FT连续谱分析⽅法进⾏实部和虚部计算,合成幅值谱和相位谱? Chirp-Z变换:最早提出于1969年,CZT是⼀种在Z平⾯上沿着螺旋线轨道计算有限时宽的Z变换⽅法?基本原理是在折叠频率范围内任意选择起始频率和频率分辨率在这有限带宽⾥对样本信号进⾏Z变换这与频谱校正⽅法中的FFT + FT 连续细化分析傅⾥叶变换法的基本原理是⼀样的?3、通过⽹络检索,对弱信号检测技术进⾏调研,分析⼀下现代弱信号检测的⽅法微弱信号检测(WeakSignalDetection)是⼀门新兴的技术学科,应⽤范围遍及光?电?磁?声?热?⽣物?⼒学?地质?环保?医学?激光?材料等领域?其仪器已成为现代科学研究中不可缺少的设备?微弱信号检测技术是采⽤电⼦学?信息论?计算机及物理学的⽅法,分析噪声产⽣的原因和规律,研究被测信号的特点与相关性,检测被噪声淹没的微弱有⽤信号?微弱信号检测的⽬的是从强噪声中提取有⽤信号,或⽤⼀些新技术和新⽅法来提⾼检测系统输出信号的信躁⽐?信号处理系统的信躁⽐改善等于输⼊(⽩)躁声带宽与系统的躁声等效带宽之⽐?因此,减少系统的躁声等效宽度便可以提⾼系统的输出信躁⽐?对于信躁⽐⼩于1的被躁声淹没的信号,只要信号处理系统的躁声等效带宽做得很⼩,就可以将信号(或信号携带的信息)从躁声中提取出来,这就是通常的微弱信号检测的指导思想之⼀?现代弱信号检测的⽅法和原理窄带滤波法: 使⽤窄带滤波器,滤掉宽带躁声只让窄带宽信号通过(仅有极少量窄带躁声通过)?窄带滤波法能减少躁声对有⽤信号的影响?滤除掉通频带以外躁声,提⾼信号的信躁⽐?但是,由于⼀般滤波器的中⼼频率不稳定,不能满⾜更⾼的滤除躁声的要求?双路消躁声法:由于信号与躁声性能完全不同,信号⼀般为⼀些变化规律已知的量,⽽躁声是⼀些随机量满⾜统计规律?当随机性的躁声从两路到达加法器时,极性正好相反,经过加法器相加后把躁声消掉?只有少数强躁声才通过阀值电路⽽产⽣本底计数,根据统计规律?本底计数时间较长时为恒定值?故可以先测出它,然后从总计数中把它减得到信号计数?这种⽅法只能检测到微弱的正弦信号是否存在,⽽不能复现信号波形?同步累积法:利⽤信号的重复性,躁声的随机性,对信号进⾏重复累积(⼏次),使SNIR提⾼,但需耗费时间?锁定接收法(频域分析法) :锁定检测法是利⽤互相关原理,使输⼊待测的周期信号与频率相同的参考相关器中实现互相关,从⽽将深埋在躁声中的周期信号携带的信息检测出来?相关检测法: 相关检测技术是应⽤信号周期性和噪声随机性的特点,通过⾃相关或互相关运算,达到去除躁声检测出信号的⼀种技术?由于信号和躁声是相互独⽴的过程,根据相关函数和互相关函数的定义,信号只与信号本⾝相关与躁声不相关??取样积分法:取样积分(或信号平均)法是将待测的重复信号逐点多次取样并进⾏同步积累,从⽽达到从噪声中恢复信号波形的⽅法?取样积分也采⽤同步相关检测的原理和⽅法,实现从噪声中提取信号,但它的参考信号只在窗⼝持续期间与被测信相关,每周相关时间很短,此外它的相移也是在很慢的变化?取样积分由单点取样积分与多点取样积分两种?4.利⽤MATLAB产⽣出⼀个线性调频信号(chirp信号),采样频率=8000Hz,持续时间1s,起始频率=500Hz,终⽌频率=1300Hz,给出其时域波形图,请利⽤短时FFT分析函数对数据进⾏时间-频率分析,观测频率随时间的变化情况分析结果:00.10.20.30.40.50.60.70.80.91-1-0.50.51时间t/s幅度线性调频信号Time F r e q u e n c y 线性调频信号的STFT 频谱图50010001500200025003000350015. 研究⼀下利⽤⾃相关实现含噪声的正弦信号检测⽅法,并利⽤MATLAB 进⾏验证:答:相关函数的应⽤很⼴,例如,噪声中信号的检测?信号中隐含周期性信号的检测,信号相关性的检测等?设信号)(n f 由正弦信号) (n x 加均值为零的⽩噪声)(n s 所组成,即)()()(n s n x n f +=;那么)(n f 的⾃相关为∑∞=++++=0)]()()][()([1)(n m n s m n x n s n x N m R=)()()()(m R m R m R m R ss sx xs xx +++其中)(m R xs 和)(m R sx 分别是正弦信号)(n x 和⽩噪声)(n s 的互相关?⽩噪声是随机的,和信号)(n x 应⽆相关性,所以)(m R xs 和)(m R sx 应趋近于零?⽩噪声)(n s 的⾃相关函数)(m R ss 主要在n=0处有值,当0||>n 时,衰减很快?由于)(n x 是周期函数,那么)(m R xx 将呈周期变化,从⽽揭⽰出隐含在)(m R xx 中的周期性?由于)(n x 总为有限长,所以这些峰值将是逐渐衰减的,且)(m R xx 的最⼤延迟应⼩于数据长度?01002003004005006007008009001000-4-224含噪声时域正弦信号01002003004005006007008009001000-0.500.5⾃相关检测出的正弦信号6. 简述⼩波滤波的原理,并利⽤MATLAB 中的⼩波⼯具进⾏⼀个⼩波滤波练习,给出计算结果,并进⾏分析答 :信号去噪是信号处理领域的⼀个经典问题,传统的去噪⽅法主要是线性滤波和⾮线性滤波,例如中值滤波和Wiener 滤波等?⼩波变换具有下列良好特性:①低熵性②多分辨率特性③去相关性④选基灵活性?⼩波在信号去噪领域已经取得越来越⼴泛的应⽤?阈值去噪的⽅法是⼀种较好的⼩波去噪法?阈值去噪⽅法的思想就是对⼩波分解后的个层系数中模⼤于和⼩于某阈值的系数进⾏处理,然后对处理完的⼩波系数再进⾏反变换,重构出经过去噪的信号?01002003004005006007008009001000-11原始信号01002003004005006007008009001000-22含噪信号01002003004005006007008009001000-202去噪后的信号。

1.总结学过的滤波器设计方法，用matlab 仿真例子分析不同设计方法的滤波器的性能及适应场合。

答：1.1模拟低通滤波器的设计方法 1.1.1 Butterworth 滤波器设计步骤： ⑴.确定阶次N① Ωc 、Ωs 和As 求Butterworth DF 阶数N② Ωc 、Ωs 和Ω=Ωp()的衰减Ap 求Butterworth DF 阶数N③ Ωp 、Ωs 和Ω=Ωp 的衰减Ap和As 求Butterworth DF 阶数N3dB p Ω≠-/10/1022(/)101,(/)101p s A A N N p c s c ΩΩ=-ΩΩ=-则：⑵.用阶次N 确定 根据公式：在左半平面的极点即为的极点，因而1.1.2 切比雪夫低通滤波器设计步骤： ⑴.确定技术指标归一化： ⑵.根据技术指标求出滤波器阶数N 及：⑶.求出归一化系统函数 其中极点由下式求出：()a H s 1,2,2N()()a a H s H s -()a H s 2,,N p Ωp αs Ωs α/1p p p λ=ΩΩ=/s s p λ=ΩΩε0.12101δε=-p δα=或者由和S 直接查表得2.数字低通滤波器的设计步骤：〔1〕 确定数字低通滤波器的技术指标：通带截止频率pω、通带最大衰减系数pα、阻带截止频率ω、阻带最小衰减系数s α。

〔2〕将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。

巴特沃斯：切比雪夫：N ()a H p /ss p λ=ΩΩ0.12101δε=-p δα=〔3〕把模拟滤波器变换成数字滤波器，即把模拟滤波器的系数)(S H 映射成数字滤波器的系统函数)(z H 。

实现系统传递函数s 域至z 域映射有脉冲响应不变法和双线性映射两种方法。

〔〕脉冲响应不变法。

按照技术要求设计一个模拟低通滤波器，得到模拟低通滤波器的传输函数()s H a 转换成数字低通滤波器的系统函数H(z)。

设模拟滤波器的传输函数为()s H a ，相应的单位冲激响应是()t h a ，()s H a =LT[()t h a ]，LT[.]代表拉氏变换，对()t h a 进展等间隔采样，采样间隔为T ，得到()nT h a ，将h(n)= ()nT h a 作为数字滤波器的单位取样响应，那么数字滤波器的系统函数H(z)便是h(n)的Z 变换。

课程大作业数字信号处理实验报告课程大作业-数字信号处理实验报告实验一信号、系统和系统响应一．实验目的1.熟悉理想采样的性质，了解信号采用前后的频谱变化，加深对采样定理的理解。

2．熟悉离散信号和系统的时域特性。

3.熟悉线性卷积的计算和编程方法：用卷积法观察和分析系统响应的时域特性。

4．掌握序列傅氏变换的计算机实现方法，利用序列的傅氏变换对离散信号、系统及系统响应进行频域分析。

二、实验原理1．连续时间信号的采样这有助于理解信号从时域到频域的变化，也有助于理解信号从时域到时域的变化。

对一个连续时间信号进行理想采样的过程可以表示为该信号和个周期冲激脉冲的乘积，即a（t）？xa（t）m（t）（1－1）x？A（T）是连续信号Xa（T）的理想采样，m（T）是周期脉冲，其中XM（T）？N（t？nt）（1－2）它也可以用傅立叶级数表示为：1.吉咪？stm（t）？？e（1－3）tn其中t为采样周期，?s?2?/t是采样角频率。

设xa(s)是连续时间信号xa(t)的双边拉氏变换，即有：xa（s）xa（t）e？stdt（1－4）a(t)的拉氏变换为此时理想采样信号x??1?(s)?x?a(t)edtxa(s?jm?s)（1－5）xa?tmst??作为拉普拉斯变换的特例，信号理想采样的傅里叶变换1xa(j?)??xa?j(??m?s)?（1－6）tm从方程（1-5）和（1-6）可以看出，信号理想采样后的频谱是原始信号频谱的周期扩展，其扩展周期等于采样频率。

根据香农采样定理，如果原始信号是带限信号，且采样频率高于原始信号最高频率分量的两倍，则采样后不会出现频率混淆。

在计算机处理时，不采用式（1－6）计算信号的频谱，而是利用序列的傅立a（t）？Xa（T）m（T），根据z变量叶变换计算信号的频谱，并定义序列x（n）？xa（新界）？根据X变换的定义，序列X（n）的Z变换可以得到：X（Z）？Nx（n）zn（1－7）以ej?代替上式中的z，就可以得到序列x(n)的傅立叶变换x(e)?j?nx(n)e???j?n（1－8）式（1-6）和式（1-8）具有以下关系：(j)x(ej)x（1－9）at由式（1－9）可知，在分析一个连续时间信号的频谱时，可以通过取样将有将相关性的计算转化为序贯傅里叶变换的计算。

现代信号处理大型作业一．试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。

滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。

（一）、分析与通常的滤波器相比，互补滤波器具有优良的结构特性和结构特性，具有较低的噪声能量和系数敏感性，其定义如下：一组滤波器H 12(),(),.......()Z H Z H Z n 如果满足下式：He Kjw k n(),==∑110<w<2π 则称这组滤波器为幅度互补滤波器;如果满足下式：He kjw k n()=∑=121, 0<w<2π则称这组滤波器为功率互补滤波器，同时互补滤波器还应该满足：Hz A z kk n()()=∑=1其中A(z)为全通函数，适当的选择全通函数，可以使两带函数具有所需要的低通和高通特性。

（二）、设计步骤(1) 对Fp 、Fr 进行预畸);();(''FsFrtg FsFptg r p ∏=Ω∏=Ω(2) 计算'''*r p c ΩΩ=Ω，判断'c Ω是否等于1，即该互补滤波器是否为互补镜像滤波器(3) 计算相关系数⎪⎩⎪⎨⎧-==+++=+-=-=ΩΩ=--=偶数)N 为(;21奇数)N 为 (;;lg /)16/1lg(;150152;1121;1;;])110)(110[(1213090500''02'''211－min1.0min1.0i i u q k N q q q q q k k q k k k k rp Ar Ap;)2cos()1(21))12(sin()1(21)1(21'2∑∑∞=∞=+-++-=Ωm mm m m m m i u Nm q u Nm q q ππ;42⎥⎦⎤⎢⎣⎡=N N;221N N N -⎥⎦⎤⎢⎣⎡=;)/1)(1(2'2'k k v i i i Ω-Ω-=12'1212,1;12N i v i i i =Ω+=--α 22'22,1;12N i v iii =Ω+=β (4) 互补镜像滤波器的数字实现;22i ii A αα+-=;22iii B ββ+-=1221,1;1)(N i ZA Z A Z H i i i =++=∏--22212,1;1)(N i ZB Z B Z Z H i i i =++=∏--- )];()([21)(21Z H Z H Z H L +=（三）、程序与结果 1． 二带滤波器组 （1） 源程序： clear; clf;Fp=1700;Fr=2300;Fs=8000; Wp=tan(pi*Fp/Fs); Wr=tan(pi*Fr/Fs); Wc=sqrt(Wp*Wr); k=Wp/Wr;k1=sqrt(sqrt(1-k^2)); q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13; N=11;N2=fix(N/4); M=fix(N/2); N1=M-N2; for jj=1:M a=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%N is odd, u=j end ab=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N); end bW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k)); endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2); endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2); endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2); endfor i=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2); endw=0:0.0001:0.5;LP=zeros(size(w));HP=zeros(size(w));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endLP(n)=abs((H1+H2)/2);HP(n)=abs((H1-H2)/2);endplot(w,LP,'b',w,HP,'r');hold on;xlabel('digital frequency');ylabel('amptitude');（2）运行结果：见图1图1 二带数字滤波器组2．四带滤波器组（1）源程序：clf;Fp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;for jj=1:Ma=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N); % N is odd, u=jendb=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));Endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfor i=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LLP=zeros(size(w));LHP=zeros(size(w));HLP=zeros(size(w));HHP=zeros(size(w));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH21=1;for i=1:N1H21=H21*(a(i)+z^(-4))/(1+a(i)*z^(-4)) ;H2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endH22=1/(z^2);for i=1:N2H22=H22*(b(i)+z^(-4))/(1+b(i)*z^(-4));endLP=((H1+H2)/2);HP=((H1-H2)/2);LLP(n)=abs((H21+H22)/2*LP);LHP(n)=abs((H21-H22)/2*LP);HHP(n)=abs((H21+H22)/2*HP);HLP(n)=abs((H21-H22)/2*HP);endplot(w,LLP,'b',w,LHP,'r',w,HLP,'k',w,HHP,'m')hold onxlabel('digital frequency');ylabel('amptitude');（2）运行结果：见图2图2 四带数字滤波器组二、根据《现代数字信号处理》第四章提供的数据，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1）Levison算法2）Burg算法3） ARMA 模型法 4） MUSIC 算法 1 Levinson 算法Levinson 算法用于求解Yule-Walker 方程，是一种按阶次进行递推的算法，即首先以AR （0）和AR （1）模型参数作为初始条件，计算AR （2）模型参数；然后根据这些参数计算AR （3）参数，等等，一直到计算出AR （p ）模型参数为止，需要的运算量数量级为2p ，其中p 为AR 模型的阶数。

２01２１1120３班ﻩﻩＡa ｒon Ｈwa ｎg ﻩﻩ学号:201２14０６１91．2 设()5cos(0.25),0,1,,15,x n n n π==为有限长序列。

(1）计算１６点DFT ，并画出幅度谱序列。

解：程序代码如下n ＝0:１5;x ＝5*cos （0.25*p ｉ＊ｎ）； figur ｅ(1); ｓｔem(ｎ,ｘ)；x ｌａbel （'ｎ＇）;ylabe ｌ（＇x(ｎ）＇)； X=fft （x)； X=ab ｓ(X ）; f ｉｇｕre （2)； ｓte ｍ(n,X);ｘlabel （'k');y ｌa ｂe ｌ(＇Ｘ(Ｋ）'）;所得图像如下 ﻩﻩﻩﻩﻩﻩﻩﻩ ﻩﻩ(2）在给序列后面补１６个零后,计算3２点DFT,并画出DF Ｔ幅度谱序列。

解:程序代码如下ｎ=0:31; n1=0:15;x1=５*c ｏｓ(0．25*ｐi*n1)； x ＝[x1 ze ｒos(1,16)］; f ｉgu ｒｅ(1);ｓt ｅm(n ，x);xlabel(＇n');y ｌa ｂe ｌ('x(n ）')； X ＝ff ｔ(ｘ); X=abs(X);fi ｇｕre （2);ｓt ｅｍ(ｎ,X);s ｔem(ｎ,Ｘｋ）;xlabel('ｋ＇)；ｙlabe ｌ(＇X(ｋ)＇);所得图像如下:ﻩﻩﻩ ﻩﻩﻩﻩﻩﻩﻩ（3）把DFT 的点数扩大为6４，然后重复(2) 解:程序代码如下n ＝０：6３； n1=０:15;x1=５*ｃos （0.２5*pi*ｎ１）; x=[x １ zer ｏs(１，４8)]； ｆi ｇｕr ｅ(１);ｓte ｍ（n ，x ）；xl ａbel （'n ＇）;yl ａbe ｌ（'ｘ(n)'); X ＝ff ｔ(x); X=ａｂｓ（X)；fi ｇu ｒe （2）；st ｅm （n,X ）;s ｔｅm （n,X);ｘla ｂe ｌ（'k');y ｌａb ｅl('X （ｋ)＇); 所得图像如下:ﻩﻩﻩﻩﻩ ﻩﻩﻩﻩﻩ(4）依据D ＴFT 与ＤFT 之间的关系,解释补零操作对DFT 的影响。

实时信号处理系统设计及实现大作业学院:班级:成员（姓名学号分工）：软件版本：Quartus II 9.0sp2 Web Edition第一题：无符号DA卷积解：原理：内积公式为假设系数c[n]为已知常数，x[n]为变量。

无符号DA系统假设变量x[n]可以表示为：其中xb[n]表示x[n]的第b位，x[n]是x的第n次采样，则内积y可以表示为：重新分布求和次序有：内积运算的实现：●用1个LUT实现映射f (c[n], xb[n])。

预先对1个2N字的LUT进行设定，输入为N位向量xb = [xb[0], xb[1], …, xb[n-1]]，输出为f (c[n], xb[n])。

●对每个LUT输出值f (c[n], xb[n]) 乘以权重2b。

●用移位-加法器实现累加运算。

定义3阶内积：设系数位宽为4，值为c[0]=2，c[1]=3，c[2]=1，c[3]=2实现f (c[n], xb[n])的LUT为：X b[3] X b[2] X b[1] X b[0] f(c[n], X b[n])0 0 0 0 2 x 0+1x0 + 3x0 + 2x0=010=000020 0 0 1 2 x 0+1x0 + 3x0 + 2x1=210=001020 0 1 0 2 x 0+1x0 + 3x1 + 2x0=310=001120 0 1 1 2 x 0+1x0 + 3x1 + 2x1=510=010120 1 0 0 2 x 0+1x1 + 3x0 + 2x0=110=000120 1 0 1 2 x 0+1x1 + 3x0 + 2x1=310=001120 1 1 0 2 x 0+1x1 + 3x1 + 2x0=410=010020 1 1 1 2 x 0+1x1 + 3x1 + 2x1=610=011021 0 0 02 x 0+1x0 + 3x0 + 2x0=210=001021 0 0 12 x 0+1x0 + 3x0 + 2x1=410=010021 0 1 02 x 0+1x0 + 3x1 + 2x0=510=010121 0 1 12 x 0+1x0 + 3x1 + 2x1=710=011121 1 0 02 x 0+1x1 + 3x0 + 2x0=310=001121 1 0 12 x 0+1x1 + 3x0 + 2x1=510=010121 1 1 02 x 0+1x1 + 3x1 + 2x0=610=010121 1 1 12 x 0+1x1 + 3x1 + 2x1=810=10002有关x[n]={ x[0]=1010=10102, x[1]=910=10012, x[2]=1210=11002, x[3]=80=10002}的内积如下：步骤t X t[3] X t[2] X t[1] X t[0] f[t]+ACC[t-1]=ACC[t]0 0 0 1 0 3x2 + 0 =31 0 0 0 1 2x21 + 3=72 0 1 0 0 1x22 + 7=113 1 1 1 1 8x23 + 11=75进行数值校验可以看到：y=<c,x>=c[0]x[0]+ c[1]x[1]+ c[2]x[2] + c[3]x[3]=2x10+3x9+1x12+2 x 8=75DA体系结构图程序代码：●顶层文件：library ieee;use ieee.std_logic_1164.ALL;use ieee.std_logic_arith.ALL;------------------------------------------------------------------------------------------------------------------- entity DAjj isport ( clk : in std_logic; --时钟reset,en: in std_logic; --复位，使能x_in0, x_in1, x_in2 ,x_in3: in std_logic_vector (3 downto 0); --四阶输入y : out integer range 0 to 200); --结果输出end DAjj;------------------------------------------------------------------------------------------------------------------- architecture con of DAjj is---------------------------------对lut表进行元件例化--------------------------------------------------- component lutbport ( table_in : in std_logic_vector (3 downto 0); --lut表输入table_out : out integer range 0 to 8); --lut表输出end component;------------------------------------------------------------------------------------------------------------------- type sta_type is (s0, s1);signal state : sta_type; --定义两个状态，用于状态机中signal x0, x1, x2,x3,table_in: std_logic_vector (3 downto 0):="0000";signal table_out : integer range 0 to 6; --定义信号与lut表输出相接begintable_in(0) <= x0(0);table_in(1) <= x1(0);table_in(2) <= x2(0);table_in(3) <= x3(0); --输入初始值process(clk,reset,en)variable m : integer range 0 to 200; --输出暂存变量variable njs : integer range 0 to 4; --计数变量beginif(reset='0') thenstate<=s0; --复位信号’0’有效elsif(en='1') then --使能信号’1’有效if(clk='1' and clk'event) then --时钟上升沿时有效case state iswhen s0 =>state <= s1;njs := 0;m := 0;x0 <= x_in0;x1 <= x_in1;x2 <= x_in2;x3 <= x_in3; --初始状态，进行初始化when s1 =>if njs = 4 theny <= m;state <= s0; --状态S1时，计数满3次时，返回初始状态（3阶）elsem := m / 2 + table_out *8; --进行移位累加运算x0(0) <= x0(1);x0(1) <= x0(2);x0(2) <= x0(3);x1(0) <= x1(1);x1(1) <= x1(2);x1(2) <= x1(3);x2(0) <= x2(1);x2(1) <= x2(2);x2(2) <= x2(3);x3(0) <= x3(1);x3(1) <= x3(2);x3(2) <= x3(3); --将下一个数送入lut表中njs := njs + 1; --计数加1，以便以后用于判断state <= s1; --进入下一个状态end if;end case;end if;end if;end process;L0: lutb port map(table_in , table_out); --进行元件例化的端口映射end con;LUT表程序：library ieee;use ieee.std_logic_1164.all;use ieee.std_logic_arith.all;------------------------------------------------------------------------------------------------------------------- entity lutb isport ( table_in : in std_logic_vector(3 downto 0);table_out : out integer range 0 to 8 );end lutb;-------------------------------------------------------------------------------------------------------------------architecture lutjs of lutb isbeginprocess (table_in)begin------------------------------------------c(0)=2,c(1)=3,c(2)=1,c(3)=2-------------------------------------------- case table_in iswhen "0000" => table_out <= 0;when "0001" => table_out <= 2;when "0010" => table_out <= 3;when "0011" => table_out <= 5;when "0100" => table_out <= 1;when "0101" => table_out <= 3;when "0110" => table_out <= 4;when "0111" => table_out <= 6;when "1000" => table_out <= 2;when "1001" => table_out <= 4;when "1010" => table_out <= 5;when "1011" => table_out <= 7;when "1100" => table_out <= 3;when "1101" => table_out <= 5;when "1110" => table_out <= 6;when "1111" => table_out <= 8;when others => table_out <= 0;------------table_out=c(0)*table_in(0)+c(1)*table_in(1)+c(2)*table_in(2)+c(3)*table_in(3)------- end case;end process;end lutjs;Quartus仿真结果由仿真图可分析的：（1）当x[n]={ x[0]=1010=10102, x[1]=910=10012, x[2]=1210=11002, x[3]=80=10002}y=<c,x>=c[0]x[0]+ c[1]x[1]+ c[2]x[2]+ c[3]x[3]=2x10+3x9+1x12+2 x8=75 与仿真结果一致；（2）当x[n]={ x[0]=510=01012, x[1]=1210=11002, x[2]=410=01002, x[3]=10=00012}y=<c,x>=c[0]x[0]+ c[1]x[1]+ c[2]x[2]+ c[3]x[3]=2x5+3x12+1x4+2 x 1=52与仿真结果一致所以仿真结果正确所用资源：使用逻辑单元为50个，组合功能14个，专用逻辑寄存器36个。

课程名称：现代信号分析考试形式：□专题研究报告□论文大作业□综合考试学生姓名：学号：评阅人：时间：年月目录一、FFT算法的特点。

......................................... - 3 -二、平稳随机过程............................................. - 8 -三、功率谱估计.............................................. - 10 -四、自适应滤波器（LMS和RLS）............................... - 16 -五、维纳滤波器.............................................. - 23 -六、FIR维纳滤波器.......................................... - 25 -七、AR/MA滤波器............................................ - 28 -八、卡尔曼滤波器............................................ - 31 -九、小波分析方法及应用...................................... - 37 -十、其他现代信号分析方法.................................... - 43 -现代信号分析- 3 -一、FFT 算法的特点。

题目：请举例说明FFT 算法的特点计算离散傅里叶变换的一种快速算法，简称FFT 。

快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。

采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N 越多，FFT 算法计算量的节省就越显著。

根据对序列分解与选取方法的不同而产生了FFT 的多种算法，基本算法是基2DIT 和基2DIF 。

现代信号处理作业2014年6月Cramer-Rao 不等式的证明定义：设),,,(21n x x x θθ=为参数θ的估计子，参数θ的估计子θ 的偏差定义为该估计子误差的期望值，即θθθθθ-=-=}{}{)(E E b def则估计子θ 称为无偏估计子。

若变差)(θ b 等于零或者θθ=}{ E ，即估计子的期望值等于真实参数，则成估计子θ 的渐进无偏估计子。

无偏性反映了参数估计量的取值在真值θ周围摆动程度。

一个参数往往具有不止一个无偏估计子，因此引入了估计量的有效性这一概念来判断偏估计的优劣，以方差大小来衡量无偏估计子的有效性。

假定参数θ存在无偏估计子1θ 和2θ ，若1θ 具有比2θ 更大的方差，即)var()var(21θθ >，则2θ 对于1θ 相对有效。

对于任意参数θ，我们自然希望得到其更为有效的无偏估计子，即要求方差值尽量小。

究竟无偏估计量的方差能够小到什么程度，是否有下界？如果有，是否有估计子能达到这一下届？这就是本次要讨论的Cramer-Rao 下界。

定义：参数θ的估计子θ 的均方误差)(2θ M 定义为该估计子与真实参数的误差平方的期望值，即)()var(}){()(222θθθθθ b E M +=-=定理：（Cramer-Rao 不等式）令),,,(21n x x x x =为样本向量。

若参数估计θ是真实参数θ的无偏估计，并且θθ∂∂)|(x f 和22)|(θθ∂∂x f 存在，则θ 的均方误差所能达到的下界成为Cramer-Rao 下界等于Fisher 信息的倒数，即)(1}){()(22θθθθJ E M ≥-= 不等式中等号成立的充分必要条件是))(()|(ln θθθθθ-=∂∂ K x f 。

其中)(θK 是θ的某个正函数，并与样本n x x x ,,,21 无关。

定义：品质函数的的方差成为Fisher 信息，用)(θJ 表示，定义为)}|(ln {})]|(ln {[)(222θθθθθx f E x f E J ∂∂-=∂∂=证明Cramer-Rao 不等式：由假设条件可得θθ=}{ E ，则可得0)|()(}{=-=-⎰∞∞-dx x f E θθθθθ（1） 对式（1）两边求偏导，得0)|()(}{=-∂∂=-∂∂⎰∞∞-dx x f E θθθθθθθ 即可得 0)|()()|(=∂∂-+-⎰⎰∞∞-∞∞-dx x f dx x f θθθθθ（2） 另，由复合函数求导法可得)|()|(ln )|(θθθθθx f x f x f ⎥⎦⎤⎢⎣⎡∂∂=∂∂ （3） 又有 ⎰∞∞-=1)|(dx x f θ （4）将式（3）和式（4）代入到式（2），得1))(|()|(ln =-⎥⎦⎤⎢⎣⎡∂∂⎰∞∞-dx x f x f θθθθθ 可改写为 []1)|()()|()|(ln =-⎥⎦⎤⎢⎣⎡∂∂⎰∞∞-dx x f x f x f θθθθθθ （5） 由Cauchy-Schwartz 不等式，可得1)()|()|()|(ln 2≥-⎥⎦⎤⎢⎣⎡∂∂⎰⎰∞∞-∞∞-dx x f dx x f x f θθθθθθ 则 ⎰⎰∞∞-∞∞-⎥⎦⎤⎢⎣⎡∂∂≥-dx x f x f dx x f )|()|(ln 1)()|(2θθθθθθ （6） 由Cauchy-Schwartz 不等式等号成立的条件可知：当且仅当)|())(()|()|(ln θθθθθθθx f K x f x f -=∂∂ ，不等式可取等号。

2012111203班AaronHwang 学号：20121406191.2 设 x(n) =5cos(0.25二n), n =0,1^|,15,为有限长序列。

(1) 计算16点DFT并画出幅度谱序列解：程序代码如下n=0:15;x=5*cos(0.25*pi* n);figure(1);stem( n, x);xlabel(' n');ylabel('x( n)');X=fft(x);X=abs(X);figure(2);stem( n,X); xlabel('k');ylabel('X(K)');所得图像如下(2) 在给序列后面补16个零后，计算32点DFT并画出DFT畐度谱序列解：程序代码如下n=0:31;n1=0:15;x1=5*cos(0.25*pi* n1);x=[x1 zeros(1,16)];figure(1);stem( n,x);xlabel(' n');ylabel('x (n)');X=fft(x);X=abs(X);figure (2) ;stem( n,X);stem( n,Xk);xlabel('k');ylabel('X(k)');所得图像如下：(4)依据DTFT 与 DFT 之间的关系，解释补零操作对 DFT 的影响。

解：通过补零操作和DFT 点数的增加可以更加准确的得出原信号的 DFS 图像。

1.5 ( 1)设H 1(z) = 1 - r o e -1扯z ，为单零点系统的传递函数，求其幅度响应HMe^)的解析式。

令r 。

=0.95「0 =0.25二，画出幅度响应的草图。

解：由已知H i (z) = 1 - r °e，可知响应为|H i (e j)|=、1 r 。

22r °cos( s •)(3)把DFT 的点数扩大为64,然后重复(2) 解：程序代码如下 n=0:63; n1=0:15; x1=5*cos(0.25*pi* n1) x=[x1 zeros(1,48)]; figure(1);stem( n,x);xlabel(' n');ylabel('x (n)'); X=fft(x); X=abs(X);figure (2) ;stem( n,X);stem( n, X);xlabel('k');ylabel('X(k)');所得图像如下：程序代码如下r0=0.95;w0=0.25*pi;w=li nspace(0,2*pi,100);H=1-rO*exp(-j*wO)*exp(-j*w);plot(w,abs(H));title('H (jw)');所得图像如下：(2)设H2(Z)= --------- ［竹为单极点系统的传递函数，求其幅度响应H2(e尬)的解析1 -r0e式。

第五章计算机仿真大作业采用计算机编程实现图1中的自适应均衡器：()h n ()s n ()x n ⊕()noise n ()y n 自适应均衡器()z n图1 信号传输的系统模型图1中()s n 为频率为10Hz 、采样频率为1000Hz 的正弦序列，假设该信号通过一个具有码间干扰特性的信道，其单位抽样响应为()[0.005,0.009,0.024,h n =--0.854,0.218,0.049,0.0323]--，经过上述信道的输出信号()x n 与高斯白噪声()noise n 叠加后作为自适应均衡器的输入信号()y n ，()z n 为自适应均衡器的输出信号。

其中图1中所示的自适应均衡器为一N=31阶FIR 滤波器，抽头系数为(),0,1,,1in i N ω=-其结构如图2所示： 1z -()y n 0ω1ω∑∑++()()d n s n =+-()z n ()e n 1z -1N ω-自适应算法图2自适应均衡器结构图按要求分析回答下列问题，并给出分析结果和波形图： 1.在一个图中用子图的形式（subplot ）画出图1中： (1)()s n 信号；(2)()s n 经信道()h n 传输后的()x n 信号；(3)当()x n 加()noise n 的信噪比SNR(dB ）为20dB 时均衡器的输入序列()y n 的波形图； 对上述波形进行对照分析和说明。

01002003004005006007008009001000-101正弦信号s(n)01002003004005006007008009001000-101x(n)序列1002003004005006007008009001000-101y(n)序列分析说明：s(n)通过具有码间干扰特性的信道h(n)，由于信道存在一定的误差和码间干扰使系统的性能下降，x(n)的波形密度减小了，但整体波形没有发生变化。

加入噪声后，y （n ）的幅值没有变化，但整个波形由于受到噪声干扰浮现“毛刺”现象，波形不在平滑。

1. Levinson-Durbin 算法1.1 Levinson-Durbin 算法P 阶AR 模型的差分方程为：1x()()()pi i n a x n i w n =+-=∑，其中()w n 是均值为0的白噪声。

AR 过程的线性预测方法为：先求得观测数据的自相关函数，然后利用Yule-Walker 方程递推求得模型参数，再根据公式求得功率谱的估计。

LD 递推算法是一种采用AR 模型的现代谱估计方法。

LD 递推算法具体计算步骤如下：（1） Yule-Walker 方程的矩阵形式如下所示：2,1,1(0)(1)(2)()(1)(0)(1)(1)0()(1)(2)(0)0xx xx xx xx k xx xx xx x k k xx xx xx xx r r r r k a r r r r k a r k r k r k r σ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥--+⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦系数矩阵Hxx xx R R =为Hermitian 矩阵，其中k=1,2……p, (i=1…k)表示k 阶时的预测系数，通过归纳递推求解出p 阶的预测系数。

（2） 当p=1时，即一阶递推为：⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-01)0()1()1()0(211,1σa R R R R x x x x求解可得：1,01,1211,1(1)1, (0)(0)(1)x x x x R a a R R a R σ==-=+（3） 当2≥p 时，递推公式如下：10,≡p a , *,1,1,kp p p k p k p a K a a ---+=, ]1[2212p p p K -=-σσ 21,-∆-==p pp p p a K σ, ∑-=--+=∆11,1)()(p k x kp x p k p R ap R矩阵x R 已知，可得到各阶AR 模型系数为：)0())1(1( ,)0()1()1(2111xx xx xx r a r r a -=-=ρ11111)()()()(--=--∑-+-=∆-=k k l xx k xx k kk l k r l a k r k a ρρ1,,2,1)()()()(*11-=-+=--k i i k a k a i a i a k k k k12))(1(--=k k k k a ρρ1.2 MATLAB 算法实现function [a_p var_p] = Levinson_Durbin(x,p) N = length(x); for ii=1:NRxx(ii) = x(1:N-ii+1)*(x(ii:N))'/N; end a(1)=1;a(2)=-Rxx(2)/Rxx(1);for k=1:p-1 % Levinson-Durbin algorithm var(k+1) = Rxx(0+1)+a(1+1:k+1)*Rxx(1+1:k+1)';reflect_coefficient(k+1+1) = -a(0+1:k+1)*(fliplr(Rxx(2:k+1+1)))'/var(k+1); var(k+1+1) = (1-(reflect_coefficient(k+1+1))^2)*var(k+1); a_temp(1) = 1; for kk=1:ka_temp(kk+1) = a(kk+1)+reflect_coefficient(k+1+1)*a(k+1-kk+1); enda_temp(k+1+1) = reflect_coefficient(k+1+1); a = a_temp; enda_p = a; % prediction coeffecients var_p = var(p+1); % prediction error power1.3 仿真结果（1） 当=2p 时，仿真结果如下：图1 p=2时仿真结果图预测系数：222[(0),(1),(2)]a a a =[1,-0.5236,0.5094]； 误差功率：0.9790（2） 当=10p 时，仿真结果如下：图2 p=10时仿真结果图预测系数：101010[(0),(1),(2),]a a a = [1,-0.4956,0.5058,-0.0057,0.0373,-0.0161,0.0276,0.0285,-0.0221,0.0212,-0.0179] 误差功率：0.9963（3） 当=50p 时，仿真结果如下：图3 p=50时仿真结果图预测系数：50505050[(0),(1),(2),,(50)]a a a a =[1,-0.4958,0.5031,…,0.0078,0.0191]误差功率：0.96761.4 结果分析从p=2的仿真结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；p=10时，采用LD 算法进行估计后，得到的功率谱产生振荡；p=50时，得到的功率谱产生较大的振荡。

现代信号处理作业班级:姓名：流水号：1、设采样周期T=250μs （采样频率fs =4kHz ），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB 边界频率为fc =1kHz 。

答：MATLAB 程序如下：[B,A]=butter(3,2*pi*1000,'s'); [num1,den1]=impinvar(B,A,4000); [h1,w]=freqz(num1,den1); [B,A]=butter(3,2/0.00025,'s'); [num2,den2]=bilinear(B,A,4000); [h2,w]=freqz(num2,den2); f=w/pi*2000;plot(f,abs(h1),'-.',f,abs(h2),'-'); grid;title(‘现代信号处理第一题三阶巴特沃兹滤波器’); xlabel('频率/Hz ’); ylabel('幅值/dB'); 程序运行结果如下：20040060080010001200140016001800200000.10.20.30.40.50.60.70.80.91现代信号处理第一题三阶巴特沃兹滤波器频率/Hz幅值/d B程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。

脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。

同时，也看到双线性变换法，在z=-1即Ω=π或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在ω=∞处的三阶传输零点通过映射形成的。

现代信号处理 大作业-1大作业-1描述：求下列假设的最佳判决方式，并求解每种错误的概率。

（提示，可以假设K 非常大，并利用中心极限定理将大量同分布随机变量之和近似为高斯。

）200211222:[](0,):[](0,)=1:[](0,)H r k N H r k N k K H r k N σσσ::L : ，，解：根据题意2022210011(|)exp [](2)2KK k p r H r k πσσ=⎡⎤=-⎢⎥⎣⎦∑ 2122211111(|)exp [](2)2K K k p r H r k πσσ=⎡⎤=-⎢⎥⎣⎦∑2222212211(|)exp [](2)2K K k p r H r k πσσ=⎡⎤=-⎢⎥⎣⎦∑ 假设判决为0H ，则应满足以下条件：00110011(|)()(|)()(|)()(|)()p r H p H p r H p H p r H p H p r H p H ≥⎧⎨≥⎩假设先验概率为0121()=()=()=3p H p H p H ，并假设0120σσσ≥≥>，代入上式得：212210012222100211exp []2211exp []22K K k KK k r k r k σσσσσσσσ==⎧⎡⎤⎛⎫⎛⎫⎪≥-⎢⎥ ⎪ ⎪⎪⎝⎭⎝⎭⎣⎦⎨⎡⎤⎛⎫⎛⎫⎪≥-⎢⎥⎪ ⎪⎪⎝⎭⎝⎭⎣⎦⎩∑∑ 进一步化简可得到：2101220122012202(ln ln )[]1122(ln ln )[]1122K k Kk K r k K r k σσσσσσσσ==-⎧≥⎪-⎪⎪⎨-⎪≥⎪-⎪⎩∑∑即 21020122220102(ln ln )(ln ln )[]max ,11112222Kk K K r k σσσσσσσσ=⎛⎫ ⎪-- ⎪≥ ⎪-- ⎪⎝⎭∑。

根据上述同样的思路，可以得到判决为1H ，需要满足：21021122221201(ln ln )(ln ln )[]11112222K k K K r k σσσσσσσσ=--≤≤--∑判决为2H ，需要满足：22021122220212(ln ln )(ln ln )[]min ,11112222Kk K K r k σσσσσσσσ=⎛⎫⎪-- ⎪≥ ⎪-- ⎪⎝⎭∑ 下面为了易于表述，记211[]()Kk r k f r K==∑。