数字信号处理实验与课程设计课案.docx

数字信号处理实验与课程设计课案.docx

《数字信号处理实验与课程设计课案.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验与课程设计课案.docx（64页珍藏版）》请在冰豆网上搜索。

数字信号处理实验与课程设计课案

数字信号处理实验与课程设计教程

戴虹编

工学部计算机与信息工程学院

2015年12月

内容简介:

数字信号处理是高校通信工程、电子信息工程、自动化、测控等专业的一门非常重要的专业基础课程，是我校的校重点课程、校精品课程与上海市教委重点课程。

本书是数字信号处理课程的实验及课程设计教材，基于matlab编程，配合理论授课内容与应用型人才培养的教学目标而编写。

全书分为上下两篇，上篇为数字信号实验指导，共有5个实验：

信号、系统及系统响应、离散信号与系统的复频域分析、DFT和FFT频谱分析、IIR数字滤波器设计、FIR数字滤波器设计。

下篇是数字信号处理课程设计教程,共3个课程设计题:

双音多频（DTMF）的产生与检测、音乐信号中啸叫噪声的消除、回音消除。

力求使学生对数字信号处理课程产生浓厚的兴趣，并掌握数字信号处理的matlab仿真实践方法。

课程设计以小组方式进行，要求学生任选一题进行设计并撰写课程设计报告。

本书适合作为普通高校电子信息工程、通信工程、测控、自动化等专业本科数字信号处理实验与课程设计教材，也可作为相关领域工程技术人员的参考书籍。

上篇数字信号处理实验指导

实验一信号、系统及系统响应

实验二离散信号与系统的复频域分析

实验三DFT和FFT频谱分析

实验四IIR数字滤波器设计

实验五FIR数字滤波器设计

下篇数字信号处理课程设计教程

课题1双音多频（DTMF）信号的产生与检测

课题2音乐信号中的啸叫噪声消除

课题3回音消除

附录AMatlab程序设计入门

附录B实验报告撰写格式

参考文献

上篇数字信号处理实验指导

实验一信号、系统及系统响应

一、实验目的

1.掌握典型序列的产生方法。

2.掌握DFT的实现方法，利用DFT对信号进行频域分析。

3.熟悉连续信号经采样前后频谱的变化，加深对时域采样定理的理解。

4.分别利用卷积和DFT分析信号及系统的时域和频域特性，验证时域卷积定理。

二、实验环境

1．Windows2000操作系统

2．MATLAB6.0

三、实验原理

1.信号采样

对连续信号xa（t）=Ae-atsin（Ω0t）u（t）进行采样，采样周期为T,采样点0≤n<50，得采样序列xa（n）=。

2.离散傅里叶变换（DFT）

设序列为x（n）,长度为N,则

X（ejωk）=DFT[x（n）]=

x（n）e-jωkn,

其中ωk=

（k=0,1,2,…,M-1）,通常M>N,以便观察频谱的细节。

|X（ejωk）|----x（n）的幅频谱。

4.连续信号采样前后频谱的变化

a（jΩ）=

即采样信号的频谱

a（jΩ）是原连续信号xa（t）的频谱Xa（jΩ）沿频率轴，以周期重复出现，幅度为原来的倍。

5.采样定理

由采样信号无失真地恢复原连续信号的条件，即采样定理为：

。

6．时域卷积定理

设离散线性时不变系统输入信号为x（n）,单位脉冲响应为h（n）,则输出信号y（n）=；由时域卷积定理，在频域中，

Y（ejω）=FT[y（n）]=。

四、实验内容

1.分析采样序列特性

（1）程序输入

产生采样序列xa（n）=Ae-anTsin（Ω0nT）u（n）（0≤n<50）,其中A=44.128,a=50

Ω0=50

采样频率fs（可变），T=1/fs。

（要求写%程序注释）

%程序shiyan11.m

clear%

clc%

A=444.128;

a=50*sqrt

（2）*pi;%

w0=50*sqrt

（2）*pi;

fs=input（'输入采样频率fs='）;

T=1/fs;

N=50;

n=0:

N-1;

xa=A*exp（-a*n*T）.*sin（w0*n*T）;%

subplot（221）;stem（n,xa,'.'）;grid;%

M=100;

[Xa,wk]=DFT（xa,M）;%

f=wk*fs/（2*pi）;%

subplot（222）;plot（f,abs（Xa））;grid;%

DFT子函数:

DFT.m

