数字信号处理实验报告

数字信号处理课程设计报告

通信11级16班

学号：52111604

姓名：薛焕福

实验一离散时间系统及离散卷积

一、实验目的

1、熟悉MATLAB软件的使用方法。

2、熟悉系统函数的零极点分布、单位脉冲响应和系统频率响应等概念。

3、利用MATLAB绘制系统函数的零极点分布图、系统频率响应和单位脉冲响应。

4、熟悉离散卷积的概念，并利用MATLAB计算离散卷积。

二、实验原理

本实验利用MATLAB的库函数，实验冲击响应零极点、幅度相位谱的绘制。熟悉MATLAB的使用。并利用离散卷积的概念编写一个计算离散卷积的程序。

三、实验内容及要求

1.离散时间系统的单位冲击响应

（1）选择一个离散时间系统；

（2）用笔进行差分方程的递推计算；

（3）编制差分方程的递推计算程序；

（4）在计算机上实现递推运算；

（5）将程序计算结果与笔算的计算结果进行比较，验证程序运行的正确性；

2.离散系统的幅频、相频的分析方法

（1）离散系统的幅频、相频的分析方法

（2）用笔计算几个特殊的幅频、相频的值，画出示意曲线图；

（3）编制离散系统的幅频、相频的分析程序；

（4）在计算机上进行离散系统的幅频、相频特性计算，并画出曲线；

（5）通过比较，验证程序的正确性；

3. 离散卷积的计算

（1）选择两个有限长序列，用笔计算其线性卷积；

（2）编制有限长序列线性卷积程序；

（3）利用计算程序对（1）选择的有限长序列进行卷积运算；

（4）比较结果验证程序的正确性。

四、实验过程

1、离散时间系统的单位脉冲响应程序:

a=[1,-1,0.91]; % 系统的系数

b=[1,2,3];

n0=0;

n=-20:120;

x=[(n-n0)==0]; % 利用filter函数进行滤波，从而求得系统的单位冲击响应

h=filter(b,a,x);

figure(1)

stem(n,h); % 利用stem绘制曲线图

title('冲击响应');

xlabel('n');

ylabel('h(n)');

h

figure(2)

[z,p,g]=tf2zp(b,a); % 求传递函数

zplane(z,p) % 求零极点

（2）实验结果：

2、离散系统的幅频、相频的分析程序：b=[0.6,1.2,0.4,0.9]; % 分子的系数

a=[1.0,-1.7,1.1,-0.5] ; % 分母的系数

m=0:length(b)-1;

l=0:length(a)-1;

K=500;

k=1:K;

w=pi*k/K;

H=(b*exp(-j*m'*w))./(a*exp(-j*l'*w)); % 求传递函数magH=abs(H); % 求传递函数的幅值

angH=angle(H); % 求传递函数的相位

figure(1)

subplot(2,1,1);

plot(w/pi,magH); % 定义坐标轴的区间grid;

axis([0,1,0,1]);

xlabel('w(pi)');

ylabel('|H|');

title('·幅度，相位响应');

subplot(2,1,2);

plot(w/pi,angH);

grid;

xlabel('w(pi)');

ylabel('angle(H)');

（2）实验结果：

3、离散卷积的计算程序：function pr13()

n=-5:50;

u1=stepseq(0,-5,50);

u2=stepseq(10,-5,50);

x=u1-u2;

h=((0.9).^n).*u1;

figure(1)

subplot(3,1,1);

stem(n,x);

axis([-5,50,0,2]);

title('输入序列');

xlabel('n');

ylabel('x(n)');

subplot(3,1,2);

stem(n,h);

axis([-5,50,0,2]);

title('冲击响应序列');

xlabel('n');

ylabel('h(n)');

[y,ny]=conv_m(x,n,h,n);

subplot(3,1,3);

stem(ny,y);

axis([-5,50,0,8]);

title('输出响应');

xlabel('n');

ylabel('y(n)');

