数字信号处理上机作业
学院： 电子工程学院
班级：021215
组员：
实验一：信号、系统及系统响应
1、实验目的
(1) 熟悉连续信号经理想采样前后的频谱变化关系， 加深对时域采样定理的理解。
(2) 熟悉时域离散系统的时域特性。
(3) 利用卷积方法观察分析系统的时域特性。
(4) 掌握序列傅里叶变换的计算机实现方法， 利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
2、实验原理与方法
(1) 时域采样。
(2) LTI系统的输入输出关系。
3、实验内容及步骤
(1) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。
(2) 编制实验用主程序及相应子程序。
① 信号产生子程序， 用于产生实验中要用到的下列信号序列：
a. xa(t)=A*e^-at *sin(Ω0t)u(t)
b. 单位脉冲序列：xb(n)=δ(n)
c. 矩形序列： xc(n)=RN(n), N=10
② 系统单位脉冲响应序列产生子程序。 本实验要用到两种FIR系统。
a. ha(n)=R10(n);
b. hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)
③ 有限长序列线性卷积子程序
用于完成两个给定长度的序列的卷积。 可以直接调用MATLAB语言中的卷积函数conv。 conv用于两个有限长度序列的卷积， 它假定两个序列都从n=0 开始。 调用格式如下：
y=conv (x, h)
4、实验结果分析
① 分析采样序列的特性。
a. 取采样频率fs=1 kHz,，即T=1 ms。
b. 改变采样频率，fs=300 Hz，观察|X(e^jω)|的变化，并做记录(打印曲线)；进一步降低采样频率，fs=200 Hz，观察频谱混叠是否明显存在，说明原因， 并记录(打印)这时的|X(e^jω)|曲线。
程序代码如下：
close all;clear all;clc;
A=50;
a=50*sqrt(2)*pi;
m=50*sqrt(2)*pi;
fs1=1000;
fs2=300;
fs3=200;
T1=1/fs1;
T2=1/fs2;
T3=1/fs3;
N=100;
n=[0:N-1];
x1=A*exp(-a*n*T1).*sin(m*n*T1);
x2=A*exp(-a*n*T2).*sin(m*n*T2);
x3=A*exp(-a*n*T3).*sin(m*n*T3);
w=linspace(-pi,pi,10000);
X1=x1*exp(-j*n'*w);
X2=x2*exp(-j*n'*w);
X3=x3*exp(-j*n'*w);
figure(1)
subplot(1,3,1)
plot(w/pi,abs(X1));
xlabel('\\omega/π');ylabel('|H(e^j^\\omega)|')
title('采样频率为1000Hz时的频谱图');
subplot(1,3,2)
plot(w/pi,abs(X2));
xlabel('\\omega/π');ylabel('|H(e^j^\\omega)|')
title('采样频率为300Hz时的频谱图');
subplot(1,3,3)
plot(w/pi,abs(X3));
xlabel('\\omega/π');ylabel('|H(e^j^\\omega)|')
title('采样频率为200Hz时的频谱图');
② 时域离散信号、 系统和系统响应分析。
a. 观察信号xb(n)和系统hb(n)的时域和频域特性； 利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n)， 比较所求响应y(n)和hb(n)的时域及频域特性， 注意它们之间有无差别， 绘图说明， 并用所学理论解释所得结果。
b. 观察系统ha(n)对信号xc(n)的响应特性。
程序代码如下：
close all;clear all;clc;
xbn=[1];
xcn=ones(1,10);
han=ones(1,10);
hbn=[1,2.5,2.5,1];
yn=conv(xbn,hbn);
n1=0:length(xbn)-1;
n2=0:length(hbn)-1;
subplot(2,1,1);stem(n1,xbn,'.')
xlabel('n');ylabel('xb(n)')
title('xb(n)的时域特性曲线')
subplot(2,1,2);stem(n2,hbn,'.')
xlabel('n');ylabel('hb(n)')
title('hb(n)的时域特性曲线')
figure(2)
subplot(2,1,1);
n1=[0:length(xbn)-1];
w=linspace(-pi,pi,10000);
Xbn=xbn*exp(-j*n1'*w);
subplot(2,1,1);
plot(w/pi,abs(Xbn));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[xb[n]的频谱');
n2=[0:length(hbn)-1];
w=linspace(-pi,pi,10000);
Hb=hbn*exp(-j*n2'*w);
subplot(2,1,2);
plot(w/pi,abs(Hb));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[hb(n)]的频谱');
figure(3)
n1=0:length(yn)-1;
n2=0:length(hbn)-1;
subplot(2,1,1);stem(n1,yn,'.')
xlabel('n');ylabel('y(n)')
title('y(n)的时域特性曲线')
subplot(2,1,2);stem(n2,hbn,'.')
xlabel('n');ylabel('hb(n)')
title('hb(n)的时域特性曲线')
figure(4)
subplot(2,1,1);
n1=[0:length(yn)-1];
w=linspace(-pi,pi,10000);
Y=yn*exp(-j*n1'*w);
subplot(2,1,1);
plot(w/pi,abs(Y));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[y(n)]的频谱');
n2=[0:length(hbn)-1];
w=linspace(-pi,pi,10000);
Hb=hbn*exp(-j*n2'*w);
subplot(2,1,2);
plot(w/pi,abs(Hb));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[hb(n)]的频谱');
zn=conv(xcn,han);
%以下为%系统ha(n)对信号xc(n)的频率响应特性的分析
figure(5)
n3=[0:length(zn)-1];
w=linspace(-pi,pi,10000);
Z=zn*exp(-j*n3'*w);
subplot(2,1,1);
plot(w/pi,abs(Z));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[zn]的幅度谱');
subplot(2,1,2);
plot(w/pi,angle(Z));
xlabel('\\omega/π');ylabel('相位')
title('DTFT[zn]的相位谱');
图一：xb(n)与hb(n)时域曲线的对比
图二：xb(n)与hb(n)频谱函数的对比
图三：y(n)与hb(n)时域曲线的对比
图四：y(n)与hb(n)频谱函数的对比
图五：系统ha(n)对信号xc(n)的频率响应特性的分析
包含z(n)幅度谱和相位谱
③ 卷积定理的验证。
（如下程序中Y1为直接线性卷积所得结果做DTFT，Y2为各自的DTFT相乘后所得的结果）
程序代码如下：
close all;clear all;clc;
x1=[1,3,4,0,2];
x2=[3,2,0,0,1,5];
n1=0:length(x1)-1;
n2=0:length(x2)-1;
y1=conv(x1,x2);
m1=0:length(y1)-1;
w=linspace(-pi,pi,10000);
X1=x1*exp(-j*n1'*w);
X2=x2*exp(-j*n2'*w);
Y1=y1*exp(-j*m1'*w);
Y2=X1.*X2;
subplot(2,1,1);
plot(w/pi,abs(Y1));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[y1(n)]的特性曲线');
subplot(2,1,2);
plot(w/pi,abs(Y2));
xlabel('\\omega/π');ylabel('幅度')
title('DTFT[y2(n)]的特性曲线');
5、思考题
(1)在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？
答：数字频率度量不相同，但他们所对应的模拟频率相同。由w=Ω*Ts得，采样间隔变化时模拟频率对应的数字频率会有相应的变化，故其度量会有所变化。
(2)在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？
答：有差别。DFT相当于序列频谱的等间隔采样，当取点少时，DFT包含序列频谱的信息会少，与序列频谱的误差会增大。取点多时，DFT会反映更多的频谱信息，误差会小，减轻了栅栏效应。二者因为误差而有差别。
实验二：用FFT作谱分析
1、实验目的
(1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法， 所以FFT的运算结果必然满足DFT的基本性质)。
(2) 熟悉FFT算法原理和FFT子程序的应用。
(3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法， 了解可能出现的分析误差及其原因， 以便在实际中正确应用FFT。
2、实验步骤
(1) 复习DFT的定义、 性质和用DFT作谱分析的有关内容。
(2) 复习FFT算法原理与编程思想， 并对照DIT-FFT运算流图和程序框图， 读懂本实验提供的FFT子程序。
(3) 编制信号产生子程序， 产生以下典型信号供谱分析用：
(4) 编写主程序。
(5) 按实验内容要求，上机实验， 并写出实验报告。
3、实验内容
(1) 对 2 中所给出的信号逐个进行谱分析。
程序代码如下：
close all;clear all;clc;
N=32;
n=0:999;
x1n=ones(1,4);
x2n=[1,2,3,4,4,3,2,1];
x3n=[4,3,2,1,1,2,3,4];
x4n=cos(0.25*pi*n);
x5n=sin(0.125*pi*n);
x6n=cos(8*pi*n)+cos(16*pi*n)+cos(20*pi*n);
X1=fft(x1n,N);
X2=fft(x2n,N);
X3=fft(x3n,N);
X4=fft(x4n,N);
X5=fft(x5n,N);
X6=fft(x6n,N);
k=0:N-1;
figure(1)
stem(k,abs(X1));
xlabel('k');ylabel('|X1(k)|');
title('DFT(x1n)的幅频特性图');
figure(2)
stem(k,abs(X2));
xlabel('k');ylabel('|X2(k)|');
title('DFT(x2n)的幅频特性图');
figure(3)
stem(k,abs(X3));
xlabel('k');ylabel('|X3(k)|');
title('DFT(x3n)的幅频特性图');
figure(4)
stem(k,abs(X4));
xlabel('k');ylabel('|X4(k)|');
title('DFT(x4n)的幅频特性图');
figure(5)
stem(k,abs(X5));
xlabel('k');ylabel('|X5(k)|');
title('DFT(x5n)的幅频特性图');
figure(6)
stem(k,abs(X6));
xlabel('k');ylabel('|X6(k)|');
title('DFT(x6n)的幅频特性图');
实验结果如下：
x1n到x6n的32点DFT幅频特性图
(2)令x(n)=x4(n)+x5(n)，用FFT计算8点和16点离散傅里叶变换。
close all;clear all;clc;
N=8;
M=16;
n=0:10000;
x4n=cos(0.25*pi*n);
x5n=sin(0.125*pi*n);
xn=x4n+x5n;
X1=fft(xn,N);
X2=fft(xn,M);
k1=0:N-1;
k2=0:M-1;
figure(1)
stem(k1,abs(X1));
xlabel('k');ylabel('|X1(k)|');
title('8点DFT[x(n)]的幅频特性图');
figure(2)
stem(k2,abs(X2));
xlabel('k');ylabel('|X2(k)|');
title('16点DFT[x(n)]的幅频特性图');
(3)令x(n)=x4(n)+j*x5(n)，用FFT计算8点和16点离散傅里叶变换。
close all;clear all;clc;
N=8;
M=16;
n=0:10000;
x4n=cos(0.25*pi*n);
x5n=sin(0.125*pi*n);
xn=x4n+j*x5n;
X1=fft(xn,N);
X2=fft(xn,M);
k1=0:N-1;
k2=0:M-1;
figure(1)
stem(k1,abs(X1));
xlabel('k');ylabel('|X1(k)|');
title('8点DFT[x(n)]的幅频特性图');
figure(2)
stem(k2,abs(X2));
xlabel('k');ylabel('|X2(k)|');
title('16点DFT[x(n)]的幅频特性图');
4、思考题
(1) N=8时，序列的幅频特性相同，因为序列x2(n)是1，2，3，4，4，3，2，1正好8个点。 x3(n)是4，3，2，1，1，2，3，4八个点，序列x2(（n－4)）8＝x3(n)，因此根据DFT循环移位，得幅频特性相同。而N=16时，通过对x2(n) 和x3(n)末端环移位特性，并且补零，两个序列不能通过移位相等，因此幅频特性不同。
(2) 如果周期信号的周期预先不知道，如何用FFT进行谱分析?
可以预先设定一个较小的N，进行N点FFT，再取2N，进行2N点FFT，如果二个FFT相等，则周期可能是N，直接取此结果作为FFT。，如果不等，取4N进行FFT，并与2N点FFT比较，依次进行，取2N 直到FFT近似相等为止，求的近似的FFT。
实验三：用窗函数法设计FIR数字滤波器
1、实验目的
(1) 掌握用窗函数法设计FIR数字滤波器的原理和方法。
(2) 熟悉线性相位FIR数字滤波器特性。
(3) 了解各种窗函数对滤波特性的影响。
2、实验内容及步骤
(1) 复习用窗函数法设计FIR数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。
(2) 编写程序。
① 编写能产生矩型窗、 升余弦窗、 改进升余弦窗和二阶升余弦窗的窗函数子程序。
可调用子函数如下程序给出
function wn=wn(n,num)
switch num
case 1
wn=(boxcar(n))'
case 2
wn=(hanning(n))'
case 3
wn=(hamming(n))'
case 4
wn=(blackman(n))'
end
end
下面用程序给出的是四种窗的单位脉冲响应波形和幅频特性曲线
以便于研究不同窗的特点
close all;clear all;clc;
N=;beita=2.5;
w1=boxcar(N);
w2=hanning(N);
w3=hamming(N);
w4=blackman(N);
wvtool(w1);
wvtool(w2);
wvtool(w3);
wvtool(w4);
②编写主程序。 其中幅度特性要求用dB表示。
在此程序中设N=10；wc=pi；且分别实验四种窗函数所得的滤波结果。
close all;clear all;clc;
n=10;
b=0:9
a=(n-1)/2;
wc=pi;
hdn=sin(wc*(b-a))./(pi*(b-a));
w1=wn(n,1);
w2=wn(n,2);
w3=wn(n,3);
w4=wn(n,4);
h1=hdn.*w1;
h2=hdn.*w2;
h3=hdn.*w3;
h4=hdn.*w4;
H1=fft(h1,300);
H2=fft(h2,300);
H3=fft(h3,300);
H4=fft(h4,300);
F1=abs(H1);
K1=20*log10(F1);
F2=abs(H2);
K2=20*log10(F2);
F3=abs(H3);
K3=20*log10(F3);
F4=abs(H4);
K4=20*log10(F4);
subplot(221);
plot(abs(K1));
title('矩形窗');
subplot(222);
plot(abs(K2));
title('升余弦窗');
subplot(223);
plot(abs(K3));
title('改进升余弦窗');
subplot(224);
plot(abs(K4));
title('二阶升余弦窗');
加四种窗结果分别是：
3、思考题
(1) 给定通带截止频率和阻带截止频率以及阻带最小衰减，用窗函数法设计线性相位低通滤波器的设计步骤：
①根据所需设计的数字滤波器类型，确定线性相位数字滤波器类型，本题为I型；
②选择合适的窗函数；
③确定理想低通数字滤波器的频率响应函数；
④计算理想低通数字滤波器的单位脉冲响应；
⑤加窗得到设计结果；
(2)用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为ω1和ω2，理想带通的单位脉冲响应hd(n)的求解过程：
