数字信号大作业.docx

数字信号大作业.docx

《数字信号大作业.docx》由会员分享，可在线阅读，更多相关《数字信号大作业.docx（38页珍藏版）》请在冰豆网上搜索。

数字信号大作业

学院：

电子工程学院

专业：

智能科学与技术

班级：

021151

学号：

02115096

姓名：

袁明

实验一：

信号、系统及系统响应

1、实验目的

（1）熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。

（2）熟悉时域离散系统的时域特性。

（3）利用卷积方法观察分析系统的时域特性。

（4）掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。

2、实验原理与方法

（1）时域采样。

（2）LTI系统的输入输出关系。

3、实验内容及步骤

（1）认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。

（2）编制实验用主程序及相应子程序。

①信号产生子程序，用于产生实验中要用到的下列信号序列：

a.xa（t）=A*e^-at*sin（Ω0t）u（t）

b.单位脉冲序列：

xb（n）=δ（n）

c.矩形序列：

xc（n）=RN（n）,N=10

②系统单位脉冲响应序列产生子程序。

本实验要用到两种FIR系统。

a.ha（n）=R10（n）;

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

③有限长序列线性卷积子程序

用于完成两个给定长度的序列的卷积。

可以直接调用MATLAB语言中的卷积函数conv。

conv用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。

调用格式如下：

y=conv（x,h）

4、实验结果记录

①分析采样序列的特性。

a.取采样频率fs=1kHz,，即T=1ms。

b.改变采样频率，fs=300Hz，观察|X（e^jω）|的变化，并做记录（打印曲线）；进一步降低采样频率，fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录打印这时的|X（e^jω）|曲线。

程序代码如下：

clearall;

A=444.128;

a=50*sqrt

（2）*pi;

m=50*sqrt

（2）*pi;

N=50;

n=0:

N-1;

w=linspace（-pi,pi,1000）;%设置w的范围

fork=1:

10%改变采样频率，显示出不同的频度谱

fs=input（'请输入采样频率'）;%输入不同的采样频率

T=1/fs;

x=A*exp（-a*n*T）.*sin（m*n*T）;

figure（k）;

subplot（2,1,1）%显示出x的时域波形

stem（n,x,'.'）;

xlabel（'n'）;ylabel（'x（n）'）;

title（'幅度谱'）;

X=x*exp（-1i*n'*w）;%对x（n）做DTFT变换

subplot（2,1,2）

plot（w/pi,abs（X））;%绘制x1（n）的幅度谱

xlabel（'\omega/π'）;ylabel（'|H（e^j^\omega）|'）

title（'频度谱'）;

c=input（'是否改变采样频率继续观察Y/N:

','s'）;

if（strcmpi（c,'N'）==1）

break;

end

end

（1）频率为1000KHZ的波形

（2）频率为300KHZ的波形

（3）频率为200KHZ的波形

由图可知，当采样频率进一步降低时，主瓣宽度逐渐变宽，频率混叠现象也逐渐严重。

存在较明显的失真现象。

②时域离散信号、系统和系统响应分析。

a.观察信号xb（n）和系统hb（n）的时域和频域特性；利用线性卷积求信号xb（n）通过系统hb（n）的响应y（n），比较所求响应y（n）和hb（n）的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。

b.观察系统ha（n）对信号xc（n）的响应特性，改变N=5，观察结果。

程序代码如下：

N=10;

xbn=1;

xcn=ones（1,N）;

han=ones（1,N）;

hbn=[1,2.5,2.5,1];

yn=conv（xbn,hbn）;%计算信号xb（n）通过系统hb（n）的响应y（n）

n1=0:

length（yn）-1;

n2=0:

length（hbn）-1;

subplot（2,1,1）;stem（n1,yn,'.'）%绘制y（n）的时域特性曲线

xlabel（'n'）;ylabel（'y（n）'）

title（'y（n）的时域特性曲线'）

subplot（2,1,2）;stem（n2,hbn,'.'）%绘制hb（n）的时域特性曲线

xlabel（'n'）;ylabel（'hb（n）'）

title（'hb（n）的时域特性曲线'）

figure

（2）%用于分析y（n）的频域特性

w=linspace（-pi,pi,1000）;%限制w的范围

Y=yn*exp（-1i*n1'*w）;%对y（n）做DTFT变换

subplot（2,1,1）;

plot（w/pi,abs（Y））;%绘制y（n）的幅度曲线

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'DTFT[y（n）]的幅度'）;

Hb=hbn*exp（-1i*n2'*w）;%对hb（n）做DTFT变换%用于分析hb（n）的频域特性

