当前位置：文档视界 › matlab与数字信号处理课程设计报告

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 minimum stopband attenuation of 10dB.

备注：题目3,4要求：实验报告中要求写出对应的滤波器H(z)，并在H(z)表达式中将共轭极点对组成二阶基本节，以极点在Z平面上分布顺序写出H(z)形式并将各二阶基本节系数以顺序列表。画出幅度频谱图的|H(ω)|及其以(dB)为单位的幅度谱图。

二：实验目的

1）熟练掌握低通滤波器的设计方法。

2）学会利用低通滤波器设计高通滤波器。

3）掌握用双线形变换法设计数字高通滤波器的方法。

4）熟悉MATLAB提供的各种滤波器设计函数。

5）掌握各种关于滤波器的幅度频谱设计函数。

三：实验原理

本题利用双线性变换法和巴特沃斯低通滤波器来设计数字高通滤波器：

双线形变换法是利用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)：返回设计的巴特沃斯滤波器的零点（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：首行中所包含的极点离原点最远，离单位圆最近；

与画图有关的函数：

（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。

（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);%零极点增益滤波器参数转换为传输函数形式

[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]);%按先后顺序对图形线条进行颜色设置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

数字高通滤波器二阶基本节系数按“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.2377

1.0000 -1.9874 0.9876 1.0000 -0.2552 0.6396

g =

0.0279

分子系数

1.0000 -0.9900 0

1.0000 -

2.0181 1.0182

1.0000 -

2.0044 1.0045

1.0000 -1.9874 0.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 为单位的幅度谱图及相位—频率图：

Normalized Frequency (?π rad/sample)

P h a s e (d e g r e e s

)

Normalized Frequency (?π rad/sample)

M a g n i t u d e (d B )

图（3）：数字高通滤波器幅度频谱图：

00.51

1.5

频率 f/HZ

归一化幅度

滤波器

八：实验总结

在本次试验中，我熟练的掌握了用模拟低通巴特沃斯滤波器设计数字高通滤波器的方法。本次设计重在设计思维的建立及其所用函数的熟练掌握。

在进行具体设计时，主要抓住“给定条件”，“目的要求”展开设计构架。根据“目的要求”：用模拟低通巴特沃斯滤波器设计数字高通滤波器；可知主要设计出传输函数H(Z)。再结合“给定条件”：数字高通滤波器技术指标；设计出系统框图。最终应用相关函数编辑出具体代码，在进行调试。

在本次试验中我遇到一个小问题：二阶基本节系数矩阵sos为什么是“4行3列”，而不是“1行3列”，最终我找到了答案：在本题中二次分式形式的传输函数H(Z)是由4个二阶基本节串联而成，因此是“4行3列”。

本次实验我掌握了利用双线性变换法设计滤波器的方式方法，相对成功的设计出要求的高通数字滤波器，达到了实验要求及实验目的。

题目七：该题目的目的是说明一个PN扩频信号在抑制正弦干扰中的有效性。现考虑下图所示的二进制通信系统，对信号发生器的输出乘上一个二进制（ 1）PN序列。同一个二进制PN序列用来与解调器输入相乘，因此消除了这个PN序列在期望信号上的影响。信道将传送信号受到一宽带加性噪声序列ω(n)和一正弦干扰序列i（n）=Asinω0n, 0<ω0<π的污损。可以假定A≥M,这里M是来自调制器的每比特的样本数。用和不用PN序列执行这个仿真系统，并在条件A≥M下，对不同的M值（如M＝50，100，500，1000）测量差错率。画出每种情况的差错率曲线，作比较并说明结论。说明这个PN序列对于正弦干扰信号的效果。由此说明为什么PN扩频系统在正弦干扰信号存在下优于常规的二进制通信系统。

二：实验目的

（1）利用matlab实现一个简单的通信系统。

（2）了解扩频序列及其实现方法。

（3）在一个简单的通信系统，信号通过一个有噪声的系统，对加上PN序列和不加PN序列后的差错率比较，说明PN序列能有效的抑制正弦干扰信号，即PN扩频系统在正弦干扰信号存在下优于常规的二进制通信系统。

三：实验原理

原始信号S(n)与PN扩频序列相乘（异或实现）后，叠加正弦信号噪声，再加一般噪声后与PN序列相乘恢复出经过系统后的信号，选择不同的调制解调器的每比特的样本数，使恢复信号与原输入信号相比较后得到差错率。根据差错率可以证明PN扩频信号是否有利于通信系统中正弦噪声的抑制。

四：实验步骤简述

（1）理解PN序列概念，掌握产生二进制PN序列的方法。

（2）编写输入信号及正弦干扰信号和噪声信号子函数。

（3）编写主函数，按照通信系统模型，通过对子函数的调用实现有PN序列的通信系统和无PN序列的通信系统：