function [x,n]=stepseq(n0,n1,n2) n=n1:n2;

x=[(n-n0)>=0];

function [y,ny]=conv_m(x,nx,h,nh) nyb=nx(1)+nh(1);

nye=nx(length(x))+nh(length(h)); ny=[nyb:nye];

y=conv(x,h);

（2）实验结果：

实验二离散傅立叶变换与快速傅立叶变换

一、实验目的：

1、加深理解离散傅立叶变换及快速傅立叶变换概念；

2、学会应用FFT对典型信号进行频谱分析的方法；

3、研究如何利用FFT程序分析确定性时间连续信号；

4、熟悉应用FFT实现两个序列的线性卷积的方法。

二、实验原理

在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为

反变换为

有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。

用FFT计算线性卷积：

用FFT可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度

N≥N1＋N2

对于长度不足N的两个序列，分别将他们补零延长到N。

当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。通常，我们选择重叠相加法。

重叠相加法：将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。

三、实验内容：

1、用离散傅立叶变换程序处理时间抽样信号，并根据实序列离散傅立叶变换的对称性，初步判定程序的正确性；

2、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列

和的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。

三角波序列：

反三角波序列：

3、已知余弦信号如下

当信号频率，采样间隔，采样长度时，对该信号进行傅立叶变换。用FFT程序分析正弦信号，分别在以下情况下进行，并且分析比较结果

（1）F=50, N=32,T=0.000625；

（2）F=50, N=32,T=0.005；

（3）F=50, N=32,T=0.0046875；

（4）F=50, N=32,T=0.004；

（5）F=50, N=64=0.000625

4、选定某一时间信号进行N=64点离散傅立叶变换，并且，对同一信号进行快速傅立

叶变换，并比较它们的速度。

四、实验要求：

1、调试实验程序，并且，给参考程序加注释；

2、完成实验内容2，并对结果进行分析。实验中的信号序列和，在单位

圆上的Z变换频谱和会相同吗？如果不同，你能说出哪一个低频分量更多一些吗？为什么？

3、完成实验内容3，并对结果进行分析；

4、利用编制的计算卷积的计算程序，分别给出一下三组函数的卷积结果。

（1）

（2）

（3）

五、实验过程：

1、三角波与反三角波序列的时域和频域特性：

（1）实验程序：

clear all;

clc;

N=64;

n=1:N;

xc=zeros(1,N);

xd=zeros(1,N);

%生成三角波和反三角波

for i=1:4

xc(i+28)=i;

end

for j=5:8

xc(j+28)=9-j;

end

for i=1:4

xd(i+28)=5-i;

end

for j=5:8

xd(j+28)=j-4;

end

figure(1);

subplot(2,1,1),stem(n,xc);

title('三角波的时域波形');

subplot(2,1,2),stem(n,xd);

title('反三角波的时域波形'); %FFT变换

XC=fft(xc,N);

XD=fft(xd,N);

AXC=abs(XC);

AXD=abs(XD);

figure(2);

subplot(2,1,1),plot(AXC);

title('三角波的幅度谱波形'); subplot(2,1,2),plot(AXD);

title('反三角波的幅度谱波形'); （2）实验结果：

2、余弦信号的傅立叶变换：

N=64;F=50;T=0.000625;

n=0:N-1;

xn=cos(2*pi*n*F*T);

X=dft(xn,N);

magX=abs(X);

subplot(2,1,1),stem(n,xn);xlabel('n');title('离散序列x(n)'); subplot(2,1,2),stem(n,magX);xlabel('(pi)');title('DTFT[X]');

离散序列x(n)

n

DTFT[X]

(pi)（1）F=50,N=32,T=0.000625;

F=50,N=64,T=0.000625;

实验程序：

clear all;

clc;

%F=50, N=64,T=0.000625

F=50;

N=64;

T=0.000625;

n=1:N;

x=cos(2*pi*F*n*T);

figure(1);