subplot（2,1,2）;

plot（w/pi,abs（Hb））;%绘制hb（n）的幅度曲线

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'DTFT[hb（n）]的幅度'）;

fork=1:

2

N=10;

xcn=ones（1,N）;

han=ones（1,N）;

figure（2+k）

zn=conv（xcn,han）;%观察系统ha（n）对信号xc（n）的响应特性

n3=0:

length（zn）-1;

w=linspace（-pi,pi,1000）;%限制w的范围

Z=zn*exp（-1i*n3'*w）;%对z（n）做DTFT变换

plot（w/pi,abs（Z））;%绘制z（n）的幅度曲线

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'DTFT[zn]的幅度'）;

End

xb（n）和系统hb（n）的时域和频域特性：

y（n）时域及频域特性：

N=10时。

Y（n）的时域频域特性曲线

N=5时。

Y（n）的时域频域特性曲线

时域离散信号、系统和系统响应分析。

实验代码：

clearall;

A=1;a=0.4;

m=2.0734;

N=50;T=1;

n1=0:

N-1;

w=linspace（-pi*2,pi*2,1000）;%设置w的范围

xn=A*exp（-a*n1*T）.*sin（m*n1*T）;

hbn=[1,2.5,2.5,1];

yn=conv（xn,hbn）;%计算信号xb（n）通过系统hb（n）的响应y（n）

n2=0:

length（hbn）-1;

n3=0:

length（yn）-1;

X=xn*exp（-1i*n1'*w）;%对x（n）做DTFT变换

H=hbn*exp（-1i*n2'*w）;%对hb（n）做DTFT变换

Y1=yn*exp（-1i*n3'*w）;%对y（n）做DTFT变换

Y2=X.*H;

figure

（1）;

subplot（2,1,1）;stem（n1,xn,'.'）%绘制y（n）的时域特性曲线

xlabel（'n'）;ylabel（'x（n）'）

title（'x（n）的时域特性曲线'）

subplot（2,1,2）;

plot（w/pi,abs（X））;%绘制y（n）的幅度曲线

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'DTFT[x（n）]的幅度'）;

figure

（2）;

subplot（2,1,1）;stem（n2,hbn,'.'）%绘制hb（n）的时域特性曲线

xlabel（'n'）;ylabel（'hb（n）'）

title（'hb（n）的时域特性曲线'）

subplot（2,1,2）;

plot（w/pi,abs（H））;%绘制y（n）的幅度曲线

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'DTFT[hb（n）]的幅度'）;

figure;

subplot（2,1,1）;%绘制z（n）的幅度和相位曲线

plot（w/pi,abs（Y1））;

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'DTFT[yn]的幅度'）;

subplot（2,1,2）;%绘制z（n）的幅度和相位曲线

plot（w/pi,abs（Y2））;

xlabel（'\omega/π'）;ylabel（'幅度'）

title（'Xn*Hn的幅度'）;

5、实验思考题分析

（1）在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？

它们所对应的模拟频率是否相同？

为什么？

答：

数字频率度量不相同，但他们所对应的模拟频率相同。

由w=Ω*Ts得，采样间隔变化时模拟频率对应的数字频率会有相应的变化，故其度量会有所变化。

（2）在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？

答：

有差异，所到的结果点数不同。

实验二：

用FFT作谱分析

1、实验目的

（1）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。

（2）熟悉FFT算法原理和FFT子程序的应用。

（3）学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。

2、实验步骤

（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。

（2）复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。

（3）编制信号产生子程序，产生以下典型信号供谱分析用：

（4）编写主程序。

（5）按实验内容要求，上机实验，并写出实验报告。

3、实验结果记录

（1）对2中所给出的信号逐个进行谱分析。

 程序代码如下：

N=[8,16,32,64];

n=0:

100;

x1=ones（1,4）;

x2=[n（1:

4）+1,8-n（5:

8）];

x3=[4-n（1:

4）,n（5:

8）-3];

x4=cos（pi*n/4）;

x5=sin（pi*n/8）;

figure

（1）;

pstem（x1）;

figure

（2）;

pstem（x2）;

figure（3）;

pstem（x3）;

figure（4）;

pstem（x4）;

figure（5）;

pstem（x5）;

f=64;

T=1/f;

fork=2:

4

n=0:

N（k）-1;

x6=cos（pi*8*n*T）+cos（pi*16*n*T）+cos（20*pi*n*T）;

X6=fft（x6,N（k））;

figure;

subplot（1,2,1）;

stem（0:

length（x6）-1,x6,'.'）;

