南邮DSPA实验报告

南京邮电大学

实验报告

实验名称：1、离散时间信号与系统的时、频域表示

2、离散傅立叶变换和z变换

3、数字滤波器的频域分析和实现

4、数字滤波器的设计

课程名称数字信号处理A(双语)

班级学号

姓名

指导老师

开课时间 2016/2017学年第一学期

实验一：离散时间信号与系统的时、频域表示

一、实验目的和任务：

●主要内容：序列的生成、序列DTFT谱计算、简单滤波器的仿真，对应实验指导书

Page1－34。

●实验要求：

●（1）按照要求产生基本序列和复杂序列，对序列进行运算，求序列的DTFT

幅度谱和相位谱。

●（2）计算简单的离散时间系统的输出，包括简单滤波器的仿真、计算单个

系统和互连系统的输出。

二、实验内容：

●具体包括

●基本序列产生和运算：Q1.1－Q1.3，Q1.23，Q1.30－Q1.33

●离散时间系统仿真：Q2.1－Q2.3

●LTI系统：Q2.19、Q2.21，Q2.28

●DTFT：Q3.1，Q3.2，Q3.4

三、实验过程与结果分析：

Q1.1 运行P1_1产生单位样本序列u[n] 的程序与显示的波形如下:

clf;

n= -10:20;

u = [zeros(1,10) 1 zeros(1,20)];

stem(n,u);

xlabel('时间序列n');ylabel('振幅');

title('单位样本序列');

axis([-10 20 0 1.2]);

结果：

时间序列 n

振幅

Q1.2

clf 清除波形 axis 设置坐标轴范围，可读比例等 title 给图形加标题 xlabel 给x 加标注 Ylabel 给y 加标注

Q1.3

clf;

n= -10:20;

u = [zeros(1,10) 1 zeros(1,20)]; stem(n+11,u);

xlabel('时间序列 n');ylabel('振幅'); title('单位样本序列'); axis([0 32 0 1.2]);

结果：

0.2

0.4

0.6

0.8

1

时间序列 n

振幅

Q1.23

n = 0:50; f= 0.08; phase = pi/2; A = 2.5;

arg =2*pi*f*n - phase; x = A*cos(arg); clf; stem(n,x); axis([0 50 -3 3]); grid;

title('正弦序列'); xlabel('时间序号n'); ylabel('振幅'); axis;

结果：

051015

20253035404550

正弦序列

时间序号n

振幅

Q1.30

s[n]是线性增加伴随着实指数缓慢衰减的图像加性噪声d[n]是均匀分布在-0.4和+0.4之间的随机序列

Q1.31

不能，因为d 是列向量，s 是行向量

Q1.32

x1是x 的延时，x2和x 相等，x3超前于x

Q1.33

产生图例说明

Q2.1 clf;

n = 0:100;

s1 = cos(2*pi*0.05*n); s2 = cos(2*pi*0.47*n); x = s1+s2;

M = input('滤波器所需的长度='); num = ones(1,M); y = filter(num,1,x)/M; subplot(2,2,1); plot(n,s1); axis([0,100,-2,2]); xlabel('时间序列n');

ylabel('振幅'); title('信号#1'); subplot(2,2,2); plot(n,s2); axis([0,100,-2,2]); xlabel('时间序列n'); ylabel('振幅'); title('信号#2'); subplot(2,2,3); plot(n,x); axis([0,100,-2,2]); xlabel('时间序列n'); ylabel('振幅'); title('输入信号'); subplot(2,2,4); plot(n,y); axis([0,100,-2,2]); xlabel('时间序列n'); ylabel('振幅'); title('输出信号'); axis;

结果：

50

100

-2-101

2时间序列n 振幅

信号#1

50

100

-2-101

2时间序列n 振幅

信号#2

50

100

-2-101

2时间序列n

振幅

输入信号

50

100

-2-101

2时间序列n

振幅

输出信号

S[n]被离散时间系统抑制

Q2.2

n = 0:100;