function[X,wk]=DFT（x,M）

N=length（x）;%

n=0:

N-1;

fork=0:

M-1

wk（k+1）=2*pi/M*k;

X（k+1）=sum（x.*exp（-j*wk（k+1）*n））;%

end

（2）实验及结果分析

a.取fs=1000（Hz）,绘出xa（n）及|Xa（ejωk）|的波形。

b.取fs=300（Hz）,绘出xa（n）及|Xa（ejωk）|的波形。

c.取fs=200（Hz）,绘出xa（n）及|Xa（ejωk）|的波形。

d.a,b,c中，哪几种情况出现了频谱混叠现象？

；出现频谱混叠的原因是

。

2.时域离散信号和系统响应分析

（1）hb（n）=δ（n）+2.5δ（n-1）+2.5δ（n-2）+δ（n-3）

程序语句为hb=[1,2.5,2.5,1];

（2）卷积语句：

y=conv（x,h）

其中x----输入序列x（n）；h----单位脉冲响应h（n）;y-------输出序列y（n）。

3.卷积定理验证

（1）编程实现y（n）=xa（n）*hb（n）,其中

xa（n）=Ae-anTsin（Ω0nT）u（n）（0≤n<50）,A=1,a=0.4,Ω0=2.0734,T=1

hb（n）=δ（n）+2.5δ（n-1）+2.5δ（n-2）+δ（n-3）

及Y（ejωk）=DFT[y（n）]（M=100），（利用DFT.m）；绘出|Y（ejωk）|波形。

（2）编程实现Xa（ejωk）=DFT[Xa（n）]（M=100）及Hb（ejωk）=DFT[hb（n）]（M=100）；

计算Y（ejωk）=Xa（ejωk）Hb（ejωk）；绘出|Y（ejωk）|波形。

问：

（1）和

（2）中的|Y（ejωk）|波形一致吗？

；为什么？

。

4.正弦序列的产生

设　正弦序列x（n）=sin（ωn）,取样频率fs=64Hz,信号频率f=5Hz,

样点n=0~N-1,（N=64）,ω=2πf/fs,编程产生x（n）并绘图。

程序：

%shiyan141.m

clear

clc

N=64;fs=64;f=5;%

n=0:

1:

N-1;%

w=2*pi*f/fs;%

x=sin（w*n）;%

stem（x,'.'）%

实验要求:

1）在%后的空格内填入注释2）运行上述程序，绘出结果波形。

3）设双音多频信号"2"为x（n）=sin（ω1n）+sin（ω2n）,其中f1=697Hz,f2=1336Hz,

取样频率fs=8000Hz,编程产生x（n）并绘出x（n）的波形。

（n=0~799）

5.拓展题

数字信号处理算法之一----时间平均

时间平均可用来消除周期信号中的随机噪声。

对m个周期的振动信号x（n）=s（n）+u（n）[其中s（n）--故障信号，为余弦序列；

u（n）---随机噪声]做时间平均,即把x（n）按周期n分段，将每个周期的对应点相加再做平均，计算时间平均前后的信噪比。

%时间平均程序　program11.m　

clear

clc

n=10;%定义n=10

m=input（'m='）;%输入周期数m

loadsip%载入sip

x=sip;

s=zeros（1,n）;

fori=1:

m%定义i从1到m

s=s+x（1+n*（i-1）:

i*n）;

end

s=s/m;

k=[0:

n-1];

subplot（321）;plot（k,x（1:

n））;grid;

subplot（322）;plot（k,s）;grid;

i=0:

1:

n-1;

s0=cos（2*pi*i/n）;

ps=sum（s0.^2）/n;

pu=1;

snr0=10*log10（ps/pu）%原信噪比

py=sum（（s-s0）.^2）/n;

snr=10*log10（ps/py）%时间平均后的信噪比

%数据文件sip的产生程序　shu.m

clear

clc

n=10;m=1000;

i=0:

1:

n-1;

s=cos（2*pi*i/n）;

x=zeros（1,n*m）;

u=randn（size（1:

n*m））;%生成1到n*m的矩阵

forj=1:

m%定义j从1到m

x（1+n*（j-1）:

j*n）=s+u（1+n*（j-1）:

j*n）;

end

savesipx-ascii%产生数据文件sip

实验要求：

1）在空格中填入注释。

2）先运行shu.m，产生数据文件sip,再运行program11.m,分别绘出当m=10,100,500,1000时，输入／输出信号的波形。

