现代信号处理作业
现代信号处理大作业王成志1
《现代信号处理》大作业姓名：王成志学号：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.做⼀个⽹络检索,简述现代信号处理技术的主要特征和技术特点,并阐述信号处理在实际⼯程中的应⽤情况代信号处理技术的主要特征和技术特点: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去噪后的信号。
现代信号处理仿真作业一(318谐波恢复)
现代信号处理仿真作业一(3.18谐波恢复)班级：自研42 姓名：李琳琳 学号：2004211068一． 谐波恢复的基本理论和方法：1． Pisarenko 谐波分解理论谐波过程可用差分方程描述，首先利用Pisarenko 谐波分解理论推导谐波过程所对应的差分方程。
对单个正弦波()sin(2)x n fn πθ=+利用三角函数恒等式，有：()2cos(2)(1)(2)0x n f x n x n π--+-=对上式作z 变换，得： 12[12cos(2)]()0f z z X z π---+=得到特征多项式：1212cos(2)0f z z π---+=。
由此可见，正弦波的频率可以由相应特征方程的一对共轭根来决定：|arctan[Im()/Re()]|/2i i i f z z π=将单个正弦波推广到多个正弦波的情形，得：如果p 个实的正弦波信号没有重复频率的话，则这p 个频率应该由特征多项式1(21)212110p p p a z a z z -----++++= (1)的根决定。
由此可得到p 个实正弦波所组成的谐波过程可以用以下的差分方程进行描述：21()()0pi i x n a xn i =+-=∑这是一个无激励的AR 过程。
2. 谐波恢复的ARMA 建模法在无激励的AR 模型差分方程21()()0pi i x n a x n i =+-=∑两边同乘()x n k -，并取数学期望，则有：21()()0,px i x i R k a R k i k =+-=∀∑ (2)正弦波过程一般是在加性白噪声中被观测的，设加性白噪声为()w n ，即观测过程为： 1()()()sin(2)()pi i i i y n x n w n A f n w n πθ==+=++∑，其中()w n 为0均值的高斯白噪声。
由于谐波过程和加性白噪声统计独立，所以有：2()()()()()y x w x w R k R k R k R k k σδ=+=+ (3)将(3)式代入(2)式得： 21()()0,2py i y i R k a R k i k p =+-=>∑..................... . (4)这就是加性白噪声中观测到的p 个正弦信号所组成的谐波信号的ARMA 过程的所服从的法方程。
现代信号处理作业
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 变换。
现代信号处理大作业
现代信号处理大型作业一．试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：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 模型的阶数。
现代信号处理例题及matlab代码实现
《现代信号处理》期末考核作业1 MATLAB仿真均值为0，方差为1的白噪声信号，信号长度N=1024，并用周期图法分别求500、1000和1500次实现的平均功率谱密度，画图。
程序代码如下：clear;clear all;N=1024;%数据长度Nfft=1024;%FFT所采用的数据长度n=0:N-1;wn=randn(1,N);%产生随机白噪声subplot(2,2,1);%绘出白噪声序列plot(n,wn);title('白噪声');%500次实现的平均功率谱密度s=zeros(1,N);for i=1:500wn=randn(1,N);%产生随机白噪声Pxx=10*log10(abs(fft(wn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转换为dbs=s+Pxx;ends=s/500;f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列subplot(222);plot(f,s);xlabel('频率/Hz');ylabel('功率谱/dB');title('500次实现的平均功率谱密度');grid on;%1000次实现的平均功率谱密度s=zeros(1,N);for i=1:1000wn=randn(1,N);%产生随机白噪声Pxx=10*log10(abs(fft(wn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转换为dbs=s+Pxx;ends=s/1000;f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列subplot(223);plot(f,s);xlabel('频率/Hz');ylabel('功率谱/dB');title('1000次实现的平均功率谱密度');grid on;%500次实现的平均功率谱密度s=zeros(1,N);for i=1:1500wn=randn(1,N);%产生随机白噪声Pxx=10*log10(abs(fft(wn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转换为dbs=s+Pxx;ends=s/1500;f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列subplot(224);plot(f,s);xlabel('频率/Hz');ylabel('功率谱/dB');title('1500次实现的平均功率谱密度');grid on;实验结果图如下：2仿真如下随机过程：sin(n+)sin(n+)23n n x V ππ=++φφ其中：V n 是均值为0，方差为1的Gaussian 白噪声过程，Φ为随机相位，在[0，2π]间服从均匀分布。
现代数字信号处理仿真作业
第三章仿真作业3.17（1）代码clear;N=32;m=[-N+1:N-1];noise=(randn(1,N)+j*randn(1,N))/sqrt(2); f1=0.15;f2=0.17;f3=0.26;SNR1=30;SNR2=30;SNR3=27;A1=10^(SNR1/20);A2=10^(SNR2/20);A3=10^(SNR3/20);signal1=A1*exp(j*2*pi*f1*(0:N-1));signal2=A2*exp(j*2*pi*f2*(0:N-1));signal3=A3*exp(j*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;uk=fft(un,2*N);sk=(1/N) *abs(uk).^2;r0=ifft(sk);r1=[r0(N+2:2*N),r0(1:N)];r=xcorr(un,N-1,'biased');figuresubplot(2,2,1)stem(m,real(r1));xlabel('m');ylabel('FFT估计r1实部');subplot(2,2,2)stem(m,imag(r1));xlabel('m');ylabel('FFT估计r1虚部');subplot(2,2,3)stem(m,real(r));xlabel('m');ylabel('平均估计r实部');subplot(2,2,4)stem(m,imag(r));xlabel('m');ylabel('平均估计r虚部');仿真结果（2）代码 clear; N=256;noise=(randn(1,N)+j*randn(1,N))/sqrt(2); f1=0.15; f2=0.17; f3=0.26; SNR1=30; SNR2=30; SNR3=27;A1=10^(SNR1/20); A2=10^(SNR2/20); A3=10^(SNR3/20);signal1=A1*exp(j*2*pi*f1*(0:N-1)); signal2=A2*exp(j*2*pi*f2*(0:N-1)); signal3=A3*exp(j*2*pi*f3*(0:N-1)); un=signal1+signal2+signal3+noise;-40-2002040-200020004000mF F T 估计r 1实部-40-2002040-2000-1000010002000mF F T 估计r 1虚部-40-2002040-200020004000m平均估计r 实部-40-2002040-2000-1000010002000m平均估计r 虚部spr=fftshift((1/NF)*abs(fft(un,NF)).^2);f1=(0:length(spr)-1)*(1/(length(spr)-1))-0.5; M=64;r=xcorr(un,M,'biased');bt=fftshift(abs(fft(r,NF)));f2=(0:length(bt)-1)*(1/(length(bt)-1))-0.5; figuresubplot(1,2,1)plot(f1,10*log10(spr/max(spr))); xlabel('w/2pi'); 仿真结果（3） 代码clear; N=1000;fai1=rand(1,1)*2*pi; fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;cx=xcorr(un,p,'biased'); rxx=cx(p+1:2*p)'; R=toeplitz(rxx); [u,s]=eig(R);w/2piw/2piww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果-4-0.5-0.4-0.3-0.2-0.100.10.20.30.40.53.20（1）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise;p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);nw=128;ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果距离单位圆最近的两个解为-0.2363-0.9717i和0.3747+0.9271i，对应的归一化频率为0.1889和-0.2880（2）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);nw=128;ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果3.21-2-1123456-3clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise; p=8;for k=1:N-pxs(:,k)=un(k+p-1:-1:k)';endrxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-p-1);rxy=xs(:,1:end-1)*xs(:,2:end)'/(N-p-1);[u,e]=svd(rxx);ev=diag(e);emin=ev(end);z=[zeros(p-1,1),eye(p-1);0,zeros(1,p-1)];cxx=rxx-emin*eye(p);cxy=rxy-emin*z;[U,E]=eig(cxx,cxy);Z=diag(E);ph=angle(Z)/(2*pi);err=abs(abs(Z)-1);仿真结果最接近单位圆的两个解分别为0.5867+0.8097i和0.0349-0.9984i，对应的归一化频率为0.1502和-0.2444。
现代信号处理期末作业
现代信号处理期末作业
1.假设一物体围绕一圆周运动，角加速度是高斯白噪声，观测噪声也是零均值
的白噪声，通过编程实现使用卡尔曼滤波器实现运动物体的跟踪
设计思路：卡尔曼滤波是采用递推的算法实现的，其基本思想是先不考虑输入信号ωk和观测噪声vk的影响，得到状态变量和输出信号(观测数据)的估计值，再用输出信号的估计误差加权后校正状态变量的估计值，使状态变量估计误差的均方值最小。
这个题目设计时，就先给出题目中给出的高斯白噪声的产生，然后通过状态方程运用，从而计算输出，然后进行校正。
此图为卡尔曼滤波器输出值与观测值显示，从图中可以看到输出值与观测值基本一致。
跟踪的实现，通过一个循环，通过for循环，不断画出卡尔曼滤波器输出值和观测值，从而实现跟踪。
图中圆周上两个点会随着时间不断变化，实现跟踪。
2.设计一个自适应滤波器，进行未知系统识别。
设计思路：先假定一未知滤波器，然后可以得到未知滤波器的输出，然后通过将参数带入自适应滤波器，通过自适应滤波器，然后将未知系统和自适应输出的时域冲击响应进行对比，然后对扶贫响应进行对比，从而识别未知系统
3.设计一个自适应滤波器，对x(n)进行分离。
只需要产生一个信号并加入噪声，通过自适应滤波器，输出信号即为预测信号。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
信号时频分析技术及matlab仿真
电路与系统王冠军 201128013926153
摘要:本文介绍了时频分析的一些基础理论,对短时傅里叶变换Wigner-Ville分布做了简单介绍，运用MATLAB语言实现了旨在构造一种时间和频率的密度函数,以揭示信号中所包含的频率分量及其演化特性的wigner-ville分布。
并对时频分析方法的优缺点进行了分析。
关键词:时频分析短时傅里叶变换wigner-ville分布
1 引言
基于Fourier变换的传统信号处理技术从信号频域表示及能量的频域分布的角度揭示了信号在频域的特征。
但Fourier变换是一种整体变换，只能为人们提供信号在时域或频域的全局特性而无法了解信号频谱随时间变化的情况。
因此，需要使用一种时间和频率的联合函数来表示信号，这种表示简称为信号，也就是信号的时频分析。
2 时频分析方法
信号时频分析主要研究非平稳信号或时变信号的频谱含量是怎样随时间变化的。
时频分析是当今信号处理领域的一个主要研究热点，目前常用时频分析方法主要有短时傅里叶变换、Gabor展开、小波变换、Wigner-Ville分布。
本文主要介绍了短时傅里叶变换和Wigner-Ville分布两种分析方法。
2．1 短时傅立叶变换STFT
从历史上看,信号的时频分析用的最多的是短时傅立叶变换,这种变换的基本思想是用一个窗函数乘时间信号,该窗函数的时宽足够窄,使取出的信号可以被看成是平稳的,然后进行的傅立叶变换可以反映该时宽中的频谱变化规律,如果让窗随时间轴移动,可以得到信号频谱随时间变化的规律。
对于时变信号,了解不同时刻附近的频域特征是至关重要的。
因此,人们采用时间—频率描述时变信号,将一维的时域信号映射到一个二维的时域平面,全面反映观测信号的时频联合特征。
短时傅立叶变换反映了这一思想,对于时变信号,采用某一滑动窗函数截取信号,并认为这些信号是准平稳的,然后,再分别对其进行傅立叶变换,构成时变信号的时变谱。
短时傅立叶变换是一种常用的时—频域分析方法,其基本思想
是在傅立叶变换的基础上实现时域的局部化。
由于傅立叶变换在时域和频域的对偶关系,所以短时傅立叶变换可从时域、频域来描述,相应的短时傅立叶变换可以从傅立叶变换以及频域滤波的观点来考虑。
短时傅立叶变换的定义(1)
(),()()jwn X m STFT n w x m w n m e ∞-=-∞=
-∑
式中:w ( n) 是一个窗函数,其作用是取出;在x ( n) 某时刻附近的一小段信号进行傅立叶变换,当n 变化时,窗函数随n 移动,从而得到信号频谱随时间n 变化的规律,此时的傅立叶变换是一个二维域( n , w) 的函数。
2.2 Wigner-Ville 分布
Wigner-Ville 分布即当核函数1),(=v τφ的Cohen 类时频分布：
(,)1(,)()v f v f ϕτδ=⎧⎨ψ=⎩
核函数 信号项：)()(),(21f f f f f t P auto -+-=δδ
交叉项： 计算公式为：
ττττπd e t z t z f t W f j z 2*)2()2(),(-∞
∞--+=⎰ 其中()z t 为解析信号
Wigner-Ville 分布在分析信号时具有很好的能量集中特性,但存在的问题是由于它是一种非线性变换,因此它在分析多分量信号时将产生交叉干扰项，这就影响了WD 作为信号能量密度的物理解释,交叉项出现在自主项中间,且呈振荡状态。
由于交叉项的存在,势必会影响到时频表示的分辨率。
交叉项是二次型时－频分布的固有结果，是多分量信号中不同信号分量之间的交叉作用的结果，信号项产生于信号的每个分量本身，它与信号所表现的物理特性是一致的；而交叉项则相反，它是一种干扰产物，它所表现的是与原信号物理特性相矛盾的。
另外，时－频分布的交叉项一般是比较严重，在大多数情况下也是有害的，如何抑制交叉项也就成了设计和使用时－频分布时一个需要认真考虑的问题。
3 wigner - ville 分布的MATLAB 设计
用wigner - ville 分布分析解析信号2j k t e π ,技术参数: k = 6 ,0 ≤T ≤4 ;信号带宽: f c = k T = 24 ;采样频率:f S = 4 f c ;采样点数N = T f S = 600 ,信号加了复高
])(2cos[)2(2),(1212t f f f f f f t P cross -+-
=πδ
斯白噪声,用WD 仿真并分析其特性。
设计的源程序简要如下:
fc = k * T ;
f s = 4 * fc ;
Ts = 1/ f s ;
N = T/ Ts ;
x = zeros (1 ,N) ;
t = 0 :N - 1 ;
x = exp (j * k * pi * (t * Ts) ^2) ;
subplot (2 ,2 ,1) ;
plot (t 3 Ts ,real (x) ) ;
X = ff t (x) ;
X = ff t shif t (X) ;
subplot (2 ,2 ,2) ;
plot ( (t - N/ 2) * f s/ N ,abs (X) ) ;
Nw = 20 ;
L = Nw/ 2 ;
Tn = (N - Nw) / L + 1 ;
nff t = 32 ;
TF = zeros ( Tn ,nff t) ;
for i = 1 : Tn
xw = x ( (i - 1) *10 + 1 ；i *10 + 10) ;
temp = ff t (xw ,nff t) ;
temp = ff t shif t (temp) ;
TF(i , :) = temp ;
end
subplot (2 ,2 ,3) ;
fnew = ( (1 :nff t) - nff t/ 2) *f s/ nfft ;
tnew = (1 : Tn) *L *Ts ;
[ F ,T] = meshgrid (fnew ,tnew) ;
mesh ( F ,T ,abs ( TF) ) ;
subplot (2 ,2 ,4) ;
contour ( F ,T ,abs ( TF) ) ;
可以看出,由于实信号的频谱包括信号中的每个物理频率对应的正负频率部分,致使其所含信号分量为有物理意义的信号分量的两倍。
如果直接采用实信号进行时频分析,我们将得到包括所有正负频率分量信号项的干扰。
当有噪声存在时,时频谱变差。
在时间相干性不强的目标及海量信息将被消弱,当信号环境更加复杂时,干扰较强的情况下,真个频带内部都会有射频干扰形成的谱峰存在,难于设置相应的频域滤波器,而且也会造成其他有用信息过多的损失。
4 时频分析的优点和缺陷的讨论
时频分析以联合时频分布的形式来表示信号的特性，具有很多的优点: 1）它克服了傅里叶分析时域和频域完全分离的缺陷，将时频两域联合起来对信号进行分析，能同时考虑到两个方面的性能。
2）弥补了信号的时间能量密度和频谱能量密度不能充分描述信号的物理特性的缺陷。
3）在时频相平面上，可以精确地定位在某一时刻出现了哪些频率分量，以及某一分量出现在哪些时刻。
4）对于不同情况的信号，通过对核函数施加一些约束条件，就可以设计符合期望性能的时频分布，来满足处理信号的要求。
5）为非平稳信号的分析提供了有效的工具，为信号的分析开辟了新的途径。
时频分析虽然具有很多优点，但同时也具有不少缺点:
1）由于双线性形式的时频分布是非线性的，使得两个信号和的时频分布已不再是两个信号各自分布的和，即存在交叉项。
2）在时频域进行去噪时计算复杂，不易实现。
3）时频分布虽然反映了信号的能量分布，但不能用信号的“瞬时能量”来解释某一时刻或某一频率处的时频分布。
参考文献：
[1]李敬勇. 对线性调频脉冲压缩雷达干扰的时频分析. 1998 年,第13卷,第三期
[2]纪跃波,秦树人,汤宝平. Wigner 分布干扰项抑制及其算法. 重庆大学学报,2001 年第4 期
[3]孙晓兵,保铮,罗琳. 时频信号分析和雷达的多目标分辨. 系统工程与电子技术,1997 年(1)
[4]张善文,甄蜀春,赵兴录. 编队目标架次的识别. 航天电子对抗,2001(1)
[5]张善文,甄蜀春,赵兴录. 强线性调频干扰下雷达目标的识别方法.航天电子对抗2000 (4)
[6]胡国胜,张国红. 一种新的信号处理方法———线性调频小波变换. 数学的实践与认识,2003 年2 月
[7]孙孔峰,皇甫堪. 国防科技大学学报. 1997 年6 月Vol . 19No. 38 丁玉美,阔永红,高新波. 数字信号处理. 2002
[9]张贤达. 现代信号处理. 清华大学出版社,1994
[10]侯慧群. 雷达对抗原理. 空军第二航空学院,2001。