s1 = cos(2*pi*0.05*n);

s2 = cos(2*pi*0.47*n);

x = s1+s2;

M = input('滤波器所需长度 = ');

num = (-1).^[0:M-1];

y = filter(num,1,x)/M;

clf;

subplot(2,2,1);

plot(n, s1);

axis([0, 100, -2, 2]);

xlabel('时间序号n'); ylabel('振幅'); title('信号 #1');

subplot(2,2,2);

plot(n, s2);

axis([0, 100, -2, 2]);

xlabel(' 时间序号n'); ylabel('振幅'); title('信号 #2');

subplot(2,2,3);

plot(n, x);

axis([0, 100, -2, 2]);

xlabel('时间序号 n'); ylabel('振幅'); title('输入信号');

subplot(2,2,4);

plot(n, y);

axis([0, 100, -2, 2]);

xlabel('时间序号 n'); ylabel('振幅'); title('输出信号');

axis;

显示的波形如下:

50

100

时间序号n 振幅

信号 #1

50

100

时间序号n 振幅

信号 #2

50

100

时间序号 n

振幅

输入信号

50

100

时间序号 n

振幅

输出信号

改变LTI 系统对输入的影响是,系统现在是一个高通滤波器。它通过高频输入组件s2来替代低频输入组件s1.

Q2.3 clf;

n = 0:100;

s1 = cos(2*pi*0.04*n); s2 = cos(2*pi*0.6*n); x = s1+s2;

M = input('滤波器所需的长度='); num = ones(1,M); y = filter(num,1,x)/M; subplot(2,2,1); plot(n,s1); axis([0,100,-2,2]); xlabel('时间序列n'); ylabel('振幅'); title('信号#1'); subplot(2,2,2); plot(n,s2); axis([0,100,-2,2]); xlabel('时间序列n'); ylabel('振幅'); title('信号#2'); subplot(2,2,3); plot(n,x); axis([0,100,-2,2]); xlabel('时间序列n');

ylabel('振幅'); title('输入信号'); subplot(2,2,4); plot(n,y); axis([0,100,-2,2]); xlabel('时间序列n'); ylabel('振幅'); title('输出信号'); axis;

结果：

50

100

-2-101

2时间序列n 振幅

信号#1

50

100

-2-101

2时间序列n 振幅

信号#2

50

100

-2-101

2时间序列n

振幅

输入信号

50

100

-2-101

2时间序列n

振幅

输出信号

Q2.19

clf; N=40;

num=[2.2403 2.4908 2.2403]; den=[1 -0.4 0.75]; y=impz(num,den,N); stem(y);

xlabel('时间序号n'); ylabel('振幅'); title('冲激响应'); grid;

51015

2025303540

-3-2-101

2

3

4

时间序号n

振幅

冲激响应

Q2.21 clf; N=40;

num=[0.9 -0.45 0.35 0.002]; den=[1.0 0.71 -0.46 -0.62]; x=[1 zeros(1,N-1)]; y=filter(num,den,x); stem(y);

xlabel('时间序号n'); ylabel('振幅'); title('冲激响应'); grid;

51015

2025303540

-1.5

-1-0.500.5

1

1.5

2

时间序号n

振幅

冲激响应

Q2.28 clf;

h=[3 2 1 -2 1 0 -4 0 3]; x=[1 -2 3 -4 3 2 1]; y=conv(h,x); n=0:14; subplot(2,1,1); stem(n,y);

xlabel('时间序号n'); ylabel('振幅');

title('用卷积得到的输出'); grid;

x1=[x zeros(1,8)]; y1=filter(h,1,x1); subplot(2,1,2); stem(n,y1); xlabel('时间序号n') ylabel('振幅');

title('由滤波生成的输出'); grid;

时间序号n 振幅

用卷积得到的输出

时间序号n

振幅

由滤波生成的输出

Q3.1

原始序列是：

Pause:暂停命令

Q3.2 clf;

