matlab与数字信号处理课程设计报告.docx

1、matlab与数字信号处理课程设计报告数字信号处理专业课程设计报告书 实验报告题目四 :Using the bilinear transformation and a lowpass analog Butterworth prototype, design a highpass digital filter operating at a rate of 20kHz and having passband extending to 5kHz with maximum passband attenuation of 0.5dB,and stopband ending at 4kHz with a

2、minimum stopband attenuation of 10dB.备注：题目3,4要求：实验报告中要求写出对应的滤波器H(z)，并在H(z)表达式中将共轭极点对组成二阶基本节，以极点在Z平面上分布顺序写出H(z)形式并将各二阶基本节系数以顺序列表。画出幅度频谱图的|H()|及其以(dB)为单位的幅度谱图。二：实验目的1）熟练掌握低通滤波器的设计方法。2）学会利用低通滤波器设计高通滤波器。3）掌握用双线形变换法设计数字高通滤波器的方法。4）熟悉MATLAB提供的各种滤波器设计函数。5）掌握各种关于滤波器的幅度频谱设计函数。三：实验原理本题利用双线性变换法和巴特沃斯低通滤波器来设计数字高通

3、滤波器：双线形变换法是利用s=2*（1-z-1）/T*（1+z-1）将s域转换到z域，从而得到系统函数H（Z）。根据所要设计要求将高通数字滤波器指标转化为低通模拟滤波器技术指标，主要利用双线性变换式=2/ Ttan(W/2)。滤波器设计中主要用到的函数：Buttord函数用来选择巴特沃斯滤波器最小阶数，调用方式如下：n,wn=buttord(Wp,Ws,rp,rs,s) ：返回符合要求性能指标的数字滤波器最小阶数n和巴特沃斯滤波器截止频率wn； n,wn=buttord(Wp,Ws,rp,rs)：（同上）此处Wp,Ws都是归一化频率。z,p,k=buttap(n)：返回设计的巴特沃斯滤波器的零

4、点（z），极点（p）和增益（k），n为滤波器阶数。b,a=zp2tf(z,p,k)：零极点增益滤波器参数转换为传输函数形式，b，a分别为传输函数的分子分母。bt,at=lp2hp(b,a,wn)：模拟低通滤波器参数转化为模拟高通滤波器参数。bd,ad=bilinear(bt,at,Fs); %模拟高通滤波器参数转化为数字高通滤波器参数。Fs为采样频率。sos,g=tf2sos(bd,ad,order)：将传递函数模型转化为二次分式模型sos。order指定sos中行的顺序： Up：首行中所包含的极点离原点最近，离单位圆最远； Down：首行中所包含的极点离原点最远，离单位圆最近；与画图有关的函

5、数：（1）Zplane(b,a)：绘制系统零极点图；（2）求解数值滤波器频率响应函数：1）freqz(b,a,n)：无输出参数，直接在当前命令窗口绘制频率响应的幅频响应（dB形式的）和相频响应。2）H,W=freqz(b,a,n)：返回数字滤波器的n点复频率响应。四：实验步骤简述 （1）先将数字高通滤波器技术指标转化为模拟低通滤波器技术指标。（2）确定滤波器最小阶数.（3）确定零极点增益。（4）由零极点增益确定模拟低通滤波器传输函数。（5）将模拟低通滤波器传输函数转换成模拟高通滤波器传输函数。（6）利用双线性变法将模拟高通滤波器转换成数字滤波器。（7）利用传递函数模型求其二次分式模型sos。（

6、8）绘制相关幅度频谱图。五：程序框图 六：源程序：clc; close; clear;Fs=20000;T=1/Fs;%数字高通技术指标转化为模拟低通技术指标 fs=4000; fp=5000;rp=0.5; rs=10;ws=fp*2*pi*T; wp=fs*2*pi*T;Wp=2/T*tan(wp/2); Ws=2/T*tan(ws/2);%*滤波器设计*N,wn=buttord(Wp,Ws,rp,rs,s)%根据技术指标选择数字滤波器最小阶数，巴特沃斯截止频率z,p,k=buttap(N);%返回设计的巴特沃斯滤波器的零点极点和增益b,a=zp2tf(z,p,k);%零极点增益滤波器参数

