数字信号处理课程设计
山东工商学院
程 设 计 报 告
课程名称: 数字信号处理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('傅里叶变换'); %标题傅立叶变换 结果:
模拟信号
1
4
3
0.8
2-45
傅里叶变换
X a (j Ω)
0t/s
0.005
0.01
0.6
x a (n )
10-1
0.4
0.2
-2
-0.01-0.005
-3-40
-20
0w/pi
2040
分析:模拟信号在[-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 结果:
fs=1000Hz的DTFT
3
X 1(j w )
2
10-1
-0.8-0.6-0.4-0.2
00.2w/pi
fs=5000Hz的DTFT
0.40.60.81
10
X 2(j w )
5
0-1
-0.8-0.6-0.4-0.2
0w/pi
0.20.40.60.81
。 和X (k )
分析:
当采样频率越大的时候,采样信号频谱越陡峭,而其失真情况也越来越小。 题目2:
已知时域信号 x (n ) = cos (0.48πn ) + cos (0.52πn ) ,求 下 面 5 种情况的 X (e j ω )
(1)取 x (n ) 的前 10 点数据,求 N = 10 点的 X (e j ω ) 和 X (k ) ,并作图。
j ω ) 和 (k ) (2)将(1)中的 x (n ) 补零至 100点,求 N = 100 点的 X (e X ,并作图。 (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
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]) 结果:
x(n) (0
x (n )
1
2
3
4
5n DTFT
6
7
8
9
10
10
x (j w )
500
1
2
3w DFT
4
5
6
x (k )
1
2
3
4
5k
6
7
8
9
10
分析:
由图可见,由于截断函数的频谱混叠作用,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
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]) 结果:
x(n) (0
x (n )
10
20
30
40
50n DTFT
60
70
80
90
100
10
x (j w )
500
1
2
3w DFT
4
5
6
x (k )
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
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]) 结果:
x(n) (0
x (n )
10
20
30
40
50n DTFT
60
70
80
90
100
50
x (j w )
00
12
3w DFT
456
x (k )
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
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]) 结果:
x(n) (0
x (n )
20
40
60n DTFT
80
100
120
x (j w )
50
012
3w DFT
456
x (k )
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
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
x (n )
5
10
15
20
25n DTFT
30
35
40
45
50
30
x (j w )
201000
1
2
3w DFT
4
5
6
x (k )
k
由图可见,截断函数的宽度正好为序列的周期,即这时的“物理”分辨率使X (K )正好能正确分辨w1=0.48*pi,w2=0.52*pi这两个频率分量。
题目3:
利用MATLAB 编程设计一个数字带通滤波器,指标要求如下: 通带边缘频率:
w p 1=0. 45π
,
w p 2=0. 65π
,通带峰值起伏:
αp ≤1[dB ]
阻带边缘频率:w s 1=0. 3π,w s 2=0. 75π最小阻带衰减:αs ≥40[dB ]
分别用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 幅度和相位频响曲线:
IIR 幅度频响曲线
1.5
1
幅度0.5
000.10.20.30.50.6
w/pi
IIR 相位频响曲线0.40.70.80.91
4
2
相位0
-2
-400.10.20.30.40.5
w/pi0.60.70.80.91
(2)
FIR 设计:
由于αs ≥40[dB ],查表可选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 幅度和相位频响曲线:
FIR 幅度频响曲线
1.5
1
幅度0.5
000.10.20.30.50.6
w/pi
FIR 相位频响曲线0.40.70.80.91
4
2
相位0
-2
-400.10.20.30.40.5
w/pi0.60.70.80.91
分析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分钟就调试结束了,期间还给我们讲解了过程,删了一些代码,又加了一些代码,运行后,非常成功,然后我们自己对老师的这种方法理解了大部分,剩下的一些题,跟老师讲的这个题是同类型的,只是参数不同,然后我们又开始调试了,这一次比较快,一个多小时,初步完成了,然后又接着写代码注释,截图,分析,在宿舍写了一晚上,最后终于基本完成了本次课程设计,我相信这次经历对以后的学习与工作都有很大的帮助,另外,衷心的感谢老师的指导,如果没有老师的指导,课设不会这么快的完成,老师,您辛苦了。