数字信号处理课程设计

山东工商学院

课程设计报告

课程名称：数字信号处理A

班级：XXXXX

姓名: XXXX

学号：XXXXX

指导教师：XXXX

时间：2016年1月1日

一、课程设计题目

题目1： （1）、已知Xa(t)=e^-1000|t|,求其傅立叶变换Xa(j Ω) ，画出模拟信号及其连续时间傅里叶变换的曲线图。 （2）、以Xa(t)为例，说明采样频率对频率响应的影响，分别采用fs=1000Hz 和fs=5000Hz ，绘出X(e^jw)曲线。 （1） 代码: close all clear;clc;

W=10;f=1000;n=-10:W-1;t=n/f; X=exp(-1000*abs(t));

subplot(1,2,1);plot(t,X); %画模拟信号曲线 xlabel('t/s');ylabel('xa(n)');

title('模拟信号'); %标题模拟信号 tf=10;N=100;dt=10/N;t=(1:N)*dt; wf=25;Nf=50;

w1=linspace(0,wf,Nf); %0-25之间分成50点 dw=wf/(Nf-1); W1=-50:50;

Xat=exp(-1000*abs(t)); %表达式 F1=Xat*exp(-1i*t'*w1)*dt; %傅立叶变换 w=[-fliplr(w1),w1(2:Nf)]; %负频率的频谱

Y1=(exp(2)-1)./(exp(2)-exp(1-1i*W1)-exp(1+1i*W1)+1); F=[fliplr(F1),F1(2:Nf)];t=[-fliplr(t),t]; subplot(1,2,2);

plot(w,F,'linewidth',1); %画傅立叶变换曲线 xlabel('w/pi');ylabel('Xa(j Ω)');

title('傅里叶变换'); %标题傅立叶变换 结果：

-0.01-0.005

00.005

0.01

0.2

0.4

0.6

0.8

1

t/s

x a (n )

模拟信号

-40

-20

02040

-3-2

-1

0123

4

-45

w/pi

X a (j Ω)

傅里叶变换

分析：模拟信号在[-0.01，0.01]区间为连续信号,其傅立叶变换曲线在[-10，10]内为连续曲线。 （2） 代码： close all clear clc

Dt=0.00005; %步长为0.00005s t=-0.005:Dt:0.005;

xa=exp(-1000*abs(t)); %取时间从-0.005s 到0.005s 这段模拟信号 Ts1=0.001;Ts2=0.0002; %周期 n=-25:1:25;

x1=exp(-1000*abs(n*Ts1)); x2=exp(-1000*abs(n*Ts2));

K=100;k=0:1:K;w=pi*k/K; %求模拟角频率 X1=x1*exp(-j*n'*w); %求其傅立叶变换 X2=x2*exp(-j*n'*w); %求其傅立叶变换 X11=real(X1);X12=real(X2);

w=[-fliplr(w),w(2:101)]; %将角频率范围扩展为从-到+ X11=[fliplr(X11),X11(2:101)]; X12=[fliplr(X12),X12(2:101)];

subplot(2,1,1);plot(w/pi,X11);%画出fs=1000Hz 的频率响应 xlabel('w/pi');ylabel('X1(jw)');

title('fs=1000Hz 的DTFT'); %标题fs=1000Hz 的DTFT subplot(2,1,2);plot(w/pi,X12);%画出fs=5000Hz 的频率响应 xlabel('w/pi');ylabel('X2(jw)');

title('fs=5000Hz 的DTFT'); %标题fs=5000Hz 的DTFT 结果：

0123

w/pi

X 1(j w )

fs=1000Hz 的

DTFT

-1

-0.8-0.6-0.4-0.2

00.20.40.60.81

05

10

w/pi

X 2(j w )

fs=5000Hz 的DTFT

分析：

当采样频率越大的时候，采样信号频谱越陡峭，而其失真情况也越来越小。题目2：

已知时域信号x(n)= cos(0.48πn)+cos(0.52πn)，求下面 5 种情况的X (e jω)和X (k )。

（1）取x(n)的前10 点数据，求N=10点的X(e jω)和X(k )，并作图。

（2）将（1）中的x(n)补零至100点，求N=100 点的X(e jω)和X(k )，并作图。（3）取x(n)的前100 点数据，求N=100 点的X(e jω)和X(k )，并作图。（4）取x(n)的前128 点数据，求N=128 点的X(e jω)和X(k )，并作图。

（5）取x(n)的前50 点数据，求N= 50 点的X(e jω)和X(k )，并作图。

讨论以上5种情况的区别。

(1)

代码：

close all

clear

clc

n=(0:1:9);

y=cos(0.48*pi*n)+cos(0.52*pi*n);

n1=(0:1:9);x=y(1:1:10);

subplot(3,1,1);

stem(n1,x); %画出x(n)曲线

title('x(n) (0<=n<=9)'); %标题0<=n<=9)