w= -4*pi:8*pi/511:4*pi; num=[2 1]; den=[1 -0.6]; h=freqz(num,den,w); subplot(2,1,1) plot(w/pi,real(h)); grid

title('H(e^{j\omega}的实部') xlabel('\omega/\pi'); ylabel('振幅'); subplot(2,1,2) plot(w/pi,imag(h)); grid

title('H(e^{j\omega}的虚部') xlabel('\omega/\pi'); ylabel('振幅'); pause subplot(2,1,1) plot(w/pi,abs(h)); grid

title('|H(e^{j\omega}|幅度谱') xlabel('\omega/\pi'); ylabel('振幅'); subplot(2,1,2) plot(w/pi,angle(h)); grid

title('相位谱arg[H(e^{j\omega}]') xlabel('\omega/\pi'); ylabel('以弧度为单位的相位'); 结果：

-4

-3

-2

-1

01

2

3

4

0246

8H(e j ω的实部

ω/π振幅

-4

-3

-2

-1

01

2

3

4

-4-202

4H(e j ω的虚部

ω/π

振幅

-4

-3

-2

-1

01

2

3

4

0246

8|H(e j ω|幅度谱

ω/π

振幅

-4

-3

-2

-1

01

2

3

4

-2-101

2相位谱arg[H(e j ω]

ω/π

以弧度为单位的相位

是w 的周期。周期2 π 实部是2π为周期偶对称的 虚部是2π为周期奇对称的 幅度是2π为周期偶对称的 相位是2π为周期奇对称的 Q3.4 clf;

w= -4*pi:8*pi/511:4*pi;

num=[1 3 5 7 9 11 13 15 17]; den=1;

h=freqz(num,den,w); subplot(2,1,1)

plot(w/pi,real(h)); grid

title('H(e^{j\omega}的实部') xlabel('\omega/\pi'); ylabel('振幅'); subplot(2,1,2)

plot(w/pi,imag(h)); grid

title('H(e^{j\omega}的虚部') xlabel('\omega/\pi'); ylabel('振幅');

pause subplot(2,1,1) plot(w/pi,abs(h)); grid

title('|H(e^{j\omega}|幅度谱') xlabel('\omega/\pi'); ylabel('振幅'); subplot(2,1,2) plot(w/pi,angle(h)); grid

title('相位谱arg[H(e^{j\omega}]') xlabel('\omega/\pi'); ylabel('以弧度为单位的相位'); 结果：

-4

-3-2-1

01234

-50050

100H(e j ω的实部

ω/π振幅

-4

-3

-2

-1

01

2

3

4

-100

-50050

100H(e j ω的虚部

ω/π

振幅

-4

-3-2-1

01234

050

100

|H(e j ω|幅度谱

ω/π

振幅

-4

-3

-2

-1

01

2

3

4

-4-202

4相位谱arg[H(e j ω]

ω/π

以弧度为单位的相位

实验二：离散傅里叶变换和z变换

一、实验目的和任务：

●实验主要内容和要求（对应实验指导书3.5和3.6节）：

●（1）DFT计算、常见性质的验证、用DFT实现快速卷积（应用于简单的

含噪信号滤波）；

●（2）实现z变换和逆z变换，（用于：根据差分方程求系统的传递函数，进

一步求冲激响应、零极点图）。

二、实验内容：

●具体包括——

●DFT：Q3.23

●DFT性质：

●圆周移位子函数（Q3.26、Q3.27）

●圆周卷积子函数（Q3.28、Q3.29 ）

●圆周移位（Q3.30－35）

●圆周卷积（Q3.36－40）

●Z变换：Q3.46、Q3.47、Q3.48

●逆Z变换：Q3.50、Q3.51

三、实验过程与结果分析：

Q3.23

clf;

N=200;

L=256;

nn=[0:N-1];

kk=[0:L-1];

xR=[0.1*(1:100) zeros(1,N-100)];

xl=[zeros(1,N)];

x=xR+i*xl;

XF=fft(x,L);

subplot(3,2,1);

grid;

plot(nn,xR);

grid;

title('实\{x[n]\}');

xlabel('时间序号n');

ylabel('振幅');

subplot(3,2,2);

plot(nn,xl);

grid;

title('虚\{x[n]\}');

xlabel('时间序号n');

ylabel('振幅');

subplot(3,2,3);

plot(kk,real(XF));

grid;

title('实\{x[k]\}');

xlabel('频率指数k');

ylabel('振幅');

subplot(3,2,4);

plot(kk,imag(XF));

grid;

title('虚\{x[k]\}');

xlabel('频率指数k');

ylabel('振幅');

xx=ifft(XF,L);

subplot(3,2,5);

plot(kk,real(xx));

grid;

title('IDFT\{x[k]\}实部'); xlabel('时间序号n');

ylabel('振幅');

subplot(3,2,6);

plot(kk,imag(xx));

grid;

title('IDFT\{x[k]\}虚部'); xlabel('时间序号n');

ylabel('振幅');

结果：

50

100150200

05

10

实{x[n]}

时间序号n 振幅

50

100150200

-10

1

虚{x[n]}

时间序号n 振幅

100200300

-1000

1000实{x[k]}

频率指数k 振幅

100200300

-500

500

虚{x[k]}

频率指数k 振幅

100200300

-100

10IDFT{x[k]}实部

时间序号n

振幅

100200300

-10

1

IDFT{x[k]}虚部

时间序号n

振幅

Q3.26在函数circshift 中，命令rem 的作用是什么？

R = rem(X,Y)求余数函数

Q3.27解释函数circshift 怎样实现圆周移位运算。

输入序列x 是循环左移M 位。如果M > 0,那么circshift 删除左边的元素向量x ，并且附加他们到剩下的元素右边来获得循环转移序列。如果如果M < 0,然后circshift 首先补充的x 的长度,最右边的长度(x)- m 样品从x 中移走并且附加在剩下的M 样本右边来得到循环转移序列。

Q3.28在函数circshift 中，运算符~=的作用是什么？

如果A 和B 不相等返回值1 如果A 和B 相等返回值0

Q3.29解释函数circonv 怎样实现圆周卷积运算。

函数circonv 操作如下:输入的是两个相等长度为L 的两个向量x1和x2.，为了理解circonv 是如何工作的,从x2的周期延拓角度来考虑很有用。让x2p 作为x2的无限长的周期延拓。从概念上讲,常规时间反转x2p 并且让x2tr 通过x2p 的时间反转等于元素1。输出向量y 元素1到L 是通过x1和一个长度L 的通过循环右移一个时间反转序列x2tr 得到的序列sh 之间的内积来获得的。对于输出样例y[n],1≤n ≤L 、正确的循环移位是n - 1点。

Q3.30 clf; M=6;

a=[0 1 2 3 4 5 6 7 8 9]; b=circshift(a,M); L=length(a)-1;

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

axis([0,L,min(a),max(a)]); title('原序列'); xlabel('时间序号n'); ylabel('a[n]'); subplot(2,1,2); stem(n,b);

axis([0,L,min(a),max(a)]);

title('通过循环移位得到的序列',num2str(M),'样本'); xlabel('时间序号n'); ylabel('b[n]');

决定时移的数量的部分是M

如果时移的数量大于序列长度，实际实现的循环时移是rem(M,length(a))点左移，相当于循环移动的M 点(不止一次),也相当于通过M 点周期延拓的左移。

Q3.31 结果：

2

468

原序列的DFT 的幅度

频率序号k

|X [k ]|

2

468

圆周位移后的DFT 的幅度

频率序号k

|Y [k ]|

2

468

原序列的DFT 的幅度频率序号k

a r g （X [k ]）

2

468

圆周位移后的DFT 相位

频率序号k

a r g （Y [k ]）

序列的长度是10，并且M = 12。这可能被解释为一个12点的循环左移(不止一次),作为一个2点循环左移,或者作为一个2的线性左移，或者序列的12点周期延拓。

Q3.32& Q3.33.

程序：