1): 有PN扩频序列时：原始信号S(n)与PN扩频序列相乘（异或实现）后，叠加正弦信号噪声，

再加一般噪声完成信号的调制和在信道中的传输；无PN扩频序列时：原始信号S(n)叠加正弦信号噪声，再加一般噪声完成信号的调制和在信道中的传输；

2): 设置合适门限，滤除噪声，通过与PN扩频序列相乘（有PN扩频序列时）和不与PN扩频序列

相乘（无PN扩频序列时）解调出信号；

3): 比较输入前后信号的差错率。

（4）画出输出不同门限下的“信噪比——差错率”图像。

五：程序框图

六：源程序

function []=diqi7()

wl=input('正弦干扰信号角频率在0~pi间且‘w0=pi*wl’，请输入w1=');

w0=pi*wl;

n=4;pg0=0;pg1=0;pg2=0;pg3=0;

while(1)

N=input('请选择解调器每比特的样本数M：\n1----M=50\n2----M=100\n3----M=500\n4----M=1000\n0----退出此题\n');

%参数设定：

if N==1

M=50; %原始信号的长度

A=60; %正弦信号的幅度

Z=50; %原始信号的幅度因子

B=15; %噪声幅度

pg0=pg0+15/(50+15);

p(n,A,M,Z,B,w0,pg0);

elseif N==2

M=100; A=110; Z=100; B=30; pg1=pg1+30/(100+30); p(n,A,M,Z,B,w0,pg1);

elseif N==3

M=500; A=550; Z=500; B=80; pg2=pg2+80/(500+80); p(n,A,M,Z,B,w0,pg2);

elseif N==4

M=1000; A=1100; Z=1000; B=130; pg3=pg3+130/(1000+130); p(n,A,M,Z,B,w0,pg3);

elseif N==0

clc; break;

end

pg=(pg0+pg1+pg2+pg3)/4 ;

%----------------------二进制数据的产生

function [er]=randout(M)

n=0:M;

er=rand(1,M);

j=1:M;

er(j)=er(j)>0.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]=sinout(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; %加入噪声

r22=x1>pg; %设置一个门限pg，滤除噪声

r2=x1>Z/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; za=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=x2>pg; %设置一个门限pg，滤除噪声

r6=x2>Z/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;

%计算信号的平均功率

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('信噪比'); ylabel('差错率');

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('信噪比与差错率曲线'); grid on;

subplot(212);

plot(sn,cc,':r'); hold on; plot(sd,cd,':b');

legend('有PN序列的差错率曲线','无PN序列的差错率曲线');

七：程序结果及图表

菜单：

正弦干扰信号角频率在0~pi间且‘w0=pi*wl’，请输入w1=

请选择解调器每比特的样本数M：

1----M=50

2----M=100

3----M=500

4----M=1000

0----退出此题

运行：

当w1=0.8

解调器每比特的样本数M：4----M=1000：

（1）判决门限为Z/3(Z：进入解调器的传输信号幅度)时：

0.10.20.3

0.4信噪比

判决门限为Z /3时的差错率

信噪比与差错率曲线

024681012

0.10.20.30.4

（2

）判决门限为pg(M 取不同值时进入解调器的传输信号幅度的平均值)时：

024

681012

0.10.20.3

0.4信噪比

判决门限为均值p g 时的差错率

信噪比与差错率曲线

0.29

0.30.310.320.330.34

当w1=0.2

解调器每比特的样本数M ：4----M=1000: (1)

0.10.20.3

0.4信噪比

判决门限为Z /3时的差错率

信噪比与差错率曲线

012345678910

0.10.20.30.4

(2):

01234

5678910

0.10.20.3

0.4信噪比

判决门限为均值p g 时的差错率

信噪比与差错率曲线

012345678910

0.29

0.30.310.32

请选择解调器每比特的样本数M ：1----M=50

0.10.20.3

0.4信噪比

判决门限为Z /3时的差错率

信噪比与差错率曲线

0510********

0.10.20.30.4

0510

15202530

0.10.20.3

0.4判决门限为均值p g 时的差错率

信噪比与差错率曲线

1015202530

0.30.310.320.330.340.35

请选择解调器每比特的样本数M ：2----M=100

0.10.20.3

0.4信噪比

判决门限为Z /3时的差错率

信噪比与差错率曲线

051015202530354045

0.10.20.30.4

051015

202530354045

0.10.20.3

0.4信噪比

判决门限为均值p g 时的差错率

信噪比与差错率曲线

0.29

0.30.310.320.330.34

请选择解调器每比特的样本数M ：3----M=500

信噪比

判决门限为Z /3时的差错率

信噪比与差错率曲线

0.10.20.3

0.4信噪比

判决门限为均值p g 时的差错率

信噪比与差错率曲线

0.27

0.280.290.30.31

结论：1：在通信系统中，PN 扩频信号能有效抑制正弦干扰信号。

2：不同的噪声门限对差错率有一定影响。噪声门限为pg 所得的曲线起伏小，但有PN 序列的图和没有

PN 序列的图对比明显；z/3所得的曲线起伏大，但有PN 序列的图和没有PN 序列的图贴的较近，对比

不明显。

3：正弦干扰信号角频率越大，同一通信系统的差错率越高。八：实验总结

在做这道题目时，我相对苦恼。因为这是关于一个通信系统仿真实验，我有通信原理课程的学习基础，但关于“PN 序列”，我是一无所知，另外又找不到相关方面的有利资料，一时难以下手。最后还是通过对现成的PN 序列代码的学习理解才得以继续此题。

本道题有题目提供的系统构架，整体思路及框图设计相对容易。重在差错率求解程序设计上，在此段程序设计时重在噪声门限的选取，不同的门限对差错率有很大影响，同时也决定着绘制出信噪比差错率曲线图的对比鲜明度。在做本题时我设置了2种门限：pg 是取各次噪声幅度占总信号幅度的均值； z/3是根据观察各次噪声幅度大约占总信号幅度的1/3所取。最终pg 所得的曲线起伏小，但有PN 序列的图和没有PN 序列的图对比明显；z/3所得的曲线起伏大，但有PN 序列的图和没有PN 序列的图贴的较近，对比不明显。

通过本次试验，我对PN 扩频信号在通信系统中的用法及作用有了一定的认识同时也加深了我对通信系统的认识。

题目八：设原N 点序列的DFT 为验证以下序列的DFT 与X(k)的关系： 1：1(2),()0x n n x n n ?=?

?为偶数，为奇数

, 2：x 2(n)=x(2n)，3：x 3(n)=x(2n+1)，

4：x 4(n)=x((n))N R 2N (n)，5：x 5(n)=[x(n)+x(n+N/2)]，6：x 6(n)=[x(n)-x(n+n/2)]n

N W

7：7(),0()0,21

x n n N

x n N n N ≤≤?=?

≤≤-?

二：实验目的

通过输出原序列x （n ）的DFT 验证题目中由x （n ）变形后得到的不同序列的DFT 和X （K ）的关系。三：实验原理

根据原序列x （n ）的DFT ——X(K)与题目中不同序列的DFT 得到不同情况下时域与频域的变化关系。按照题目要求确定新序列及其长度，利用快速傅里叶变换函数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)=');