xlabel('n');ylabel('x(n)');

axis([0,10,-2.5,2.5]); %axis([xmin xmax ymin ymax])

w=linspace(0,2*pi,length(x)); %0-2*pi区域分为10点

xw=x*exp(-j*[1:length(x)]'*w);

magx=abs(xw); %对xw取绝对值

subplot(3,1,2);

plot(w,magx); %画出x(jw)曲线

title('DTFT'); %标题DTFT

xlabel('w');ylabel('x(jw)');

axis([0,2*pi,0,10]); %axis([xmin xmax ymin ymax])

subplot(3,1,3);

x1=fft(x); %对x进行傅立叶变换

magx1=abs(x1); %对x1取绝对值

stem(n1,abs(magx1)); %画出x(k)曲线

title('DFT'); %标题DFT

axis([0,10,0,10]); %axis([xmin xmax ymin ymax]) 结果：

1

2

3

4

56

7

8

9

10

x(n) (0<=n<=9)

n x (n )

0510

DTFT

w x (j w )

DFT

k

x (k )

分析：

由图可见，由于截断函数的频谱混叠作用，X （K ）不能正确分辨w1=0.48*pi,w2=0.52*pi 这两个频率分量。 (2) 代码： close all clear clc

n=(0:1:9);

y=cos(0.48*pi*n)+cos(0.52*pi*n); n1=(0:1:99);

x=[y(1:1:10),zeros(1,90)]; %第10位到100位的数据都为0 subplot(3,1,1);

stem(n1,x); %画出x(n)曲线

title('x(n) (0<=n<=9+90zeros)'); %标题0<=n<=9+90zeros) xlabel('n');ylabel('x(n)');

axis([0,100,-2.5,2.5]); %axis([xmin xmax ymin ymax]) w=linspace(0,2*pi,length(x)); %0-2*pi 区域分为100点 xw=x*exp(-j*[1:length(x)]'*w);

magx=abs(xw); %对xw 取绝对值 subplot(3,1,2);

plot(w,magx); %画出x(jw)曲线 title('DTFT'); %标题DTFT

axis([0,2*pi,0,10]); %axis([xmin xmax ymin ymax]) subplot(3,1,3);

x1=fft(x); %对x 进行傅立叶变换 magx1=abs(x1); %对x1取绝对值 stem(n1,abs(magx1)); %画出x(k)曲线 title('DFT'); %标题DFT xlabel('k');ylabel('x(k)');

axis([0,100,0,10]); %axis([xmin xmax ymin ymax]) 结果：

10

20

30

40

5060

70

80

90

100

x(n) (0<=n<=9+90zeros)

n x (n )

0510

DTFT

w x (j w )

DFT

k

x (k )

分析：

由图可见，虽然x(n)补零至100点，X （K ）的密度，截断函数的频谱混叠作用没有改变，这时的物理分辨率使X(K)仍不能正确分辨w1=0.48*pi,w2=0.52*pi 这两个频率分量。 (3) 代码： close all clear clc

n=(0:1:99);

y=cos(0.48*pi*n)+cos(0.52*pi*n); n1=(0:1:99);x=y(1:1:100); subplot(3,1,1);

stem(n1,x); %画出x(n)曲线 title('x(n) (0<=n<=99)'); %标题0<=n<=99) xlabel('n');ylabel('x(n)');