实验二离散信号与系统的复频域分析

一、实验目的

1.掌握求解离散时间系统差分方程的两种方法：

迭代法和filter函数法。

2.利用Z变换对系统进行复频域分析。

3.掌握系统零极点的绘制方式及Z反变换的求解方法。

4.熟悉Z变换的应用。

二、实验环境

1.Windowsxp操作系统

2.MATLAB2007a软件

三、实验原理

3.系统的输出y（n）与输入x（n）及单位脉冲响应h（n）的关系:

4.一个N阶常系数线性差分方程表达式:

5.系统函数H（z）与X（z）,Y（z）的关系为：

H（z）=。

4.H（z）零点定义为：

；H（z）极点定义为：

。

5.本实验涉及的Matlab函数

1）求解差分方程：

y=filter（b,a,x）

其中：

x---输入信号；y---差分方程的解。

b,a----输入/输出信号前的系数。

2）画零极点图:

zplane（b,a,n）

其中：

b,a----H（z）分子/分母系数。

3）求z反变换

[r,p,k]=residuez（b,a）

b---H（z）的分子多项式系数；a---H（z）的分母多项式系数；

输出：

r,p,k的含义如下：

“residuez”可将H（z）分解为简单有理分式之和：

式<1>中p

（1）,p

（2）…p（n）的集合为列向量p;r

（1）,r

（2）…r（n）的集合为列向量r;

k

（1）,k

（2）…的集合为行向量k;

已知r,p,k即可手工写出h（n）的表达式：

h（n）={r

（1）[p

（1）]n+r

（2）[p

（2）]n+…r（n）[p（n）]n}u（n）+k

（1）δ（n）+k

（2）δ（n-1）+……

四、实验内容

1.已知：

系统的差分方程为：

y（n）-y（n-1）+0.9y（n-2）=x（n）

（1）分别利用迭代法和filter函数法求此系统的单位脉冲响应h（n）,并绘图。

[1]迭代法：

%程序：

program11.m

clear

clc

N=100;

x=[1zeros（1,N-1）];%产生x（n）=δ（n）

y

（1）=0;%y（-2）=0

y

（2）=0;%y（-1）=0

y（3）=1;%y（0）=1

forn=4:

N%迭代法求解该差分方程,得y（n）

y（n）=x（n）+y（n-1）-0.9*y（n-2）;

end

k=-2:

N-3;%时间轴范围k

stem（k,y,'.'）;%画y（n）=h（n）

title（'系统的单位脉冲响应h（n）'）;%写标题

h（n）波形：

[2]filter函数法：

%程序：

program12.m

clear

clc

N=100;

x=[1zeros（1,N-1）];%产生x（n）=δ（n）

b=[1];%差分方程系数b,a

a=[1,-1,0.9];

h=filter（b,a,x）;%用filter函数求解系统的单位脉冲响应h（n）

k=-2:

N-3;

stem（k,h,'.'）;title（'系统的单位脉冲响应h（n）'）;

h（n）波形：

2）修改输入信号为x（n）=u（n）,编程实现：

求系统的阶跃响应y1（n）。

[u（n）由：

u=ones（1,n）产生]。

%程序：

program13.m

y1（n）波形：

3）利用由1）产生的h（n）,利用卷积法产生系统的阶跃响应：

y2（n）=u（n）*h（n）。

%程序：

program14.m

y2（n）波形：

问：

y1（n）和y2（n）是否相同？

2.求解差分方程（迭代法及filter函数法），实现一个数字振荡器（正弦波信号）：

h（n）=sin（ωkn）,其中ωk=2πfk/fs,频率fk=1000Hz,取样频率fs=10kHz。

[提示：

h（n）对应的差分方程为：

h（n）-ah（n-1）+h（n-2）=bδ（n-1）<1>

a=2cos（ωk）=1.618,b=sin（ωk）=0.588]

画出h（n）并用“zplane”画该系统的零极点图。

%迭代法实现程序program21.m

h（n）波形：

%filter函数法实现程序program22.m

h（n）波形：

问：

极点位置处于单位圆（内/上/外）？

3.求下列序列的z变换并用”zplane”画出零极点分布图。

1）x（n）={1,2,1,3}，该序列的z变换为：

X（z）=1+2z-1+z-2+3z-3

%程序program31.m

clear

clc

b=[1213];%H（z）分子系数

a=[1000];%H（z）分母系数