title（'时域波形图'）;

subplot（1,2,2）;

stem（（0:

N（k）-1）*2/N（k）,abs（X6）,'.'）;

title（'频域波形图'）;

end

fork=1:

2

n=0:

N（k）-1;

x4=cos（pi*n/4）;

x5=sin（pi*n/8）;

x7=x4+x5;

X7=fft（x7,N（k））;

figure;

subplot（1,2,1）;

stem（0:

length（x7）-1,x7,'.'）;

title（'时域波形图'）;

subplot（1,2,2）;

stem（（0:

length（X7）-1）*2/length（X7）,abs（X7）,'.'）;

title（'频域波形图'）;

end

x8=x4+x5*j;

X81=fft（x8,N

（1））;

X82=fft（x8,N

（2））;

figure;

subplot（3,1,1）;

stem（0:

length（x8）-1,abs（x8）,'.'）;

title（'时域波形图'）;

subplot（3,1,2）;

stem（（0:

length（X81）-1）*2/length（X81）,abs（X81）,'.'）;

axis（[0,2,0,10]）;

title（'频域波形图（8）'）;

subplot（3,1,3）;

stem（（0:

length（X82）-1）*2/length（X82）,abs（X82）,'.'）;

title（'频域波形图（16）'）;

FFT变换及画图小程序

functionpstem（x）

N=[8,16];

if（length（x）>N

（1））

X1=fft（x（1:

N

（1））,N

（1））;

else

X1=fft（x,N

（1））;

end

if（length（x）>N

（2））

X2=fft（x（1:

N

（2））,N

（2））;

else

X2=fft（x,N

（2））;

end

subplot（2,2,1）;

if（length（x）>N

（1））

stem（0:

N

（1）-1,x（1:

N

（1））,'.'）;

else

stem（0:

length（x）-1,x,'.'）;

end

title（'x（n）时域波形图'）;

subplot（2,2,2）;

stem（（0:

N

（1）-1）*2/N

（1）,abs（X1）,'.'）;

axis（[0,2,0,max（abs（X1））]）;

title（'X频域波形图（8）'）;

subplot（2,2,3）;

if（length（x）>N

（2））

stem（0:

N

（2）-1,x（1:

N

（2））,'.'）;

else

stem（0:

length（x）-1,x,'.'）;

end

title（'x（n）时域波形图'）;

subplot（2,2,4）;

stem（（0:

N

（2）-1）*2/N

（2）,abs（X2）,'.'）;

title（'X频域波形图（16）'）;

（1）X1（n）和8点，16点TTF变换

（2）X2（n）和8点，16点TTF变换

（3）X3（n）和8点，16点TTF变换

（4）X4（n）和8点，16点TTF变换

（5）X5（n）和8点，16点TTF变换

（6）X6（n）和16,32,64点TTF变换

（7）X（7n）=x4（n）+x5（n）的8点和16点fft变换

（8）X8=x4（n）+i*x5（n）的8点和16点fft变换

4、实验思考题分析

（1）在N=8时，x2（n）和x3（n）的幅频特性会相同吗?

为什么?

N=16呢?

答：

N=8时幅频特性一样，N=16时幅频特性不一样。

（2）如果周期信号的周期预先不知道，如何用FFT进行谱分析?

答：

设一个定长的m值,先取2m,看2m/m的误差是否大，如大的话再取4m,看4m/2m的误差是否大，如不大，4m（4倍的m值）则可近似原来点的谱分析。

实验三用双线性变换法设计IIR数字滤波器

一、实验目的

1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法；

2、掌握数字滤波器的计算机仿真方法；

3、通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。

二、实验内容及原理

1、用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。

设计指标参数为：

在通带内截止频率低于0.2

时，最大衰减小于1dB；在阻带内[0.3

]频率区间上，最小衰减大于15dB。

2、以0.02

为采样间隔，打印出数字滤波器在频率区间[0,

/2]上的幅频响应特性曲线。

3、用所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

三、实验步骤

（1）：

复习有关巴特沃斯模拟滤波器的设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数H（z）;

（2）:

编写滤波器仿真程序，计算H（z）对心电图信号采样序列x（n）的相应序列y（n）,

（3）：

在通用计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验

（2）和（3）

四、心电图信号采样序列x（n）

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。

下面给出一实际心电图信号采样序列样本x（n），其中存在高频干扰。

在实验中，以x（n）作为输入序列，滤除其中的干扰成分。