subplot(2,2,1);

plot(n,x);

title('x(n)F=50, N=64,T=0.000625');

xlabel('n');

X=fft(x);

magX=abs(X);

subplot(2,2,3);

plot(n,X);

title('FFT|X|');

xlabel('f(pi)');

%F=50, N=32,T=0.000625

F=50;

N=32;

T=0.000625; n=1:N;

x=cos(2*pi*F*n*T);

subplot(2,2,2); plot(n,x);

title('x(n)F=50, N=32,T=0.000625'); xlabel('n');

X=fft(x);

magX=abs(X); subplot(2,2,4); plot(n,X); title('FFT|X|'); xlabel('f(pi)');

实验结果：

（2）F=50,N=32,T=0.005;

F=50,N=32,T0.0046875; 实验程序：

F=50; N=32; T=0.005; n=1:N;

0 20

40 60 80

01x(n)F=50, N=64,T=0.000625

n 0

20

40 60

80

0FFT|X|

f(pi)

0 10

20 30 40

x(n)F=50, N=32,T=0.000625

n 0

10

20 30

40

FFT|X|

f(pi)

x=cos(2*pi*F*n*T);

figure(2);

subplot(2,2,1);

plot(n,x);

title('x(n)F=50, N=32,T=0.005'); xlabel('n');

X=fft(x);

magX=abs(X);

subplot(2,2,3);

plot(n,X);

title('FFT|X|');

xlabel('f(pi)');

%F=50, N=32,T=0.0046875

F=50;

N=32;

T=0.0046875;

n=1:N;

x=cos(2*pi*F*n*T);

subplot(2,2,2);

plot(n,x);

title('x(n)F=50, N=32,T=0.0046875'); xlabel('n');

X=fft(x);

magX=abs(X);

subplot(2,2,4);

plot(n,X);

title('FFT|X|');

xlabel('f(pi)');

实验结果：

（3）F=50,N=32,T=0.004; 实验程序：

F=50; N=32; T=0.004; n=1:N;

x=cos(2*pi*F*n*T); figure(3);

subplot(2,1,1); plot(n,x);

title('x(n)F=50, N=32,T=0.004'); xlabel('n');

X=fft(x);

magX=abs(X); subplot(2,1,2); plot(n,X); title('FFT|X|'); xlabel('f(pi)');

实验结果：

10

20 30 40

01x(n)F=50, N=32,T=0.005

n 0

10

20 30 40

012-14

FFT|X|

f(pi)

10

20 30

40

x(n)F=50, N=32,T=0.0046875

n 0

10

20 30

40 FFT|X|

f(pi)

x(n)F=50, N=32,T=0.004

05101520253035

n

FFT|X|

05101520253035

f(pi)

3、三组卷积结果：

（1）N1=32;

n1=1:14;

x1n=ones(1,14);

h1n=(4/5).^n1;

figure(1);

X1=fft(x1n,N1);

H1=fft(h1n,N1);

Y1=X1.*H1;

y1=ifft(Y1);

stem(abs(y1));

title('第一组y(n)=x(n)*h(n)) N=32点的卷积');

0.511.522.533.54第一组 y(n)=x(n)*h(n)) N=32点的卷积

（2）

N2=32;

x2n([1:10])=1; n2=1:20;

h2n=0.5*sin(0.5*n2);

figure(2);

X2=fft(x2n,N2); H2=fft(h2n,N2); Y2=X2.*H2; y2=ifft(Y2);

stem(abs(y2));

title('第二组 y(n)=x(n)*h(n)) N=32点的卷积');

00.20.40.60.811.21.41.61.82第二组 y(n)=x(n)*h(n)) N=32点的卷积

（3） N1=32; n1=1:10;

x1n=(1-0.1*n1); h1n=0.1*n1;

figure(1);

X1=fft(x1n,N1); H1=fft(h1n,N1); Y1=X1.*H1; y1=ifft(Y1);