7、转换为传输函数形式bt,at=lp2hp(b,a,wn);%模拟低通转化为模拟高通bd,ad=bilinear(bt,at,Fs); %模拟高通转化为数字高通%模拟低通和数字高通二阶基本节系数顺序表sos,g=tf2sos(b,a,up)sos,g=tf2sos(bd,ad,up)%*相关设计图*%模拟低通和数字高通零极点图figure(1); subplot(121);zplane(b,a); gtext(模拟低通零极点图);set(gcf,defaultAxesLinestyleOrder,-|:);set(gcf,defaultAxesColorOrder,1 0 0;0 0 1);%按

8、先后顺序对图形线条进行颜色设置subplot(122);zplane(bd,ad);gtext(数字高通零极点图);pause; figure(2);%滤波器幅频特性图freqz(bd,ad,100);axis(0 1 -40 0);pause; figure(3);H,W=freqz(bd,ad,100);%求出数字滤波器幅频特性plot(W*Fs/2/pi,abs(H);axis(0 8000 0 1.5);xlabel( 频率 f/HZ ); ylabel(归一化幅度 );title(滤波器); grid on;pause;close all;七：程序结果及图表滤波器阶数N=7数字高通滤

9、波器二阶基本节系数按“up（首行中所包含的极点离原点最近，离单位圆最远）”排列:sos = 1.0000 -0.9900 0 1.0000 -0.0783 0 1.0000 -2.0181 1.0182 1.0000 -0.1647 0.0582 1.0000 -2.0044 1.0045 1.0000 -0.1927 0.23771.0000 -1.9874 0.9876 1.0000 -0.2552 0.6396g =0.0279分子系数1.0000 -0.9900 0 1.0000 -2.0181 1.0182 1.0000 -2.0044 1.00451.0000 -1.9874 0.

10、9876分母系数:1.0000 -0.0783 0 1.0000 -0.1647 0.0582 1.0000 -0.1927 0.2377 1.0000 -0.2552 0.6396图（1）：图（2）：以dB为单位的幅度谱图及相位频率图： 图（3）：数字高通滤波器幅度频谱图：八：实验总结 在本次试验中，我熟练的掌握了用模拟低通巴特沃斯滤波器设计数字高通滤波器的方法。本次设计重在设计思维的建立及其所用函数的熟练掌握。在进行具体设计时，主要抓住“给定条件”，“目的要求”展开设计构架。根据“目的要求”：用模拟低通巴特沃斯滤波器设计数字高通滤波器；可知主要设计出传输函数H(Z)。再结合“给定条件”：数

11、字高通滤波器技术指标；设计出系统框图。最终应用相关函数编辑出具体代码，在进行调试。在本次试验中我遇到一个小问题：二阶基本节系数矩阵sos为什么是“4行3列”，而不是“1行3列”，最终我找到了答案：在本题中二次分式形式的传输函数H(Z)是由4个二阶基本节串联而成，因此是“4行3列”。本次实验我掌握了利用双线性变换法设计滤波器的方式方法，相对成功的设计出要求的高通数字滤波器，达到了实验要求及实验目的。题目七：该题目的目的是说明一个PN扩频信号在抑制正弦干扰中的有效性。现考虑下图所示的二进制通信系统，对信号发生器的输出乘上一个二进制（1）PN序列。同一个二进制PN序列用来与解调器输入相乘，因此消除了

12、这个PN序列在期望信号上的影响。信道将传送信号受到一宽带加性噪声序列(n)和一正弦干扰序列i（n）=Asin0n, 000.5; %使每个元素的值大于0.5%-PN序列的产生function p_s=pn(n,M) a=zeros(1,n); a(n)=1; for i=1:1000 y(i)=a(n); temp=mod(a(n-1)+a(n),2); for j=n:-1:2 a(j)=a(j-1); %改变a(n)序列 end a(1)=temp; %改变a(n)序列 end p_s=y(1:M); %取a(n)序列的最后一位元素生成pn序列%-正弦信号的产生function in=si

13、nout(A,w0,M)j=1:M;in=A*sin(w0*j);%-随机噪声的产生function ao=aoout(M)ao=rand(1,M); %-计算序列点数时信噪比与差错率的关系 function =p(n,A,M,Z,B,w0,pg)er=randout(M);r=pn(n,M);for k=1:25 %循环增加Z（信号的幅度） Z=Z+B; %使调制后的序列幅度Z大于噪声幅度B g=xor(er,r)*Z; %调制（逻辑异或）差分调制 q=sinout(A,w0,M); x=g+q; %调制后序列加正弦信号 re=B*aoout(M); %生成噪声信号 x1=re+x; %加入

14、噪声 r22=x1pg; %设置一个门限pg，滤除噪声 r2=x1Z/3; %设置一个门限Z/3，滤除噪声 r3=xor(r2,r); %解调 r33=xor(r22,r); i=1:M; r4=er(i)=r3(i); %检测 r44=er(i)=r33(i); r5=length(find(r4); %计算错误的个数（计算r4中元素取1的个数） r55=length(find(r44); cc(k)=r5/M; %错误的概率 ccc(k)=r55/M; %*计算信号的平均功率 i=1:M; xg=0; xg=xg+g(i).2; xgp=xg/M; %*计算噪声的平均功率 j=1:M; z

15、a=0; za=za+q(j).2+re(j).2; zap=za/M; %*计算信噪比 sn(k)=xgp/zap; %-不加PN序列 h=er*Z; %原始信号乘幅度因子 x2=h+re+q; %信号序列加正弦信号和噪声信号（不与PN序列调制） r66=x2pg; %设置一个门限pg，滤除噪声 r6=x2Z/3; %设置一个门限z/3，滤除噪声 %误码的概率 i=1:M; r77=er(i)=r66(i); r7=er(i)=r6(i); r88=length(find(r77); r8=length(find(r7); cd(k)=r8/M; cdd(k)=r88/M; %计算信号的平均

16、功率 i=1:M; xgo=0; xgo=xgo+h(i).2; xgog=xgo/M; %计算噪声的平均功率 j=1:M; os=0; os=os+q(j).2+re(j).2; osst=os/M; sd(k)=xgog/osst; %计算信噪比end%*画信噪比与差错率曲线figure(1); h1=figure(1); set(h1,name,判决门限为均值pg时的差错率,color,0 1 0.5);subplot(211);stem(sn,ccc,vk); hold on; stem(sd,cdd,.m); %输出有PN序列的图和输出没有PN序列的图xlabel(信噪比); yla

17、bel(差错率);title(信噪比与差错率曲线); grid on;subplot(212);plot(sn,ccc,:k); hold on; plot(sd,cdd,:m);legend(有PN序列的差错率曲线,无PN序列的差错率曲线);figure(2);h2=figure(2); set(h2,name,判决门限为Z/3时的差错率,color,0 1 0.5);subplot(211);stem(sn,cc,or);hold on; stem(sd,cd,xb); %输出有PN序列的图和输出没有PN序列的图xlabel(信噪比); ylabel(差错率);title(信噪比与差错率曲

18、线); grid on;subplot(212);plot(sn,cc,:r); hold on; plot(sd,cd,:b);legend(有PN序列的差错率曲线,无PN序列的差错率曲线);七：程序结果及图表菜单：正弦干扰信号角频率在0pi间且w0=pi*wl，请输入w1=请选择解调器每比特的样本数M：1-M=502-M=1003-M=5004-M=10000-退出此题运行：当w1=0.8解调器每比特的样本数M：4-M=1000：（1）判决门限为Z/3(Z：进入解调器的传输信号幅度)时：（2）判决门限为pg(M取不同值时进入解调器的传输信号幅度的平均值)时：当w1=0.2解调器每比特的样本

19、数M：4-M=1000:(1)(2):请选择解调器每比特的样本数M：1-M=50请选择解调器每比特的样本数M：2-M=100请选择解调器每比特的样本数M：3-M=500结论：1： 在通信系统中，PN扩频信号能有效抑制正弦干扰信号。 2：不同的噪声门限对差错率有一定影响。噪声门限为pg所得的曲线起伏小，但有PN序列的图和没有PN序列的图对比明显；z/3所得的曲线起伏大，但有PN序列的图和没有PN序列的图贴的较近，对比不明显。 3：正弦干扰信号角频率越大，同一通信系统的差错率越高。八：实验总结在做这道题目时，我相对苦恼。因为这是关于一个通信系统仿真实验，我有通信原理课程的学习基础，但关于“PN序列

20、”，我是一无所知，另外又找不到相关方面的有利资料，一时难以下手。最后还是通过对现成的PN序列代码的学习理解才得以继续此题。本道题有题目提供的系统构架，整体思路及框图设计相对容易。重在差错率求解程序设计上，在此段程序设计时重在噪声门限的选取，不同的门限对差错率有很大影响，同时也决定着绘制出信噪比差错率曲线图的对比鲜明度。在做本题时我设置了2种门限：pg是取各次噪声幅度占总信号幅度的均值； z/3是根据观察各次噪声幅度大约占总信号幅度的1/3所取。最终pg所得的曲线起伏小，但有PN序列的图和没有PN序列的图对比明显；z/3所得的曲线起伏大，但有PN序列的图和没有PN序列的图贴的较近，对比不明显。通

21、过本次试验，我对PN扩频信号在通信系统中的用法及作用有了一定的认识同时也加深了我对通信系统的认识。题目八：设原N点序列的DFT为验证以下序列的DFT与X(k)的关系：1：, 2：x2(n)=x(2n)，3：x3(n)=x(2n+1)，4：x4(n)=x(n)NR2N(n)，5：x5(n)=x(n)+x(n+N/2)，6：x6(n)=x(n)-x(n+n/2) 7： 二：实验目的通过输出原序列x（n）的DFT验证题目中由x（n）变形后得到的不同序列的DFT和X（K）的关系。三：实验原理 根据原序列x（n）的DFTX(K)与题目中不同序列的DFT得到不同情况下时域与频域的变化关系。按照题目要求确定

22、新序列及其长度，利用快速傅里叶变换函数fft,得到其离散傅里叶变化。四：实验步骤简述（1）：确定原N点序列x（n），利用fft函数求其DFT；（2）：根据题目要求确定新序列y（n）及其长度 “N“，利用fft函数求其DFT；（3）：根据新序列长度，利用画图函数做出相应序列时域、频域图。五：程序框图六：源程序function =main()clear all; clc;while(1)Q=input(请选择输入序列类型M：n1-函数序列n2-任意值序列n0-退出此题n); %参数设定：if Q=1 N=input(输入序列长度N=); n=1:N; x=input(输入函数序列x(n)=); j

23、h(x,N);elseif Q=2 x=input(输入任意值序列x(x1,x2,x3,x4.):); N=length(x) jh(x,N);elseif Q=0 break;endendfunction =jh(x,N)m=fix(N/2);mm= fix(N-1)/2);%*生成序列x(n)并算出其DFT n=1:N; xk=fft(x(n),N); xxk %*生成序列x1(n)并算出其DFTy1=zeros(1,N); for i=1:Nif (mod(i,2)=0) y1(i)=x(i/2); else y1(i)=0; end endyk1=fft(y1,N); y1yk1%*生成序列x2(n)并算出其DFTy2=zeros(1,m); for i=1:m y2(i)=x(2*i); end n=1:m; yk2=fft(y2,m); y2yk2%*生成序列x3(n)并算出其DFTy3=zeros(1,mm); for i=1:mm y3(i)=x(2*i+1);end n=1:mm ; yk3=fft(y3,mm); y3yk3figure(2); %*生成序列x4(n)并算出其DFTy4=zero