jh(x,N);

elseif Q==2

x=input('输入任意值序列x([x1,x2,x3,x4...]):');

N=length(x)

jh(x,N);

elseif Q==0

break;

end

function []=jh(x,N)

m=fix(N/2);mm= fix((N-1)/2);

%*********生成序列x(n)并算出其DFT

n=1:N;

xk=fft(x(n),N);

%********生成序列x1(n)并算出其DFT

y1=zeros(1,N);

for i=1:N

if (mod(i,2)==0)

y1(i)=x(i/2);

else

y1(i)=0;

end

yk1=fft(y1,N);

yk1

%********生成序列x2(n)并算出其DFT

y2=zeros(1,m);

for i=1:m

y2(i)=x(2*i);

end

n=1:m;

yk2=fft(y2,m);

yk2

%********生成序列x3(n)并算出其DFT

y3=zeros(1,mm);

for i=1:mm

y3(i)=x(2*i+1);

end

n=1:mm ;

yk3=fft(y3,mm);

yk3

figure(2);

%*********生成序列x4(n)并算出其DFT

y4=zeros(1,2*N); for i=1:N

y4(i)=x(i);

end

for i=N+1:2*N y4(i)=x(i-N);

end

n=1:2*N;

yk4=fft(y4,2*N);

yk4

%*********生成序列x5(n)并算出其DFT

y5=zeros(1,m);

for i=1:m

y5(i)=x(i)+x(i+m);

end

n=1:m;

yk5=fft(y5,m);

yk5

%**********生成序列x6(n)并算出其DFT

y6=zeros(1,m);

for i=1:m

y6(i)=(x(i)-x(i+m))*(exp(-j*2*pi/N*i)); end

yk6=fft(y6,m); n=1:m;

yk6

%*********生成序列x7(n)并算出其DFT画图

y7=zeros(1,2*N); for i=1:N

y7(i)=x(i);

end

for i=N+1:2*N y7(i)=0;

end

n=1:2*N;

yk7=fft(y7,2*N);

yk7

%%%%%%%%%%%%%%%%%%%%%%%%% DFT图

figure(1);

k=1:N; subplot(2,2,1);stem(k,xk);title('\fontsize{16}原序列x（n)的DFT'); ht=get(gca,'title');set(ht,'color',[1 0 1]);

k=1:N; subplot(2,2,2);stem(k,yk1);title('\fontsize{16}y1(n)的DFT');