axis([0,100,-2.5,2.5]); %axis([xmin xmax ymin ymax]) w=linspace(0,2*pi,length(x)); %0-2*pi 区域分为100点 xw=x*exp(-j*[1:length(x)]'*w);

magx=abs(xw); %对xw 取绝对值 subplot(3,1,2);

plot(w,magx); %画出x(jw)曲线 title('DTFT'); %标题DTFT xlabel('w');ylabel('x(jw)');

axis([0,2*pi,0,54]); %axis([xmin xmax ymin ymax]) subplot(3,1,3);

x1=fft(x); %对x 进行傅立叶变换 magx1=abs(x1); %对x1取绝对值 stem(n1,abs(magx1)); %画出x(k)曲线 title('DFT'); %标题DFT xlabel('k');ylabel('x(k)');

axis([0,100,0,54]); %axis([xmin xmax ymin ymax]) 结果：

10

20

30

40

5060

70

80

90

100

x(n) (0<=n<=99)

n x (n )

12

3456

050

DTFT

w x (j w )

DFT

k

x (k )

分析：

由图可见，截断函数的加宽且为周期序列的整数倍，改变了频谱混叠作用，提高了“物理”分辨率使X （K ）能正确分辨w1=0.48*pi,w2=0.52*pi 这两个频率分量。 (4) 代码： close all clear clc

n=(0:1:127);

y=cos(0.48*pi*n)+cos(0.52*pi*n); n1=(0:1:127);x=y(1:1:128); subplot(3,1,1);

stem(n1,x); %画出x(n)曲线 title('x(n) (0<=n<=127)'); %标题0<=n<=127) xlabel('n');ylabel('x(n)');