zplane（b,a）;%画零极点分布图

结果图形：

2）x（n）=0.2n[u（n）-u（n-5）],则:

X（z）=

%程序program32.m

结果图形：

4.z反变换

用”residuez”函数求：

%程序program4.m

clear

clc

b=[4-21180];%H（z）分子系数

a=[1-611-6];%H（z）分母系数

[r,p,k]=residuez（b,a）%求z反变换中的系数

运行该程序，得：

r=

p=

k=

则h（n）=

实验三DFT和FFT频谱分析

一、实验目的

1.掌握DFT频谱分析的原理与编程方法。

2.理解FFT算法的编程思想。

2.熟练掌握利用FFT对信号作频谱分析，包括正确地进行参数选择、画频谱

及读频谱图。

3.利用FFT频谱分析进行快速卷积和太阳黑子周期性检测。

二、实验环境

1.Windowsxp以上操作系统

2.安装MATLAB2007a软件

三、实验原理

1.离散傅里叶变换（DFT）

设序列为x（n）,长度为N,则

X（ejωk）=DFT[x（n）]=

x（n）e-jωkn,

其中ωk=

（k=0,1,2,…,M-1）,通常M>N,以便观察频谱的细节。

|X（ejωk）|----x（n）的幅频谱。

2.谱分析参数选择

1）设信号x（t）最高频率为fc,对其进行取样得x（n）,根据取样定理，取样频率fs必须满足：

。

2）设谱分辨率为F,则最小记录时间tpmin=；取样点数

N≥;为使用快速傅里叶变换（FFT）进行谱分析，N还须满足：

N=。

3.用FFT计算信号x（n）的频谱。

[设x（n）为实信号]

快速傅里叶变换（FFT）是DFT的一种快速算法，其使得DFT的运算速度大为加快。

1）对信号x（n）作N点FFT,得频谱X（k）（k=0~N-1）

X（k）=XR（k）+jXI（k）（k=0~N/2-1），XR（k）—X（k）的实部；XI（k）—X（k）的虚部。

Matlab语句：

Y=fft（x,N）

其中：

x----x（n）;Y----X（k）

2）幅频谱:

|X（k）|=，由于x（n）为实信号，因此|X（k）|对称，

Matlab语句：

abs（Y）

iii）功率谱:

PSD（k）=|X（k）|2/N=X（k）X*（k）/N

Matlab语句：

PSD=Y.*conj（Y）/N

其中：

conj（Y）--X*（k）[X（k）的共轭]

4.读频谱图

频谱图中任意频率点k对应实际频率为:

fk=。

5.用FFT实现线性卷积运算

用FFT实现y（n）=x（n）*h（n）的步骤为：

1）设x（n）及h（n）的长度分别为N1和N2。

为使循环卷积等于线性卷积，用补0的方

法使x（n）,h（n）长度均为N，则N须满足N≥；为用FFT计算DFT,则N

还须满足N=。

2）用FFT计算X（k）,H（k）（N点）。

3）Y（k）=;y（n）=。

四、实验内容

1.根据公式设计DFT原理程序，并计算：

x（n）=[1,1,1,1]的4,16,64点DFT并绘图。

%DFT/IDFT程序DFT.m

clc

clear

xn=input（'x（n）='）;%输入序列x（n）=[1111]

M=length（xn）;%x（n）的长度M

N=input（'变换区间N='）;%变换区间N

xn=[xnzeros（1,N-M）];%补0,使xn长度为N

n=0:

N-1;k=0:

N-1;

nk=n'*k;

wn=exp（-j*2*pi/N）;%旋转因子wn

wnK=wn.^nk;

xk=xn*wnK%作x（n）的DFT=xk

subplot（211）;stem（k,abs（xk）,'.'）;gridon;

%显示xk的幅频谱（离散曲线）

subplot（212）;plot（k,abs（xk））;gridon;

%显示xk的幅频谱（连续曲线）

运行结果：

问：

由此得出怎样的结论？

2.理解DIT-FFT算法原理程序，并用它计算X（k）=FFT[R4（n）],分别取N=4,8,16

和64,绘出幅频谱|X（k）|。

%程序DIT.m

clear

clc

x=input（'x='）;%

N=input（'N='）;%

x（length（x）+1:

N）=zeros（1,N-length（x））;%补0x（1:

N）

l=log2（N）;

x1=zeros（1,N）;

forj1=1:

N%倒序

