现代数字信号处理仿真作业

实验1:理想采样信号的序列,幅度谱,相位谱,以及改变参数后的图像。

源程序: clc;n=0:50;A=444.128;a=50*sqrt(2.0*pi;T=0.001;w0=50*sqrt(2.0*pi;x=A*exp(-a*n*T.*sin(w0*n*T;close allsubplot(3,2,1;stem(x,’.’;title('理想采样信号序列';k=-25:25;W=(pi/12.5*k;X=x*(exp(-j*pi/12.5.^(n'*k;magX=abs(X;s ubplot(3,2,2;stem(magX,’.’;title('理想采样信号序列的幅度谱';angX=angle(X;subplot(3,2,3;stem(angX;title('理想采样信号序列的相位谱'n=0:50;A=1;a=0.4,w0=2.0734;T=1; x=A*exp(-a*n*T.*sin(w0*n*T;subplot(3,2,4;stem(x,’.’; title('理想采样信号序列'; k=-25:25; W=(pi/12.5*k;X=x*(exp(-j*pi/12.5.^(n'*k; magX=abs(X; subplot(3,2,5; stem(magX,’.’title('理想采样信号序列的幅度谱';0204060-2000200理想采样信号序列020406005001000理想采样信号序列的幅度谱0204060-505理想采样信号序列的相位谱0204060-11理想采样信号序列020406012理想采样信号序列的幅度谱上机实验答案:分析理想采样信号序列的特性产生在不同采样频率时的理想采样信号序列Xa(n,并记录各自的幅频特性,观察频谱‚混淆‛现象是否明显存在,说明原因。

源程序:A=444.128;a=50*pi*sqrt(2.0;W0=50*pi*sqrt(2.0;n=-50:1:50; T1=1/1000;Xa=A*(exp(a*n*T1.*(sin(W0*n*T1;subplot(3,3,1;plot(n,Xa;title('Xa序列';xlabel('n';ylabel('Xa';k=-25:25;X1=Xa*(exp(-j*pi/12.5.^(n'*k;subplot(3,3,2; stem(k,abs(X1,'.';title('Xa的幅度谱';xlabel('k';ylabel('〃幅度';subplot(3,3,3;stem(k,angle(X1,'.';title('Xa的相位谱';xlabel('k';ylabel('相位';T2=1/300;Xb=A*(exp(a*n*T2.*(sin(W0*n*T2;subplot(3,3,4;plot(n,Xb;title('Xb序列';xlabel('n';ylabel('相位';k=-25:25;X2=Xb*(exp(-j*pi/12.5.^(n'*k;subplot(3,3,5; stem(k,abs(X2,'.'; title('Xb 的幅度谱';xlabel('k';ylabel('〃幅度';subplot(3,3,6;stem(k,angle(X2,'.'; title(' Xb 的相位谱';xlabel('k';ylabel('相位';T3=1/200;Xc=A*(exp(a*n*T3.*(sin(W0*n*T3; subplot(3,3,7;plot(n,Xc;title('Xc 序列'; xlabel('n';ylabel('Xc';k=-25:25;X3=Xc*(exp(-j*pi/12.5.^(n'*k;subplot(3,3,8; stem(k,abs(X3,'.'; title('Xc 的幅度谱'; xlabel('k';ylabel('幅度';subplot(3,3,9;stem(k,angle(X3,'.'; title('Xc 的相位谱';xlabel('k';ylabel('相位';-50050-5057X a 序列n X a-500500128X a 的幅度谱k 幅度-50050-55X a 的相位谱k相位-50050-50518X b 序列n 相位-50050051018X b 的幅度谱k 幅度-50050-55X b 的相位谱k相位-50050-505x 1026X c 序列nX c-500500510x 1026X c 的幅度谱k幅度-50050-505X c 的相位谱k相位由图可以看出:当采样频率为1000Hz时,采样序列在折叠频率附近处,无明显混叠。

研究生“现代信号处理”课程大型作业（以下四个题目任选三题做）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），并画出学习曲线。

姓名：潘晓丹学号：班级：A1403492作业1LD 算法实现AR 过程估计1.1 AR 模型p 阶AR 模型的差分方程为：)()()(1n w i nx a n x pi i ，其中)(n w 是均值为0的白噪声。

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

Yule-Walker 方程可写成矩阵形式：1.2 LD 算法介绍Levinson-Durbin算法可求解上述问题，其一般步骤为：1) 计算观测值各自相关系数pjjrxx,,1,0),(；)0(0xxr；i=1；2) 利用以下递推公式运算：3) i=i+1，若i>p，则算法结束；否则，返回(2)。

1.3 matlab编程实现以AR模型：为例，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时，仿真结果图如下预测系数：误差功率：var_p=1.01942）p=20时，仿真结果图如下预测系数：误差功率：var_p=0.99983）p=50时，仿真结果图如下预测系数：误差功率：var_p=0.99551.5 结果分析由不同阶数（P值）得到的仿真结果可得：当P的阶数较低时，L-D算法估计AR模型对功率谱估计的分辨率较低，有平滑的效果，从P=2的仿真结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；随着阶数增大，采用L-D算法进行估计后，得到的功率谱会产生振荡，从仿真可以看到，当阶数P较高为50时，估计得到的功率谱与原始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。

1-2习题1-2图所示为一个理想采样—恢复系统，采样频率Ωs =8π,采样后经过理想低通G jΩ 还原。

解：(1)根据余弦函数傅里叶变换知：)]2()2([)]2[cos(πδπδππ-Ω++Ω=t F ，)]6()6([)]6[cos(πδπδππ-Ω++Ω=t F 。

又根据抽样后频谱公式：∑∞-∞=∧Ω-Ω=Ωk s a a jk j X T j X )(1)(，得到14T= ∑∞-∞=∧--Ω+-+Ω=Ωk a k k j X )]82()82([4)(1ππδππδπ∑∞-∞=∧--Ω+-+Ω=Ωk a k k j X )]86()86([4)(2ππδππδπ所以，)(1t x a ∧频谱如下所示)(2t x a ∧频谱如下所示（2）)(1t y a 是由)(1t x a ∧经过理想低通滤波器)(Ωj G 得到，)]2()2([)()()]([11πδπδπ-Ω++Ω=ΩΩ=∧j G j X t y F a a ，故)2cos()(1t t y a π=(4π) (4π) (4π)(4π)(4π) (4π) Ω-6π-10π-2π 2π0 6π10π)(1Ω∧j X a Ω10π-10π -6π-2π 0 2π6π-14π 14π(4π)(4π) (4π)(4π) (4π) (4π)(4π) (4π))(2Ω∧j X a同理，)]2()2([)()()]([22πδπδπ-Ω++Ω=ΩΩ=∧j G j X t y F a a 故)2cos()(2t t y a π=（3）由题（2）可知，无失真，有失真。

原因是根据采样定理，采样频率满足信号)(1t x a 的采样率，而不满足)(2t x a 的，发生了频谱混叠。

1-3判断下列序列是否为周期序列，对周期序列确定其周期。

（1）()5cos 86x n A ππ⎛⎫=+ ⎪⎝⎭（2）()8n j x n eπ⎛⎫- ⎪⎝⎭=（3）()3sin 43x n A ππ⎛⎫=+ ⎪⎝⎭解：（1）85πω=，5162=ωπ为有理数，是周期序列，.16=N （2）πωπω162,81==，为无理数，是非周期序列； （3）382,43==ωππω，为有理数，是周期序列，8=N 。

实验6 数字滤波器的网络结构一、实验目的：1、加深对数字滤波器分类与结构的了解。

2、明确数字滤波器的基本结构及其相互间的转换方法。

3、掌握用MA TLAB 语言进行数字滤波器结构间相互转换的子函数及程序编写方法。

二、实验原理：1、数字滤波器的分类离散LSI 系统对信号的响应过程实际上就是对信号进行滤波的过程。

因此，离散LSI 系统又称为数字滤波器。

数字滤波器从滤波功能上可以分为低通、高通、带通、带阻以及全通滤波器；根据单位脉冲响应的特性，又可以分为有限长单位脉冲响应滤波器（FIR ）和无限长单位脉冲响应滤波器（IIR ）。

一个离散LSI 系统可以用系统函数来表示：M-m-1-2-m mm=0012m N -1-2-k-k12k k k=1bz b +b z +b z ++b z Y(z)b(z)H(z)====X(z)a(z)1+a z +a z ++a z1+a z ∑∑L L 也可以用差分方程来表示：N Mk m k=1m=0y(n)+a y(n-k)=b x(n-m)∑∑以上两个公式中，当a k 至少有一个不为0时，则在有限Z 平面上存在极点，表达的是以一个IIR 数字滤波器；当a k 全都为0时，系统不存在极点，表达的是一个FIR 数字滤波器。

FIR 数字滤波器可以看成是IIR 数字滤波器的a k 全都为0时的一个特例。

IIR 数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、直接Ⅲ型、级联型和并联型。

FIR 数字滤波器的基本结构分为横截型（又称直接型或卷积型）、级联型、线性相位型及频率采样型等。

本实验对线性相位型及频率采样型不做讨论，见实验10、12。

另外，滤波器的一种新型结构——格型结构也逐步投入应用，有全零点FIR 系统格型结构、全极点IIR 系统格型结构以及全零极点IIR 系统格型结构。

2、IIR 数字滤波器的基本结构与实现 （1）直接型与级联型、并联型的转换 例6-1 已知一个系统的传递函数为-1-2-3-1-2-38-4z +11z -2z H(z)=1-1.25z +0.75z -0.125z将其从直接型（其信号流图如图6-1所示）转换为级联型和并联型。

数字信号处理大作业班级学号姓名陈志豪一．要求本次作业要求对一段音乐进行处理，该音乐包含了蜂鸣噪声，根据该段音乐，我们需处理以下问题：1. 利用matlab软件对audio1211.wav音频信号进行数字信号采样，分别对采样后的信号进行时/频域分析，并提供仿真图和分析说明；2. 设计合理的数字滤波器，滤去音频信号中的蜂鸣音，给出详细设计流程，并提供频域仿真图和分析说明；3. 将数字滤波后的数字信号转换成wav格式音频文件二．分析(1）通过播放所给音乐文件，很明显能听出wav文件中包含蜂鸣噪音，所以我们应该先分析频谱。

在matlab下可以用函数wavread/audioread读入语音信号进行采样，我们可以通过wavread得到声音数据变量x和采样频率fs、采样精度nbits，在读取声音信号之后，利用读出的采样频率作为参数，这段音频读出的采样精度为16，fs为44100hz，所以我们将此后采集时间、fft的参数设置为fs，也就是44100hz。

最后我们通过plot函数绘制出了音频信号与时间的关系图pic1，使用fft函数进行fft处理。

处理后的信号频谱pic2，如下所示图1.音频信号与时间的关系图从图1横坐标我们看到t在9-10s之间截止，与我们在音乐播放器中显示的时间一致。

图2.fft之后得到的频域分析结果图3.噪声读取图2为运用fft后得到的处理结果，可以从中读取到，在293.7hz、4671hz附近幅值突然增大，可以确定为噪声干扰。

所以我们应该针对频率附近进行滤波。

如果针对性进行滤波处理，应该使用低通滤波器进行处理，去除这部分的噪音。

之后需要选定滤波器并进行程序设计，在4671hz附近进行滤波,去除蜂鸣杂音。

（2）我们需要对蜂鸣音进行除去，自然需要用到滤波器。

所以第二步我们需要设计滤波器并给出详细流程。

在第一问的频谱分析中，通过FFT我们已经知道噪音所在，所以我们需要针对这个问题设计参数。

在这里我们选用巴特沃斯低通滤波器进行处理，我们需要设定好的参数有通带边界频率、阻带边界频率、通带最大衰减和通过阻带的最小分贝数（由buttord 在matlab定义得）。

7.1 给定一理想低通FIR 滤波器的频率特性⎪⎩⎪⎨⎧≤<≤=πππw w e H jwd 4/04/1)(现希望用窗函数（矩形窗和汉明窗）法设计该滤波器，要求具有线性相位。

假定滤波器系数的长度为29点，即M/2=14。

由公式：)2()2sin(21)(2M n w M n dw een h c w w jwnMw jd cv--==⎰-+-ππ所以：)14(]25.0*)14sin[()(--=n n n h d ππ在MA TLAB 中输入指令：hd=sin((n-14).*0.25*pi)./((n-14)*pi) >> n=[0:14];>> hd=sin((n-14).*0.25*pi)./((n-14)*pi) hd =Columns 1 through 10-0.0227 -0.0173 0.0000 0.0205 0.0318 0.0250 -0.0000 -0.0322 -0.0531 -0.0450Columns 11 through 150.0000 0.0750 0.1592 0.2251 NaN 即：)(n h d25.0)14(0.2251 )15()13(0.1592)16()12(0.0750)17()11(0000.0)18()10( 0.0450)19()9( -0.0531)20()8( -0.0322)21()7( -0.0000)22()6(0.0250)23()5( 0.0318)24()4(0.0205)25()3( 0.0000)26()2( -0.0173)27()1( -0.0227)28()0(=============================d d d d d d d d d d d d d d d d d d d d d d d d d d d d d h h h h h h h h h h h h h h h h h h h h h h h h h h h h h1、当选用窗函数为矩形窗时，h(n)= hd(n)所以滤波器的系数为[-0.0227 -0.0173 0.0000 0.0205 0.03180.0250 -0.0000 -0.0322 -0.0531 -0.0450 0.0000 0.0750 0.15920.2251 0.2500 0.2251 0.1592 0.0750 0.0000 -0.0450 -0.0531-0.0322 -0.0000 0.0250 0.0318 0.0205 0.0000 -0.0173 -0.0227 ]MATLAB指令：N=29Wn=0.5*pi/2;window=boxcar(N+1);b=fir1(N,Wn/pi,'low',window);freqz(b,1,512)可得幅频和相频响应曲线：2、当选用窗函数为汉宁窗时（书305页公式），matlab 的算法：wn=0.5-(0.5).*cos((2*pi*n)/29)； hn=hd.*wn>> window=hanning(N+1);3、当选用窗函数为汉明窗时（书306页公式）， h(n)= hd(n)*w(n)： 其中,1,........3,2,1,0)2cos(*46.054.0)(-=-=N n Nnn w πmatlab 的算法：wn=0.54+(0.46).*cos((2*pi*n)/29)； hn=hd.*wn>> window=hamming(N+1);7.2给定一理想带阻滤波器的频率特性：现希望用窗函数（矩形窗、汉明窗和布莱克曼窗）法设计该滤波器，要求具有线性相位。

MATLAB在音频信号处理中的应用 1.引言 MATLAB是美国Math Works公司推出的一种面向工程和科学计算的交互式计算软件、它以矩阵运算为基础，把计算、可视化、程序设计融合在一个简单易用的交互式工作环境中，是一款数据分析和处理功能都非常强大的工程实用软件。MATLAB软件除可以对文本数据进行操作外，还可以对图像文件、WAV类型的音频文件和AVI类型的视频文件进行读写.借助windows系统自带的音频录制播放设备和MATLAB软件，可以完成音频信号的综合分析处理。 音频是信号的一种，处理数字音频信号也是一种数字信号分析与处理。人能听到的声音频率范围为20Hz到20kHz，语音信号频率为300Hz到3.4kH。木文首先对音频信号的采集与频谱分析作介绍，然后介绍了用MATLAB处理音频信号的整体流程，并且对我们平常所听的音乐文件进行一系列的处理，然后用我们的耳朵亲自去听，去感受不同处理方法的效果。

2.音频信号的采集与频谱分析 进行音频的频谱分析时，首先要对声音信号进行采。MATLAB 的数据采集工具箱提供了一整套命令和函数,通过调用这些函数和命令,可直接控制声卡进行数据采集。Windows 自带的录音机程序也可驱动声卡来采集语音信号，并能保存为WAV 格式文件供MATLAB 相关函数直接读取、写入或播放。函数[y, fs, bits] =wavread(‘Blip’, N)。用于读取音频，对音频信号进行采集，采样值放在参数y中，fs表示每秒采样点数，即采样频率，bits表示每个采样点在编码时所占位数。N表示采样点总数。参数’Blip’为音频所在地址，如:'C:\yinpinl'。调用函数fft可对己采集音频信号进行时频转换，通过函数abs()和angle()可分别得到信号频谱的幅频图和相频图。对放在C盘目录下

手机短信铃声”test.wav”进行采样与频谱分析。 其体代码如下: clear all; [y0,Fs0,nbits0]=wavread('test.wav'); a=round(length(y0)); y01=fft(y0);y011=abs(y01); t0=linspace(0,a/Fs0,a); figure(1);subplot(2,1,1);plot(t0,y0); title('(a)','fonts',10.5 ,'position' ,[25,-1.5]); xlabel('时间/s','fonts',10.5 ,'position',[45,-1.4]); ylabel('幅值n','fonts', 10.5,'position',[-3,0.7]);