清风Get西电数字信号处理大作业.docx
《清风Get西电数字信号处理大作业.docx》由会员分享，可在线阅读，更多相关《清风Get西电数字信号处理大作业.docx（26页珍藏版）》请在冰豆网上搜索。
清风Get西电数字信号处理大作业
西电-数字信号处理大作业
数字信号处理上机大作业
实验一：
信号、系统及系统响应
（1）简述实验目的及实验原理。
1.实验目的
●熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。
●熟悉时域离散系统的时域特性。
●利用卷积方法观察分析系统的时域特性。
●掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
2.实验原理与方法
●时域采样。
●LTI系统的输入输出关系。
（2）按实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时域和幅频特性曲线，并对所得结果进行分析和解释。
Matlab源程序如下：
A=1;
T1=1/1000;
T2=1/300;
T3=1/200;
a=25*pi;
w0=30*pi;
n=0:
99;
x1=A*exp（-a*n*T1）.*sin（w0*n*T1）;
x2=A*exp（-a*n*T2）.*sin（w0*n*T2）;
x3=A*exp（-a*n*T3）.*sin（w0*n*T3）;
m=linspace（-pi,pi,10000）;
X1=x1*exp（-j*n'*m）;%n'与m构造矩阵，xi向量与矩阵每一列相乘对应元素相加，构成DTFT后的矩阵
X2=x2*exp（-j*n'*m）;
X3=x3*exp（-j*n'*m）;
figure
（1）;
subplot（3,2,1）
plot（m/pi,abs（X1））;
xlabel（'\omega/π'）;ylabel（'|H（e^j^\omega）|'）;
title（'采样频率为1000Hz时的幅度谱'）;
subplot（3,2,3）
plot（m/pi,abs（X2））;
xlabel（'\omega/π'）;ylabel（'|H（e^j^\omega）|'）;
title（'采样频率为300Hz时的幅度谱'）;
subplot（3,2,5）
plot（m/pi,abs（X3））;
xlabel（'\omega/π'）;ylabel（'|H（e^j^\omega）|'）;
title（'采样频率为200Hz时的幅度谱'）;
subplot（3,2,2）
plot（n,abs（x1））;
xlabel（'n'）;ylabel（'x1（t）'）;
title（'采样频率为1000Hz时的时域波形'）;
subplot（3,2,4）
plot（n,abs（x2））;
xlabel（'n'）;ylabel（'x2（t）'）;
title（'采样频率为300Hz时的时域波形'）;
subplot（3,2,6）
plot（n,abs（x3））;
xlabel（'n'）;ylabel（'x3（t）'）;
title（'采样频率为200Hz时的时域波形'）;
波形图如下：
②时域离散信号、系统和系统响应分析。
Matlab源程序如下:
xb=[1];
xc=ones（1,10）;
ha=ones（1,10）;
hb=[1,2.5,2.5,1];
y=conv（xb,hb）;
n1=0:
length（y）-1;
n2=0:
length（hb）-1;
figure
（1）
subplot（2,1,1）;
stem（n1,y,'filled'）;
xlabel（'n'）;ylabel（'y（n）'）;
title（'y（n）的时域响应'）;
subplot（2,1,2）;
stem（n1,hb,'filled'）;
xlabel（'n'）;ylabel（'hb（n）'）;
title（'hb（n）的时域相应'）;
w=linspace（-pi,pi,10000）;
Y=y*exp（-j*n1'*w）;
Hb=hb*exp（-j*n2'*w）;
figure
（2）
subplot（2,2,1）;
plot（w/pi,abs（Y））;
xlabel（'\omega/π'）;ylabel（'幅度'）;
title（'DTFT[y（n）]的幅度'）;
subplot（2,2,2）;
plot（w/pi,angle（Y））;
xlabel（'\omega/π'）;ylabel（'相位'）;
title（'DTFT[y（n）]的相位'）;
subplot（2,2,3）;
plot（w/pi,abs（Hb））;
xlabel（'\omega/π'）;ylabel（'幅度'）;
title（'DTFT[Hb（n）]的幅度'）;
subplot（2,2,4）;
plot（w/pi,angle（Hb））;
xlabel（'\omega/π'）;ylabel（'相位'）;
title（'DTFT[Hb（n）]的相位'）;
z=conv（xc,ha）;
n3=0:
length（z）-1;
Z=z*exp（-j*n3'*w）;
figure（3）;
subplot（2,1,1）;
plot（w/pi,abs（Z））;
xlabel（'\omega/π'）;ylabel（'幅度'）;
title（'DTFT[z（n）]的幅度'）;
subplot（2,1,2）;
plot（w/pi,angle（Z））;
xlabel（'\omega/π'）;ylabel（'相位'）;
title（'DTFT[z（n）]的相位'）;
波形图如下：
系统ha（n）对信号xc（n）的响应特性：
③由题目②中图2即可验证卷积定理。
（3）实验中的主要结论。
1.时域采样定理
只有当采样频率fs.max大于信号中最高频率fmax的2倍时（fs.max>2fmax），即满足奈奎斯特定律的时候，采样之后的数字信号才能够不发生混叠，保留原有信号，不失真。
2.卷积定理：
时域卷积等于频域相乘。
3.任何函数和单位脉冲函数卷积得到的都是它本身。
4.当所取得N不同时，卷积出来的结果也不同。
（4）思考题
●在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？
它们所对应的模拟频率是否相同？
为什么？
答：
数字频率度量不相同，但他们所对应的模拟频率相同。
由w=Ω*Ts得，采样间隔变化时模拟频率对应的数字频率会有相应的变化，故其度量会有所变化。
●在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？
答：
有差异，所得到的结果点数不同
实验二：
用FFT作谱分析
（1）简述实验目的及实验原理。
实验目的
●进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。
●熟悉FFT算法原理和FFT子程序的应用。
●学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
（2）结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。
1）Matlab源程序如下
N=64;
n=0:
999;
fs=50;
T=1/fs;
x1=ones（1,4）;
x2=[1,2,3,4,4,3,2,1];
x3=[4,3,2,1,1,2,3,4];
x4=cos（0.25*pi*n）;
x5=sin（0.125*pi*n）;
x6=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
X10=fft（x1,N）;
X20=fft（x2,N）;
X30=fft（x3,N）;
X40=fft（x4,N）;
X50=fft（x5,N）;
X60=fft（x6,N）;
X1=fftshift（X10）;
X2=fftshift（X20）;
X3=fftshift（X30）;
X4=fftshift（X40）;
X5=fftshift（X50）;
X6=fftshift（X60）;
k1=0:
31;
k2=32:
63;
w1=2*pi*fs*k1/N;
w2=（2*pi*k2/N-2*pi）*fs;
w=[w2w1];
figure
（1）
subplot（2,1,1）;
n1=0:
length（x1）-1;
stem（n1,x1）;
xlabel（'n'）;ylabel（'x1（n）'）;
title（'x1（n）的时域波形'）;
subplot（2,1,2）;
stem（w,abs（X1））;
xlabel（'w'）;ylabel（'|X1（k）|'）;
title（'DFT（x1）幅频特性'）
figure
（2）
subplot（2,1,1）;
n2=0:
length（x2）-1;
stem（n2,x2）;
xlabel（'n'）;ylabel（'x2（n）'）;
title（'x2（n）的时域波形'）;
subplot（2,1,2）;
stem（w,abs（X2））;
xlabel（'w'）;ylabel（'|X2（k）|'）;
title（'DFT（x2）幅频特性'）
figure（3）
subplot（2,1,1）;
n3=0:
length（x3）-1;
stem（n3,x3）;
xlabel（'n'）;ylabel（'x3（n）'）;
title（'x3（n）的时域波形'）;
subplot（2,1,2）;
stem（w,abs（X3））;
xlabel（'w'）;ylabel（'|X3（k）|'）;
title（'DFT（x3）幅频特性'）
figure（4）
stem（w,abs（X4））;
xlabel（'w'）;ylabel（'|X4（k）|'）;
title（'DFT（x4）幅频特性'）
figure（5）
stem（w,abs（X5））;
xlabel（'w'）;ylabel（'|X5（k）|'）;
title（'DFT（x5）幅频特性'）
figure（6）
stem（w,abs（X6））;
xlabel（'w'）;ylabel（'|X6（k）|'）;
title（'DFT（x6）幅频特性'）
波形图如下：
（2）x（n）=x4（n）+x5（n）
8点离散傅立叶变换Matlab源程序如下：
N=8;
n=0:
7;
x4=cos（0.25*pi*n）;
x5=sin（0.125*pi*n）;
xn=x4+x5;
X10=fft（xn,N）;
X1=fftshift（X10）;
figure
（1）
stem（abs（X1））;axisauto;
xlabel（'k'）;ylabel（'|Xn（k）|'）;
title（'8点DFT[x（n）]幅频特性'）
波形图如下：
16点离散傅立叶变换Matlab源程序如下：
M=16;
n=0:
15;
x4=cos（0.25*pi*n）;
x5=sin（0.125*pi*n）;
xn=x4+x5;
X20=fft（xn,M）;
X2=fftshift（X20）;
stem（abs（X2））;axisauto;
xlabel（'k'）;ylabel（'|Xn（k）|'）;
title（'16点DFT[x（n）]幅频特性'）
波形图如下：
（3）x（n）=x4（n）+jx5（n）
8点离散傅立叶变换Matlab源程序如下：
N=8;
n=0:
7;
x4=cos（0.25*pi*n*T）;
x5=sin（0.125*pi*n*T）;
xn=x4+j*x5;
X10=fft（xn,N）;
X1=fftshift（X10）;
figure
（1）
stem（abs（X1））;
xlabel（'w'）;ylabel（'|Xn（k）|'）;
title（'8点DFT[x（n）]幅频特性'）
波形图如下：
16点离散傅立叶变换Matlab源程序如下：
M=16;
n=0:
15;
x4=cos（0.25*pi*n）;
x5=sin（0.125*pi*n）;
xn=x4+j*x5;
X20=fft（xn,M）;
X2=fftshift（X20）;
figure
（1）
stem（abs（X2））;axisauto;
xlabel（'k'）;ylabel（'|Xn（k）|'）;
title（'16点DFT[x（n）]幅频特性'）
波形图如下：
（3）总结实验所得主要结论。
序列的N点离散傅里叶变换是其z变换在单位圆上的N点等间隔采样，也是其离散时间傅里叶变换在[0,2π]区间上的N点等间隔采样。
DFT的变换区间长度N不同，表示对X（z）的单位圆上，或对其离散时间傅里叶变换在[0,2π]区间上的采样间隔和采样点数不同，因而DFT的变换结果不同。
信号的DFT变换分为共轭对称部分合共个反对称部分，即信号的实部对应离散傅里叶变换的共轭对称部分，虚部对应离散信号的共轭反对称部分。
（4）思考题
●在N=8时，x2（n）和x3（n）的幅频特性会相同吗?
为什么?
N=16呢?
答：
N=8时幅频特性一样，N=16时幅频特性不一样。
●如果周期信号的周期预先不知道，如何用FFT进行谱分析?
答：
设一个定长的m值,先取2m,看2m/m的误差是否大，如大的话再取4m,看4m/2m的误差是否大，如不大，4m（4倍的m值）则可近似原来点的谱分析。
实验三：
用窗函数法设计FIR数字滤波器
（1）简述实验目的及实验原理。
实验目的
（1）掌握用窗函数法设计FIR数字滤波器的原理和方法。
（2）熟悉线性相位FIR数字滤波器特性。
（3）了解各种窗函数对滤波特性的影响。
（2）实验结果
●矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数及频谱特
Matlab源程序如下：
N=31;
wdrec=boxcar（N）;
[hrec,w]=freqz（wdrec,1）;
figure
（1）
subplot（1,2,1）
stem（wdrec）;
xlabel（'n'）;ylabel（'\omega（n）'）;
title（'矩形窗函数'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（hrec）/abs（hrec
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'矩形窗函数幅度特性'）
wdhan=hanning（N）;
[hhan,w]=freqz（wdhan,1）;
figure
（2）
subplot（1,2,1）
stem（wdhan）;
xlabel（'n'）;ylabel（'\omega（n）'）;
title（'汉宁窗函数'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（hhan）/abs（hhan
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'汉宁窗函数幅度特性'）
wdham=hamming（N）;
[hham,w]=freqz（wdham,1）;
figure（3）
subplot（1,2,1）
stem（wdham）;
xlabel（'n'）;ylabel（'\omega（n）'）;
title（'汉明窗函数'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（hham）/abs（hham
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'汉明窗函数幅度特性'）
wdblack=blackman（N）;
[hblack,w]=freqz（wdblack,1）;
figure（4）
subplot（1,2,1）
stem（wdblack）;
xlabel（'n'）;ylabel（'\omega（n）'）;
title（'布莱克曼窗函数'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（hblack）/abs（hblack
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'布莱克曼窗函数幅度特性'）
波形图如下：
●利用不同窗函数与hd（n）构成系统，其单位脉冲相应及频谱特性
Matlab源程序如下：
N=31;
a=（N-1）/2;
fp=120;
fs=150;
Fs=1000;%以上数据自由选取
wc=（fp+fs）/Fs;
h1=fir1（N-1,wc,boxcar（N））;
h2=fir1（N-1,wc,hanning（N））;
h3=fir1（N-1,wc,hamming（N））;
h4=fir1（N-1,wc,blackman（N））;
[H1,w]=freqz（h1,1）;
[H2,w]=freqz（h2,1）;
[H3,w]=freqz（h3,1）;
[H4,w]=freqz（h4,1）;
figure
（1）
subplot（1,2,1）
stem（h1）;
xlabel（'n'）;ylabel（'h1（n）'）;
title（'矩形窗数字滤波器单位脉冲响应'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（H1）/abs（H1
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'矩形窗数字滤波器幅频特性'）
figure
（2）
subplot（1,2,1）
stem（h2）;
xlabel（'n'）;ylabel（'h2（n）'）;
title（'汉宁窗数字滤波器单位脉冲响应'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（H2）/abs（H2
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'汉宁窗数字滤波器幅频特性'）
figure（3）
subplot（1,2,1）
stem（h3）;
xlabel（'n'）;ylabel（'h3（n）'）;
title（'汉明窗数字滤波器单位脉冲响应'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（H3）/abs（H3
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'汉明窗数字滤波器幅频特性'）
figure（4）
subplot（1,2,1）
stem（h4）;
xlabel（'n'）;ylabel（'h4（n）'）;
title（'布莱克曼窗数字滤波器单位脉冲响应'）
subplot（1,2,2）;
plot（w/pi,10*log（abs（H4）/abs（H4
（1））））;
xlabel（'ω/Π'）;ylabel（'W（\omega）/dB'）
title（'布莱克曼窗数字滤波器幅频特性'）
波形图如下：
（3）总结用窗函数法设计FIR滤波器的主要特点。
FIR滤波器目前常用的设计方法有窗函数法和频率采样法,窗函数法是从时域进行设计,而频率采样法是从频域进行设计.窗函数法由于简单、物理意义清晰,因而得到了较为广泛的应用.
窗函数法设计的基本思想是：
首先根据技术指标要求,选取合适的阶数N和窗函数的类型w（n）,使其幅频特性逼近理想滤波器幅频特性.其次,因为理想滤波器的hd（n）是无限长的,所以需要对hd（n）进行截断,数学上称这种方法为窗函数法.简而言之,用窗函数法设计FIR滤波器是在时域进行的,先用傅里叶变换求出理想滤波器单位抽样相应hd（n）,然后加时间窗w（n）对其进行截断,以求得FIR滤波器的单位抽样响应h（n）.
（3）思考题
（1）如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器?
写出设计步骤。
答：
①根据所需设计的数字滤波器类型，确定线性相位数字滤波器类型，
②选择合适的窗函数；
③确定理想低通数字滤波器的频率响应函数；
④计算理想低通数字滤波器的单位脉冲响应；
⑤加窗得到设计结果；
如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd（n）。
答：
hd（n）=1/2π∫w1w0e-jwaejnwdw+1/2π∫w2w3e-jwaejnwdw（其中w0=-w0-wc,w1=-w0+wc,w2=w0-wc,w3=w0+wc）
计算整理后可得：
hd（n）=2/（（n-a）*π）*sin[（n-a）wc]*cos[（n-a）w0]
=2wc/π*sa[（n-a）wc]*cos[（n-a）w0]
升级会员 