stem(abs(y1));

title('第三组 y(n)=x(n)*h(n)) N=32点的卷积');

第三组y(n)=x(n)*h(n)) N=32点的卷积

实验三IIR数字滤波器设计

一、实验目的

1、学习模拟－数字变换滤波器的设计方法；

2、掌握双线性变换数字滤波器设计方法；

3、掌握实现数字滤波器的具体方法。

二、实验内容：

1、设计一个巴特沃思数字低通滤波器，设计指标如下：

通带内幅度衰减不大于1dB;

阻带内幅度衰减不小于15dB；

2、编制计算设计的数字滤波器幅度特性和相位特性的程序，并进行实验验证。

编制实现该数字滤波器程序并且实现数字滤波

（1）分别让满足所设计的滤波器的通带、过渡带、阻带频率特性的正弦波通过滤波器，验证滤波器性能；

（2）改变正弦抽样时间，验证数字低通滤波器的模拟截止频率实抽样时间的函数。

三、实验要求：

1、编制实验内容要求的程序，并给程序加注释。

2、根据实验结果，给出自己设计的数字滤波器的幅度特性和相位特性。

3、用所设计的滤波器对不同频率的正弦波信号进行滤波，并说明其特性。

4、fp=0.2KHz,Rp=1dB,fs=0.3KHz,As=25dB,T=1ms;分别用脉冲响应不变法及双极性变换法设计一巴特沃斯数字低通滤波器，观察所设计的数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。

四、实验原理

1、脉冲响应不变法：

用数字滤波器的单位脉冲响应序列模仿模拟滤波器的冲激响应，让正好等于的采样值，即，其中T为采样间隔，如果以及分

别表示的拉氏变换及的Z变换，则

2、双线性变换法：

s平面与z平面之间满足以下映射关系：s平面的虚轴单值地映射于z 平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存

在混叠问题。双线性变换时一种非线性变换，这种非线性引起的幅频特性畸变可通过预畸而得到校正。

以低通数字滤波器为例，将设计步骤归纳如下：

（1）确定数字滤波器的性能指标：通带临界频率fp、阻带临界频率fs；通带内的最大衰减Rp；阻带内的最小衰减As；采样周期T；

=2πfrT；

（2）确定相应的数字角频率，ωp=2πfpT；ω

s

（3）计算经过预畸的相应模拟低通原型的频率，根据Ω和Ωs计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数

Ha(s)；

（4）用上面的双线性变换公式代入Ha(s)，求出所设计的传递函数H(z)；

（5）分析滤波器特性，检查其指标是否满足要求。

五、实验过程

function pr7()

%31 IIR filter

wp=0.2*pi;

ws=0.3*pi;

Rp=1;

As=15;

T=1;

Fs=1/T;

OmegaP=(2/T)*tan(wp/2);

OmegaS=(2/T)*tan(ws/2);

ep=sqrt(10^(Rp/10)-1);

Ripple=sqrt(1/(1+ep.^2));

Attn=1/10^(As/20);

N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS)));

OmegaC=OmegaP/((10.^(Rp/10)-1).^(1/(2*N)));

[cs,ds]=u_buttap(N,OmegaC);

[b,a]=bilinear(cs,ds,Fs);

[mag,db,pha,w]=freqz_m(b,a);

subplot(3,1,1);

plot(w/pi,mag);

title('幅度响应');

xlabel('w(pi)');

ylabel('H');

axis([0,1,0,1.1]);

set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);

set(gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]);

grid;

subplot(3,1,2);

plot(w/pi,db);

title('幅度响应(dB)');

xlabel('w(pi)');

ylabel('H');

axis([0,1,-40,5]);

set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);

set(gca,'YTickmode','manual','YTick',[-50,-15,-1,0]);

grid;

subplot(3,1,3);

plot(w/pi,pha);

title('相位响应');

xlabel('w(pi)');

ylabel('pi unit');

%axis([0,1,0,1.1]);