x1（j1）=x（bin2dec（fliplr（dec2bin（j1-1,l）））+1）;

end

%%FFT（DIT）%%

M=2;

while（M<=N）

W=exp（-2*j*pi/M）;%旋转因子W

V=1;

fork=0:

1:

M/2-1%k为每级蝶形运算旋转因子的个数

fori=0:

M:

N-1%i为各群的首序号

p=k+i;

q=p+M/2;

A=x1（p+1）;

B=x1（q+1）*V;

x1（p+1）=A+B;%本级蝶形运算，x1最终存放X（k）

x1（q+1）=A-B;

end

V=V*W;%旋转因子W的变化

end

M=2*M;%第M级

end

%%%%%%

subplot（211）;stem（x,'.'）;gridon;%

title（'x（n）'）;%

subplot（212）;stem（abs（x1）,'.'）;gridon;%

title（'|X（k）|'）;%

3.FFT谱分析

设信号为x（t）=sin（2πf1t）+sin（2πf2t）+随机噪声，f1=50Hz,f2=120Hz,以取样频率

fs=1kHz对x（t）进行取样,样本长度tp=0.25s,得x（n）,对x（n）作256点FFT,得频谱X（k）,画原信号x（n）,幅频谱|X（k）|以及功率谱PSD（k）,对信号进行谱分析。

%程序pufenxi.m

clear

clc

fs=1000;

t=0:

1/fs:

0.25;%

N=256;%

f1=50;f2=120;%

s=sin（2*pi*f1*t）+sin（2*pi*f2*t）;%

x=s+randn（size（t））;%信号+噪声x（n）

Y=fft（x,N）;%

PSD=Y.*conj（Y）/N;%

f=fs/N*（0:

N/2-1）;%

subplot（311）;plot（x）;%

subplot（312）;plot（f,abs（Y（1:

N/2）））;%

subplot（313）;plot（f,PSD（1:

N/2））;%

1）画出图形窗口显示的图形，并注名每个图形的含义。

2）回答下列问题：

i）观察幅频谱图，可以发现，信号x（n）含有的两个频率分量分别是

Hz和Hz。

ii）在该程序中的“f=fs/N*（0:

N/2-1）;”下添加“k=0:

N/2-1;”，

“plot（f,abs（Y（1:

N/2）））;”改为“plot（k,abs（Y（1:

N/2）））;”

重新运行该程序并观察幅频谱图，图中两峰值对应的下标分别是

和。

它们的含义为。

再将该程序中的N改为512,重新运行该程序并观察幅频谱图，这时图中

两峰值对应的下标分别是和。

结果是否和上面的相同？

为什么？

。

iii）本例的频谱分辨率F是Hz,改变f2=60Hz,问：

在幅频谱中，能否分辨f1和f2

对应的频率分量？

。

为什么？

。

再改变f2=52Hz,问：

在幅频谱中，能否分辨f1和f2对应的频率分量？

。

为什么？

。

再改变f2=600Hz，在幅频谱中，f2对应的频率分量出现在Hz；

问：

在fs=1000Hz的情况下，能否正确检测f2对应的频率分量？

。

为什么？

。

为了正确检测f2对应的频率分量，则fs至少取多少Hz?

Hz。

在该程序中改变fs，验证你的结论。

iv）比较幅频谱和功率谱，可以发现功率谱具有的特性。

4.FFT实现任意两个序列的快速卷积。

%程序fftjuanji.m

clear

clc

x1=input（'x1='）;x2=input（'x2='）;%

N1=length（x1）;N2=length（x2）;%序列x1（n）,x2（n）的长度

E=ceil（log2（N1+N2-1））;%ceil---向+∞方向取整

N=2^E;%

x1=[x1,zeros（1,N-N1）];%

x2=[x2,zeros（1,N-N2）];

X1=fft（x1,N）;%

X2=fft（x2,N）;

Y=X1.*X2;%

y=ifft（Y,N）%

结果分析：

1）回到MATLAB窗口，键入：

x1=[111],x2=[12],回车。

结果:

y=

2）问：

可用Matlab中的什么函数验算上述卷积结果？

5.利用谱分析观察太阳黑子周期性。

以100年中记录到的太阳黑子出现次数为信号x（n）,对x（n）作功率谱，从中观察太阳黑子周期性。

%程序taiyangheizi.m

clear

clc

x=[101826635317209215412585683823102483...

132131118906760474121166471434454348...

4228108201512143