现代信号处理大作业题目+答案
现代信号处理大作业王成志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、已知0()2cos(2)a x t f t π=式中0f =100HZ,以采样频率s f =400Hz 对()a x t 进行采样，得到采样信号ˆ()a xt 和时域离散信号()x n ，试完成下面各题： （1）写出()a x t 的傅里叶变换表示式()a X j Ω； （2）写出()a x t 和()x n 的表达式；（3）分别求出()a x t 的傅里叶变换和()x n 的傅里叶变换。
解：（1）000()()2cos()()j tj ta a j t j t j t X j x t edt t edte e e dt∞∞-Ω-Ω-∞-∞∞Ω-Ω-Ω-∞Ω==Ω=+⎰⎰⎰上式中指数函数和傅里叶变换不存在，引入奇异函数δ函数，它的傅里叶变换可以表示成：00()2[()()]a X j πδδΩ=Ω-Ω+Ω+Ω（2）0ˆ()()()2cos()()()2cos(),a an n xt x t t nT nT t nT x n nT n δδ∞∞=-∞=-∞=-=Ω-=Ω-∞<<∞∑∑2、用微处理器对实数序列作谱分析，要求谱分辨率50F Hz ≤，信号最高频率1KHz,是确定以下各参数： （1）最小记录时间min p T （2）最大取样时间max T （3）最少采样点数min N（4）在频带宽度不变的情况下将频率分辨率提高一倍的N 值。
解：（1）已知50F Hz ≤min 110.0250p T s F === （2） max 3min max 1110.52210s T ms f f ====⨯ （3） min 30.02400.510p T s N T s-===⨯ （4）频带宽度不变就意味着采样间隔T 不变，应该使记录时间扩大一倍为0.04s 实频率分辩率提高1倍（F 变成原来的12）min 30.04800.510p T s N T s-===⨯3、在时域对一有限长的模拟信号以4KHZ 采样，然后对采到的N 个抽样做N 点DFT ，所得离散谱线的间距相当于模拟频率100HZ 。
信号处理原理与应用答案
信号处理原理与应用答案【篇一：《现代信号处理理论与应用》作业题】txt>（2）纸质版（手写），认真，不得抄袭。
《现代信号处理理论与应用》作业题（1）推导《现代信号处理》（张贤达版）p37页公式（2-5-8）。
（2）推导矢量参数估计的cramer-rao不等式，并讨论等号成立条件。
（3）令观测样本为xi?s??i，(i?1,?,n)其中{?i}是一高斯白噪声，其均值为零，方差为1。
证明：s的极大似然估计是无偏的和一致的。
（4）若信号满足：s1(t)?a1co?st?a2co2s?t???apcosp?ts2(t)?b1sin?t?b2sin2?t???bpsinp?t观测信号为x(t)?s1(t)?s2(t)?n(t)，n(t)是均值为0，均方差为1的高斯白噪声。
计算a1,a2,?,ap,b1,b2,?,bp的最小二乘估计。
(5) 输入信号x(t)为高斯-马尔可夫信号s(t)和噪声n(t)的叠加，信号和噪声假定不相关，其功率谱分别为ss(?)?3和sn(?)?1。
计算?分别取0，+1和-1时，物理不可实现21??维纳滤波器的冲激响应和最小均方误差。
(6) 设信号满足的状态方程及观测方程分别为：?11?xk?1???xk?uk 01 ??yk?1?[1,0]xk?1?nk?1{uk,k?0}和{nk?1,k?0}是均值为零的高斯白噪且与初始状态x0独立，并有?50??10?，var{n}?r?2?(?1)k?1，初始状态的方程矩阵。
vx0??var{uk}?q??k?1k?1???010??01?计算卡尔曼滤波增益m(k)。
(7)接收信号为x(t)?s(t)?n(t)，其中n(t)是功率谱密度为n0的高斯白噪。
信号为 2t???2s(t)??2e,t?0。
??0,t?0求匹配滤波器传输函数及其脉冲响应，该匹配滤波器是物理可实现吗？有无可能将它变为物理可实现？若可能，求出滤波器的传递函数及其输出信噪比最大时刻，若不可能，说明理由。
现代信号处理试题(习题教学)
1、已知0()2cos(2)a x t f t π=式中0f =100HZ,以采样频率s f =400Hz 对()a x t 进行采样，得到采样信号ˆ()a xt 和时域离散信号()x n ，试完成下面各题： （1）写出()a x t 的傅里叶变换表示式()a X j Ω；（2）写出()a x t 和()x n 的表达式；（3）分别求出()a x t 的傅里叶变换和()x n 的傅里叶变换。
解：（1）000()()2cos()()j t j t a a j t j t j t X j x t e dt t e dt e e e dt ∞∞-Ω-Ω-∞-∞∞Ω-Ω-Ω-∞Ω==Ω=+⎰⎰⎰上式中指数函数和傅里叶变换不存在，引入奇异函数δ函数，它的傅里叶变换可以表示成：00()2[()()]a X j πδδΩ=Ω-Ω+Ω+Ω（2）00ˆ()()()2cos()()()2cos(),a a n n xt x t t nT nT t nT x n nT n δδ∞∞=-∞=-∞=-=Ω-=Ω-∞<<∞∑∑2、用微处理器对实数序列作谱分析，要求谱分辨率50F Hz ≤，信号最高频率1KHz,是确定以下各参数：（1）最小记录时间min p T（2）最大取样时间max T（3）最少采样点数min N（4）在频带宽度不变的情况下将频率分辨率提高一倍的N 值。
解：（1）已知50F Hz ≤min 110.0250p T s F === （2） max 3min max 1110.52210s T ms f f ====⨯ （3） min 30.02400.510p T s N T s-===⨯ （4）频带宽度不变就意味着采样间隔T 不变，应该使记录时间扩大一倍为0.04s 实频率分辩率提高1倍（F 变成原来的12）min 30.04800.510p T s N T s -===⨯ 3、在时域对一有限长的模拟信号以4KHZ 采样，然后对采到的N 个抽样做N 点DFT ，所得离散谱线的间距相当于模拟频率100HZ 。
现代信号处理大型作业题目+答案
研究生“现代信号处理”课程大型作业（以下四个题目任选三题做）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、求证：,()(,)x i j x i j xi xj R t t C t t m m =+。
证明：(,)(,)(,,,)x i j i j iji j i j i j R t t E x x x xp x x t t dx dx ==⎰⎰(,)[(),()](),()(,,,)()(,,,)(,)(,)i j i j j i i j i j j i i j i jx i j i x j x i x j x i j i j i ji j i x j x x x i j i j i j x i j x x x x x x x i j x x C t t E x m x m x m x m p x x t t dx dx x x x m x m m m p x x t t dx dx R t t m m m m m m R t t m m =--=--=--+=--+=-⎰⎰⎰⎰ 2、令()x n 和()y n 不是相关的随机信号，试证：若()()()w n x n y n =+，则w x y m m m=+和222w x y σσσ=+。
证明：(1)[()][()()][()][()]x ym E n E x n y n E x n E y n m m ωω==+=+=+ (2)2222222222[(())]{[()()()]}[(())(())][(())][(())]2[(())(())]2[]x y x y x y x y x y x y x y x y x y x yE n m E x n y n m m E x n m y n m E x n m E y n m E x n m y n m m m m m m m m m ωωσωσσσσ=-=+-+=-+-=-+-+--=++--+=+即222x y ωσσσ=+3、试证明平稳随机信号自相关函数的极限性质，即证明： ①当0τ=时，2(0),(0)x x x x R D C σ==； ②当τ=∞时,2(),()0x x x R m C ∞=∞=。
现代信号处理大作业
姓名：潘晓丹 学号：班级：A1403492作业1LD 算法实现AR 过程估计1.1 AR 模型p 阶AR 模型的差分方程为：)()()(1n w i n x a n x pii =-+∑=，其中)(n w 是均值为0的白噪声。
AR 过程的线性预测方法为：先求得观测数据的自相关函数，然后利用Yule -Walker 方程递推求得模型参数，再根据公式求得功率谱的估计。
Yule -Walker 方程可写成矩阵形式：1.2 LD 算法介绍Levinson-Durbin算法可求解上述问题，其一般步骤为：1) 计算观测值各自相关系数pjjrxx,,1,0),(K=；)0(0xxr=ρ；i=1；2) 利用以下递推公式运算：3) i=i+1，若i>p，则算法结束；否则，返回(2)。
1.3 matlab编程实现以AR模型：x(n)=12x(n−1)−12x(n−2)+w(n)为例，Matlab 程序代码如下：clear; clc;var = 1;noise = var*randn(1,10000);p = 2;coefficient = [1 -0.5 0.5];x = filter(1,coefficient,noise);divide = linspace(-pi,pi,200);for ii = 1:200w = divide(ii);S1(ii) = var/(abs(1+coefficient(2:3)*exp(-j*w*(1:2))'))^2; end[a_p var_p]=Levinson_Durbin(x,p);for ii = 1:200w = divide(ii);Sxx(ii) = var_p/(abs(1+a_p(2:p+1)*exp(-j*w*(1:p))'))^2; endfigure;subplot(2,2,1);plot(divide,S1,'b');grid onxlabel('w');ylabel('功率');title('AR 功率谱');subplot(2,2,2);plot(divide,Sxx,'r-');grid onxlabel('w');ylabel('功率');title('L-D算法估计'); subplot(2,2,3);plot(divide,S1,'b');hold onplot(divide,Sxx,'r--');hold offgrid onxlabel('w');ylabel('功率');title('AR功率谱和算法比较');子函数：Levinson_Durbin.mfunction [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;enda(1)=1;a(2)=-Rxx(2)/Rxx(1);for k=1:p-1 % Levinson-Durbin algorithmvar(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 coeffecientsvar_p = var(p+1); % prediction error power1.4 仿真结果1）p=2时，仿真结果图如下预测系数：[a2(0),a2(1),a2(2)]=[1,−0.5068,0.5031]误差功率：var_p=1.01942）p=20时，仿真结果图如下预测系数：[a2(0),a2(1),a2(2),a2(3),a2(4),……]=[1,−0.5098,0.4999,−0.0066,0.0060,−0.0179,0.0193,……]误差功率：var_p=0.99983）p=50时，仿真结果图如下预测系数：[a2(0),a2(1),a2(2),a2(3),a2(4),……]=[1,−0.4951,0.5178,−0.0145,0.0117,−0.0169,0.0141,……]误差功率：var_p=0.99551.5 结果分析由不同阶数（P值）得到的仿真结果可得：当P的阶数较低时，L-D算法估计AR模型对功率谱估计的分辨率较低，有平滑的效果，从P=2的仿真结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；随着阶数增大，采用L-D算法进行估计后，得到的功率谱会产生振荡，从仿真可以看到，当阶数P较高为50时，估计得到的功率谱与原始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。
现代信号处理大作业
现代信号处理大型作业一．试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：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 模型的阶数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
研究生“现代信号处理”课程大型作业（以下四个题目任选三题做）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），并画出学习曲线。
其X T中，非线性函数采用S型Logistic函数。
1、原理：反向传播(BP)算法：(1)、多层感知器的中间隐层不直接与外界连接，其误差无法估计。
(2)、反向传播算法：从后向前（反向）逐层“传播”输出层的误差，以间接算出隐层误差。
分两个阶段：正向过程：从输入层经隐层逐层正向计算各单元的输出反向过程：由输出层误差逐层反向计算隐层各单元的误差，并用此误差修正前层的权值。
2、流程图：j3、程序：%使用了3层结构，第二层隐藏层4个单元。
2，3层都使用Logisitic函数。
%训练xor数据。
function mlp()f= fopen('XOR.txt');A = fscanf(f, '%g',[3 inf]);A = A;p = A(1:2, :)';%训练输入数据t = A(3, :)';%desire out[train_num , input_scale]= size(p) ;%规模fclose(f);accumulate_error=zeros(1,3001);alpha = 0.5;%学习率threshold = 0.005;% 收敛条件∑e^2 < thresholdwd1=0; wd2=0;bd1=0; bd2=0;circle_time =0;hidden_unitnum = 4; %隐藏层的单元数w1 = rand(hidden_unitnum,2);%4个神经元，每个神经元接受2个输入w2 = rand(1,hidden_unitnum);%一个神经元，每个神经元接受4个输入b1 = rand(hidden_unitnum,1);b2 = rand(1,1);while 1temp=0;circle_time = circle_time +1;for i=1:train_num%前向传播a0 = double ( p(i,:)' );%第i行数据n1 = w1*a0+b1;a1 = Logistic(n1);%第一个的输出n2 = w2*a1+b2;a2 = Logistic(n2);%第二个的输出a = a2;%后向传播敏感性e = t(i,:)-a;accumulate_error(circle_time) = temp + abs(e)^2;temp=accumulate_error(circle_time);s2 = F(a2)*e; %输出层delta值s1 = F(a1)*w2'*s2;%隐层delta值%修改权值wd1 = alpha .* s1*a0';wd2 = alpha .* s2*a1';w1 = w1 + wd1;w2 = w2 + wd2;bd1 = alpha .* s1;bd2 = alpha .* s2;b1 = b1 + bd1;b2 = b2 + bd2;end;%end of forif accumulate_error(circle_time) <= threshold| circle_time>3001 %then break;end;%end of ifend;%end of whileplot(accumulate_error,'m');grid;xlabel('学习次数')ylabel('误差')disp(['计算误差= ',num2str(accumulate_error(circle_time))] ) ;disp(['迭代次数= ',num2str(circle_time)]);%测试a0 = double ([0 0]');n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(['0 0 = ',num2str(a)]);a0 = double ([0 1]');n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(['0 1 = ',num2str(a)]);a0 = double ([1 0]');n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(['1 0 = ',num2str(a)]);a0 = double ([1 1]');n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(['1 1 = ',num2str(a)]);m=0;%---------------------------------------------------------- function [a]= Logistic(n)a = 1./(1+exp(-n));%---------------------------------------------------------- function [result]= F(a)[r,c] = size(a);result = zeros(r,r);for i =1:rresult(i,i) = (1-a(i))*a(i);end;4、实验结果：计算误差= 0.0049993迭代次数= 27060 0 = 0.0231820 1 = 0.9631101 0 = 0.9653901 1 = 0.0433745、学习曲线图：图1.MLP二、试用用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。
1、设计步骤：(1)对Fp 、Fr 进行预畸);();(''Fs Fr tg Fs Fp tg r p ∏=Ω∏=Ω (2)计算'''*r p c ΩΩ=Ω，判断'c Ω是否等于1，即该互补滤波器是否为互补镜像滤波器(3)计算相关系数⎪⎩⎪⎨⎧-==+++=+-=-=ΩΩ=--=偶数)N 为(;21 奇数)N 为 (;;lg /)16/1lg(;150152;1121;1;;])110)(110[(1213090500''02'''211－min 1.0min 1.0i i u q k N q q q q q k k q k k k k r pAr Ap ;)2cos()1(21))12(sin()1(210)1(21'2∑∑∞=∞=+-++-=Ωm m m m m m m i u Nm q u N m 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 i i i =Ω+=β (4)互补镜像滤波器的数字实现;22i i i A αα+-= ;22ii i B ββ+-= 1221,1;1)(N i Z A Z A Z H i i i =++=∏-- 22212,1;1)(N i Z B Z B Z Z H i ii =++=∏--- )];()([21)(21Z H Z H Z H L += 2、程序：function filter2()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=j endab=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endbW(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,'k',w,HP,'m');%hold on;xlabel('数字频率');ylabel('幅度');3、实验结果：图2.两带滤波器4、四带滤波器组程序：function filterfourFp=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=j endb=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)) ;endH2=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,'k',w,LHP,'m',w,HLP,'g',w,HHP,'b');xlabel('数字频率');ylabel('幅度');5、实验结果：图3.四带滤波器组三、根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1）Levinson算法2）Burg算法3）ARMA模型法4）MUSIC算法1、Levinson算法分析：(1)、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：∑--=-≤+=m N n xxN m n m x n x Nm R 10*1),()(1)(ˆ(2)、根据估计所得到的自相关函数，用下面的迭代公式估算AR 模型参数：)()0(*1,2i R a R ki i k k∑=+=σ∑==-+=ki k i k k a i k R a D 00,,0),1(21kkk D σγ=+22121)1(k k k σγσ++-=k i a a a i k k k i k i k ,,2,1,*1,1,,1 =-=-+++γ11,1+++-=k k k a γ(3)、对于AR （p ）模型，按以上述递推公式迭代计算直到p k =+1时为止。