axis([0,128,-2.5,2.5]); %axis([xmin xmax ymin ymax]) w=linspace(0,2*pi,length(x)); %0-2*pi 区域分为128点 xw=x*exp(-j*[1:length(x)]'*w);

magx=abs(xw); %对xw 取绝对值 subplot(3,1,2);

plot(w,magx); %画出x(jw)曲线 title('DTFT'); %标题DTFT xlabel('w');ylabel('x(jw)');

axis([0,2*pi,0,65]); %axis([xmin xmax ymin ymax]) subplot(3,1,3);

x1=fft(x); %对x 进行傅立叶变换 magx1=abs(x1); %对x1取绝对值 stem(n1,abs(magx1)); %画出x(k)曲线 title('DFT'); %标题DFT xlabel('k');ylabel('x(k)');

axis([0,128,0,65]); %axis([xmin xmax ymin ymax]) 结果：

20

40

6080

100

120

x(n) (0<=n<=127）

n x (n )

50

DTFT

w x (j w )

DFT

k

x (k )

分析：

由图可见，截断函数虽进一步加宽，但不是周期序列的整数倍，所以尽管 X （K ）能正确分辨w1=0.48*pi,w2=0.52*pi 这两个频率分量，但还是呈现除了频

(5) 代码： close all clear clc n=(0:1:49);

y=cos(0.48*pi*n)+cos(0.52*pi*n); n1=(0:1:49); x=y(1:1:50);

subplot(3,1,1);stem(n1,x); %画出x(n)曲线 title('x(n) (0<=n<=49)'); %标题0<=n<=49) xlabel('n');ylabel('x(n)');

axis([0,50,-2.5,2.5]); %axis([xmin xmax ymin ymax]) w=linspace(0,2*pi,length(x)); %0-2*pi 区域分为50点 xw=x*exp(-j*[1:length(x)]'*w);

magx=abs(xw); %对xw 取绝对值 subplot(3,1,2);plot(w,magx); %画出x(jw)曲线 title('DTFT'); %标题DTFT xlabel('w');ylabel('x(jw)');

axis([0,2*pi,0,30]); %axis([xmin xmax ymin ymax]) subplot(3,1,3);

x1=fft(x); %对x 进行傅立叶变换 magx1=abs(x1); %对x1取绝对值 stem(n1,abs(magx1)); %画出x(k)曲线 title('DFT'); %标题DFT xlabel('k');ylabel('x(k)');

axis([0,50,0,30]); %axis([xmin xmax ymin ymax]) 结果：

x(n) (0<=n<=49）

n x (n )

0102030

DTFT

w x (j w )

DFT

k

x (k )

由图可见，截断函数的宽度正好为序列的周期，即这时的“物理”分辨率使X（K）正好能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量。

题目3：

利用MATLAB编程设计一个数字带通滤波器，指标要求如下：

通带边缘频率：

π

45

.0

1

=

p

w

，

π

65

.0

2

=

p

w

，通带峰值起伏：

]

[1dB

p

≤

α

阻带边缘频率：

π3.0

1

=

s

w，π

75

.0

2

=

s

w最小阻带衰减：]

[

40dB

s

≥

α

分别用IIR和FIR两种数字滤波器类型进行设计。

实验要求：给出IIR数字滤波器参数和FIR数字滤波器的冲激响应，绘出它们的

幅度和相位频响曲线，讨论它们各自的实现形式和特点。

（1）

代码：

close all

clear

clc

ws=[0.3,0.75];

wp=[0.45,0.65];

Rp=1;Rs=40;

[N,wc]=buttord(wp,ws,Rp,Rs); %求阶数N及频率参数wc

[B,A]=butter(N,wc); %设计巴特沃斯滤带通波器

[H,w]=freqz(B,A); %求频响

subplot(2,1,1);

plot(w/pi,abs(H)); %画幅频曲线

title('IIR幅度频响曲线'); %标题IIR幅度频响曲线

xlabel('w/pi');

ylabel('幅度');

subplot(2,1,2);

plot(w/pi,angle(H)); %画相频曲线

title('IIR相位频响曲线'); %标题IIR相位频响曲线

xlabel('w/pi');

ylabel('相位');

IIR数字滤波器参数：

N = 7

wc = 0.4284 0.6693

B = 0.0003 0 -0.0019 0 0.0057 0 -0.0095

0 0.0095 0 -0.0057 0 0.0019 0 -0.0003

A =1.0000 1.7451 4.9282 6.1195 9.8134 9.2245 10.4323

7.5154 6.4091 3.4595 2.2601 0.8470 0.4167 0.0856

0.0299

IIR幅度和相位频响曲线：

0.51

1.5IIR 幅度频响曲线

w/pi

幅度

0.1

0.2

0.3

0.4

0.50.6

0.7

0.8

0.9

1

-4-202

4IIR 相位频响曲线

w/pi

相位

(2)

FIR 设计：

由于

]

[40dB s ≥α，查表可选hamming 窗，其阻带最小衰减-53dB 满足要求。

代码：

close all clear clc

ws1=0.3*pi;ws2=0.75*pi; wp1=0.45*pi;wp2=0.65*pi;

wb=(wp1-ws1+ws2-wp2)/2; %计算过渡带宽

wc=[wp1/pi-(wp1-ws1)/(2*pi); %设置理想带通截止频率 wp2/pi+(ws2-wp2)/(2*pi)];

N=ceil(3.3*2*pi/wb); %计算窗口长度 b=fir1(N,wc); %设计滤波器系数 n=0:N;

hn=b(n+1); %求冲激响应 [H,w]=freqz(b,1); %求频率响应 subplot(2,1,1);

plot(w/pi,abs(H)); %画幅频曲线

title('FIR 幅度频响曲线'); %标题FIR 幅度频响曲线 xlabel('w/pi');ylabel('幅度'); subplot(2,1,2);

plot(w/pi,angle(H)); %画相频曲线

title('FIR 相位频响曲线'); %标题FIR 相位频响曲线 xlabel('w/pi');ylabel('相位'); 冲激响应hn =b(n+1)：

0.0011 0.0005 0.0001 0.0006 -0.0031 -0.0008 0.0051 -0.0002

-0.0002 -0.0016 -0.0109 0.0087 0.0149 -0.0098 -0.0018 -0.0119

-0.0145 0.0463 0.0118 -0.0445 0.0003 -0.0340 0.0250 0.1591

-0.1126 -0.2397 0.2135 0.2135 -0.2397 -0.1126 0.1591 0.0250

-0.0340 0.0003 -0.0445 0.0118 0.0463 -0.0145 -0.0119 -0.0018

-0.0098 0.0149 0.0087 -0.0109 -0.0016 -0.0002 -0.0002 0.0051

-0.0008 -0.0031 0.0006 0.0001 0.0005 0.0011 FIR 幅度和相位频响曲线：

00.10.20.3

0.40.50.60.70.80.91

0.51

1.5FIR 幅度频响曲线

w/pi

幅度

0.1

0.2

0.3

0.4

0.50.6

0.7

0.8

0.9

1

-4-202

4FIR 相位频响曲线

w/pi

相位

分析IIR 与FIR 的实现形式与特点： 1、IIR 滤波器阶数比FIR 少。

IIR 滤波器存在着输出对输入的反馈，因此可以用比FIR 滤波器少的阶数来满足技术指标。

2、FIR 滤波器可得到严格的线性相位，而IIR 滤波器则做不到这一点。 IIR 滤波器的选频特性越好，则相位的非线性就越严重，在需要严格线性相位

的情况下应该选择FIR滤波器。

3、IIR滤波器必须采用递归结构实现，FIR滤波器主要采用非递归结构。

IIR滤波器必须采用递归结构实现，只有当所有极点都在单位圆内时滤波器才是稳定的。但实际中由于存在有限字长效应，滤波器有可能变得不稳定。而FIR 滤波器主要采用非递归结构，因而从理论上以及从实际的有限精度的运算中，都是稳定的。

4、IIR滤波器有现成的设计公式、数据和表格，FIR滤波器没有现成的设计公式。

IIR滤波器可利用模拟滤波器现成的设计公式、数据和表格，因而计算工作量较小，对计算工具要求不高。FIR滤波器没有现成的设计公式，窗函数法只给出窗函数的计算公式，但计算通带和阻带衰减仍无显式表达式。

5、IIR滤波器主要是设计规格化的、频率特性为分段常数的标准低通、高通、带通、带阻和全通滤波器，而FIR滤波器可设计出理想正交变换器、理想微分器、线性调频器等各种网络，适应性较广。

6、FIR滤波器可以采用快速傅立叶变换（FFT）来实现，在相同阶数下，运算速度可以快得多

三、自我评价和总结

在这将近三天的课设里，我学习了一些与MATLAB相关的基础知识，增强的自己的知识，这段日子，时间过的真的好快，解决一个问题要花费很长时间，这是因为自己的能力缘故，但是每当解决一个问题之后，很开心，真的很开心，收获得很大，刚开始的时候，没有头绪，不知道从何做起，后来自己翻课本，在网上查阅相关的资料，再加上老师的讲解，慢慢的有了头绪，然后自己还挺高兴的，就这样到了最后一天下午，老师看到截图后说周期有问题，我们的缺失了一个周期，必须要调试，于是我们走上了漫漫的调试之旅，很长时间后，问题都还没有解决，然后老师就开始慢慢的指导我们了，真的很佩服老师，感觉老师真的很厉害，不到5分钟就调试结束了，期间还给我们讲解了过程，删了一些代码，又加了一些代码，运行后，非常成功，然后我们自己对老师的这种方法理解了大部分，剩下的一些题，跟老师讲的这个题是同类型的，只是参数不同，然后我们又开始调试了，这一次比较快，一个多小时，初步完成了，然后又接着写代码注释，截图，分析，在宿舍写了一晚上，最后终于基本完成了本次课程设计，我相信这次经历对以后的学习与工作都有很大的帮助，另外，衷心的感谢老师的指导，如果没有老师的指导，课设不会这么快的完成，老师，您辛苦了。