x（n）={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,

-4,-4,-6,-6,-2,6,12,8,0,-16,

-38,-60,-84,-90,-66,-32,-4,-2,-4,8,

12,12,10,6,6,6,4,0,0,0,

0,0,-2,-4,0,0,0,-2,-2,0,

0,-2,-2,-2,-2,0}

五、实验结果及代码：

clearall;

%输入并显示滤波前心电图信号采样序列

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,...

0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,...

4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

figure;

subplot（1,2,1）;

n=0:

55;

stem（n,x,'.'）;

axis（[0,56,-100,50]）;

xlabel（'n'）;

ylabel（'x（n）'）;

title（'滤波前心电图信号采样序列'）;

gridon;

%利用所设计出的滤波器的差分方程，对心电图信号采样序列进行滤波处理

%绘制滤波后心电图采样序列波形

B=[0.09036,0.09036*2,0.09036];

A1=[1,-1.2686,0.7051];

A2=[1,-1.0106,0.3583];

A3=[1,-0.9044,0.2155];

y1=filter（B,A1,x）;

y2=filter（B,A2,y1）;

y3=filter（B,A3,y2）;

subplot（1,2,2）;

n=0:

55;

stem（n,y3,'.'）;

axis（[0,56,-100,50]）;

xlabel（'n'）;

ylabel（'y（n）'）;

gridon;

title（'滤波后心电图信号采样序列'）;

%绘制滤波器的幅频响应曲线

m=200;

n=0:

m-1;

[h1,w]=freqz（B,A1,m）;

[h2,w]=freqz（B,A2,m）;

[h3,w]=freqz（B,A3,m）;

h4=h1.*h2;

h4=h4.*h3;

figure;

subplot（1,2,1）

plot（n/m,abs（h4））;

title（'幅频响应特性曲线'）;

xlabel（'\omega\\pi'）;

ylabel（'幅度'）;

gridon;

axis（[0,1,0,1.1]）;

subplot（1,2,2）

plot（n/m,20*log10（abs（h4）））;

xlabel（'\omega\\pi'）;

ylabel（'幅度/db'）;

title（'损耗函数曲线'）;

axis（[0,0.5,-50,10]）;

gridon;

六、思考题

用双线性变换法设计数字滤波器过程中，变换公式

中T的取值，对设计结果有无影响？

为什么？

答：

没有影响。

因为双线性变换法不存在频率混叠现象，故T可以任意选择，但一般选T等于1，便于结果计算

实验四：

用窗函数法设计FIR数字滤波器

一、实验题目

用窗函数法设计FIR数字滤波器

二、实验目的

（1）掌握用窗函数法设计FIR数字滤波器的原理和方法。

（2）熟悉线性相位FIR数字滤波器特性。

（3）了解各种窗函数对滤波特性的影响。

三、实验原理

要设计出的滤波器的理想频率响应函数为

，则对应的单位脉冲响应为

窗函数设计法的基本原理是用有限长单位脉冲响应序列h（n）逼近hd（n）。

由于hd（n）往往是无限长序列，且是非因果的，所以用窗函数ω（n）将hd（n）截断，并进行加权处理，得到

h（n）就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数

为

式中，N为所选窗函数ω（n）的长度。

用窗函数法设计的滤波器性能取决于窗函数ω（n）的类型及窗口长度N的取值。

设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。

各种类型的窗函数可达到的阻带最小衰减和过渡带宽度可参见课本。

这样选定窗函数类型和长度N后，求出单位脉冲响应

，并求出

。

是否满足要求，要进行验算。

一般在h（n）尾部加零使长度满足2的整数次幂，以便用FFT计算

。

如果要观察细节，补零点数增多即可。

如果

不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。

如果要求线性相位特性，则h（n）还必须满足：

h（n）=+/-h（N-1-n），根据式中的正负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。

要根据所设计的滤波特性正确选择其中一类。

四、实验内容及步骤

（1）复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。

（2）编写程序

①编写能产生矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序。

②编写主程序。

主程序框图如下图所示，仅供参考。

其中幅度特性要求用dB表示。

（3）上机实验内容

a、用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率Wc=rad。

窗口N=15，33。

要求在两种窗口长度情况下，分别求出h（n），打印出相应的幅频特性和相频曲线。

观察3db带宽和20db带宽，总结窗口N对滤波特性的影响。

设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数即

b、N=33，Wc=rad，用四种窗函数设计线性相位低通滤波器。

绘制相应的幅频特性曲线，观察3dB带宽和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波特性的影响。

五、实验程序及波形

a:

N=15,33

N=15;wc=pi/4;

a=（N-1）/2;

n=0:

（N-1）;

m=n-a+eps;

h