南邮 数字信号处理实验报告带问题答案小结.docx
《南邮 数字信号处理实验报告带问题答案小结.docx》由会员分享，可在线阅读，更多相关《南邮 数字信号处理实验报告带问题答案小结.docx（39页珍藏版）》请在冰豆网上搜索。
南邮数字信号处理实验报告带问题答案小结
南京邮电大学
实验报告
实验名称熟悉MATLAB环境
快速傅里叶变换（FFT）及其应用
IIR数字滤波器的设计
FIR数字滤波器的设计
课程名称数字信号处理A
班级学号_12006311____
姓名_______张文欣_____________
开课时间2014/2015学年，第二学期
实验一熟悉MATLAB环境
一、实验目的
（1）熟悉MATLAB的主要操作命令。
（2）学会简单的矩阵输入和数据读写。
（3）掌握简单的绘图命令。
（4）用MATLAB编程并学会创建函数。
（5）观察离散系统的频率响应。
二、实验内容
（1）数组的加、减、乘、除和乘方运算。
输入A=[1234]，B=[3,4,5,6]，求
C=A+B，D=A-B，E=A.*B,F=A./B,G=A.^B。
并用stem语句画出A、B、C、D、
E、F、G。
n=0:
1:
3;
A=[1234];
subplot（4,2,1）
stem（n,A）
xlabel（'n'）
ylabel（'A'）
B=[3,4,5,6];
subplot（4,2,2）
stem（n,B）
xlabel（'n'）
ylabel（'B'）
C=A+B;
subplot（4,2,3）
stem（n,C）
xlabel（'n'）
ylabel（'C'）
D=A-B
Subplot（4,2,4）
stem（n,D）
xlabel（'n'）
ylabel（'D'）
E=A.*B
subplot（4,2,5）
stem（n,E）
xlabel（'n'）
ylabel（'E'）
F=A./B
subplot（4,2,6）
stem（n,F）
xlabel（'n'）
ylabel（'F'）
G=A.^B
subplot（4,2,7）
stem（n,G）
xlabel（'n'）
ylabel（'G'）
（2）用MATLAB实现下列序列：
a）
n=0:
1:
15;
x1=0.8.^n;
stem（n,x1）
xlabel（'n'）
ylabel（'x（n）'）
title（'2（a）'）
b）
n=0:
1:
15;
i=sqrt（-1）;
a=0.2+3*i;
x2=exp（a*n）;
figure
subplot（1,2,1）
stem（n,real（x2））
xlabel（'n'）
ylabel（'x（n）实部'）
subplot（1,2,2）
stem（n,imag（x2））
xlabel（'n'）
ylabel（'x（n）虚部'）
c）
n=0:
1:
15;
x3=3*cos（0.125*pi*n+0.2*pi）+2*sin（0.25*pi*n+0.1*pi）;
stem（n,x3）
xlabel（'n'）
ylabel（'x（n）'）
title（'2（c）'）
（4）绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注：
a）
t=0:
0.001:
10;
x=sin（2*pi*t）;
plot（t,x,'r-'）
xlabel（'t'）,ylabel（'x（t）'）,title（'sin（2\pit）'）
b）
t=0:
0.001:
4;
x=cos（100*pi*t）.*sin（pi*t）;
plot（t,x,'b-'）
xlabel（'t'）,ylabel（'x（t）'）,title（'cos（100\pit）*sin（\pit）'）
（6）给定一因果系统
，求出并绘制H（z）的幅频响应和相频响应。
b=[1,sqrt
（2）,1];a=[1,-0.67,0.9];
[h,w]=freqz（b,a）;
am=20*log10（abs（h））;%am=20*log10（abs（h））为幅频响应取dB
subplot（2,1,1）;
plot（w,abs（h））;xlabel（'w'）;ylabel（'幅频响应'）;title（'系统响应'）
ph=angle（h）;
subplot（2,1,2）;
plot（w,ph）;xlabel（'w'）;ylabel（'相频响应'）;
（7）计算序列{8-2-123}和序列{23-1-3}的离散卷积，并作图表示卷积结果。
N=50;
a=[1-2];
b=[10.1-0.06];
x=[1zeros（1,N-1）];
k=0:
1:
N-1;
y=filter（a,b,x）;
stem（k,y）;
xlabel（'n'）;ylabel（'幅度'）;
（8）求以下差分方程所描述系统的单位脉冲响应h（n）,
N=50;
a=[1-2];
b=[10.1-0.06];
x=[1zeros（1,N-1）];
k=0:
1:
N-1;
y=filter（a,b,x）;
stem（k,y）;
xlabel（'n'）;ylabel（'幅度'）;
每一小题均给出实验过程与结果（含实验程序、运行的数据结果和图形）；
实验二快速傅里叶变换（FFT）及其应用
一、实验目的
（1）在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。
（2）应用FFT对典型信号进行频谱分析。
（3）了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。
（4）应用FFT实现序列的线性卷积和相关。
二、实验内容
实验中用到的信号序列
a） 高斯序列
b） 衰减正弦序列
c） 三角波序列
d） 反三角波序列
（1）观察高斯序列的时域和幅频特性，固定信号
中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？
记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。
n=0:
15;
p=8;
q=2;
xa=exp（-（n-p）.^2/q）;
subplot（3,2,1）;
stem（n,xa,'.'）;
title（'q=2高斯序列时域特性'）;
subplot（3,2,2）;
stem（fft（xa）,'.'）;
title（'q=2高斯序列幅频特性'）;
n=0:
15;
p=8;
q=4;
xa=exp（-（n-p）.^2/q）;
subplot（3,2,3）;
stem（n,xa,'.'）;
title（'q=4高斯序列时域特性'）;
subplot（3,2,4）;
stem（fft（xa）,'.'）;
title（'q=4高斯序列幅频特性'）;
n=0:
15;
p=8;
q=8;
xa=exp（-（n-p）.^2/q）;
subplot（3,2,5）;
stem（n,xa,'.'）;
title（'q=8高斯序列时域特性'）;
subplot（3,2,6）;
stem（fft（xa）,'.'）;
title（'q=8高斯序列幅频特性'）;
由上图比较可知，当保持P不变时，随q的增大时域波形展宽且变得平滑，而频域波形且频域分量变少，不容易发生混叠，可见高斯序列中的q表示时域波形的陡峭程度。
n=0:
15;
p=8;
q=8;
xa=exp（-（n-p）.^2/q）;
subplot（3,2,1）;
stem（n,xa,'.'）;
title（'p=8高斯序列时域特性'）;
subplot（3,2,2）;
stem（fft（xa）,'.'）;
title（'p=8高斯序列幅频特性'）;
n=0:
15;
p=13;
q=8;
xa=exp（-（n-p）.^2/q）;
subplot（3,2,3）;
stem（n,xa,'.'）;
title（'p=13高斯序列时域特性'）;
subplot（3,2,4）;
stem（fft（xa）,'.'）;
title（'p=13高斯序列幅频特性'）;
n=0:
15;
p=14;
q=8;
xa=exp（-（n-p）.^2/q）;
subplot（3,2,5）;
stem（n,xa,'.'）;
title（'p=14高斯序列时域特性'）;
subplot（3,2,6）;
stem（fft（xa）,'.'）;
title（'p=14高斯序列幅频特性'）;
由上图比较可知当保持q不变时，随p的增大时域波形整体向右移动，由上图可见时域的窗口固定时在p=13,p=14有几点被泄露了，可见p值表示时域波形峰值的位置，p值越大泄露越多，而频域波形随p的增大频域分量会增多，容易产生混叠。
（3）观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列
和
的幅频特性，观察两者的序列形状和频谱曲线有什么异同？
绘出两序列及其幅频特性曲线。
在
和
末尾补零，用N=32点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？
两种情况的FFT频谱还有相同之处吗？
这些变化说明了什么？
n=[0:
3];k=[1:
8];
Xc（n+1）=n;Xc（n+5）=4-n;
Xd（n+1）=4-n;Xd（n+5）=n;
%{三角波特性}%
subplot（2,2,1）;plot（k-1,Xc）;
xlabel（'n'）;ylabel（'时域特性'）;text（1,3,'三角波'）;
subplot（2,2,2）;plot（k-1,abs（fft（Xc）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'三角波'）;
%{反三角波特性}%
subplot（2,2,3）;plot（k-1,Xd）;
xlabel（'n'）;ylabel（'时域特性'）;text（3,3,'反三角波'）;
subplot（2,2,4）;plot（k-1,abs（fft（Xd）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'反三角波'）;
legend（'N=8'）
由图可知，N=8时，正反三角的频域图形是相同的，因为做DFT时要先周期延拓，做完后取主值部分，而正反三角形周期延拓后是相同的，只差一个相位，因此得到的频域图形也是相同的。
n=[0:
3];k=[1:
32];
Xc（n+1）=n;Xc（n+5）=4-n;
Xd（n+1）=4-n;Xd（n+5）=n;
Xc（9:
32）=0;Xd（9:
32）=0;k=1:
32;figure;
%{三角波特性}%
subplot（2,2,1）;plot（k-1,Xc）;
xlabel（'n'）;ylabel（'时域特性'）;text（1,3,'三角波'）;
subplot（2,2,2）;plot（k-1,abs（fft（Xc）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'三角波'）;
%{反三角波特性}%
subplot（2,2,3）;plot（k-1,Xd）;
xlabel（'n'）;ylabel（'时域特性'）;text（3,3,'反三角波'）;
subplot（2,2,4）;plot（k-1,abs（fft（Xd）））;
xlabel（'k'）;ylabel（'幅频特性'）;text（4,10,'反三角波'）;
legend（'N=32'）
N=32时，两者的频谱是不同的，因为此时再做周期延拓就不相同了。
在后面补零对于正三角形在n=8时是连续的，而反三角波在n=8时有一个突变，时域中出现了陡峭的地方，在频域中频谱分量会增多。
通过N=8和N=32比较得，通过在原序列的末端补零，增加了采样的点数，是谱线增多，弱化了栅栏效应，但增多后的谱线形状是与时域信号的形状有关的。
（5）用FFT分别实现
（p＝8，q＝2）和
（a＝0.1，f＝0.0625）的16点循环卷积和线性卷积。
n1=0:
1:
15;
p1=8;
q1=2;
x=exp（-（n1-p1）.^2/q1）;
n2=0:
1:
15;
a=0.1;
f2=0.0625;
y=（exp（-a*n2））.*sin（2*pi*f2*n2）;
N=length（x）;
n=0:
N-1;
n3=0:
30;
X=fft（x）;
Y=fft（y）;
x32=[xzeros（1,16）];
y32=[yzeros（1,16）];
X32=fft（x32）;
Y32=fft（y32）;
z16=ifft（X.*Y）;
z32=ifft（X32.*Y32）;
subplot（2,2,1）;
plot（n,z16,'-*'）;
xlabel（'n'）;
ylabel（'z（n）'）;
title（'循环卷积结果'）;
subplot（2,2,2）;
plot（n3,z32（1:
2*N-1）,'-o'）;
xlabel（'n'）;
ylabel（'z（n）'）;
title（'线性卷积结果'）;
rm16=real（ifft（conj（X）.*Y））;
rm32_0=real（ifft（conj（X32）.*Y32））;
rm32=[rm32_0（N+2:
2*N）rm32_0（1:
N）];
m=n;
subplot（2,2,3）;
plot（m,rm16,'--'）;
xlabel（'m'）;
ylabel（'rm'）;
title（'循环相关结果'）;
m=-（N-1）:
N-1;
subplot（2,2,4）;
plot（m,rm32,'--o'）;
xlabel（'m'）;
ylabel（'rm'）;title（'线性相关结果'）;
每一小题均给出实验过程与结果（含实验程序、运行的数据结果和图形）；
实验三IIR数字滤波器的设计
一、实验目的
（1）掌握双线性变换法及脉冲响应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。
（2）观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。
（3）熟悉巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器的频率特性。
二、实验内容
（1）P162例4.4
设采样周期T=250
（采样频率
=4kHz），分别用脉冲响应不变法和双线性变换法设计一个三阶巴特沃思低通滤波器，其3dB边界频率为
=1kHz。
%脉冲响应不变法
fc=1000;
fs=4000;
[B,A]=butter（3,2*pi*fc,'s'）;
[num1,den1]=impinvar（B,A,fs）;
[h1,w]=freqz（num1,den1）;
f=w/pi*fs/2;
%双线性变换法
fc=1000;
fs=4000;
[B,A]=butter（3,2*fs*tan（pi*fc/fs）,'s'）;
[num2,den2]=bilinear（B,A,fs）;
[h2,w]=freqz（num2,den2）;
f=w/pi*fs/2;
plot（f,abs（h1）,'-.',f,abs（h2）,'-'）;
gridon
legend（'脉冲响应不变法','双线性变换法'）
xlabel（'频率/Hz'）
ylabel（'幅度'）
（2）
=0.2kHz,
=1dB,
=0.3kHz,At=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。
比较这两种方法的优缺点。
clear;
fc=200;
Rp=1;
fr=300;
Rs=25;
T=0.001;
fs=1/T;
%脉冲响应不变法
Wp=2*pi*fc;
Ws=2*pi*fr;
[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）
[B,A]=butter（N,Wn,'s'）;
[num1,den1]=impinvar（B,A,fs）;
[h1,w]=freqz（num1,den1）;
%双线性变换法
Wp=2*fs*tan（pi*fc/fs）;
Ws=2*fs*tan（pi*fr/fs）;
[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）
[B,A]=butter（N,Wn,'s'）;
[num2,den2]=bilinear（B,A,fs）;
[h2,w]=freqz（num2,den2）;
f=w/pi*fs/2;
plot（f,abs（h1）,'-.',f,abs（h2）,'-'）;
gridon
legend（'脉冲响应不变法','双线性变换法'）
xlabel（'频率/Hz'）
ylabel（'幅度'）
N=
9
Wn=
1.3693e+003
N=
6
Wn=
1.7043e+003
由图可以看出，利用脉冲响应不变法得到的滤波器在高频部分（阻带）衰减特性不好；而双线性变换法将模拟频率的无穷大映射到数字频率的pi,这样使得其对应的滤波器在高频部分衰减很快。
（3）利用双线性变换法设计满足下列指标的巴特沃思数字低通滤波器，并作图验证设计结果：
=1.2kHz,
=2kHz,At
=8kHz。
clear;
fc=1200;
Rp=0.5;
fr=2000;
Rs=40;
fs=8000;
Wp=2*fs*tan（pi*fc/fs）;
Ws=2*fs*tan（pi*fr/fs）;
[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）
[B,A]=butter（N,Wn,'s'）;
[num2,den2]=bilinear（B,A,fs）;
[h2,w]=freqz（num2,den2）;
f=w/pi*fs/2;
plot（f,abs（h2）,'-'）;
gridon
legend（'双线性变换法'）
xlabel（'频率/Hz'）
ylabel（'幅度'）
N=
9
Wn=
9.5918e+003
每一小题均给出实验过程与结果（含实验程序、运行的数据结果和图形）；
实验四FIR数字滤波器的设计
一、实验目的
（1）掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉相应的计算机编程；
（2）熟悉线性相位FIR滤波器的幅频特性和相频特性；
（3）了解各种不同窗函数对滤波器性能的影响。
二、实验内容
（1）N=45，计算并画出矩形窗、汉明窗、布莱克曼窗的归一化的幅度谱。
N=45;
n=0:
1:
N-1;
win=boxcar（N）;%矩形窗
[H,w]=freqz（win,1）;
H=H/H
（1）;
stem（n,win）
figure
plot（w,abs（H））;
N=45;
n=0:
1:
N-1;
win=hamming（N）;%汉明窗
[H,w]=freqz（win,1）;
H=H/H
（1）;
stem（n,win）
figure
plot（w,abs（H））;
N=45;
n=0:
1:
N-1;
win=blackman（N）;%布莱克曼窗
[H,w]=freqz（win,1）;
H=H/H
（1）;
stem（n,win）
figure
plot（w,abs（H））;
（2）用矩形窗设计一个21阶的线性相位低通FIR数字滤波器，截止频率Wc＝0.25π，求出滤波器系数，并绘出滤波器的幅频特性。
修改程序，分别得到阶次为N＝41，61的滤波器，并显示其各自的幅频曲线。
N=21
矩形窗
Wc=0.25*pi;
N=21;
M=（N-1）/2;%位移量
forn=0:
（N-1）
if（n==fix（M））%中间的点单独算
hd（n+1）=Wc/pi;
else
hd（n+1）=sin（Wc*（n-M））/（pi*（n-M））;
end;
end;
win=boxcar（N）;
%%%不同窗函数
h=hd.*win';
[H,w]=freqz（h,1）;
n=0:
1:
N-1;
subplot（3,1,1）;
stem（n,h）
xlabel（'n'）,ylabel（'h（n）'）
subplot（3,1,2）;
plot（w,abs（H））;
xlabel（'\omega'）,ylabel（'幅度谱'）
subplot（3,1,3）;
plot（w,angle（H））;
xlabel（'\omega'）,ylabel（'相位谱'）
legend（'N=21矩形窗'）
汉明窗
Wc=0.25*pi;
N=21;
M=（N-1）/2;%位移量
forn=0:
（N-1）
if（n==fix（M））%中间的点单独算
hd（n+1）=Wc/pi;
else
hd（n+1）=sin（Wc*（n-M））/（pi*（n-M））;
end;
end;
win=hamming（N）;%汉明窗
h=hd.*win';
[H,w]=freqz（h,1）;
n=0:
1:
N-1;
subplot（3,1,1）;
stem（n,h）
xlabel（'n'）,ylabel（'h（n）'）
subplot（3,1,2）;
plot（w,abs（H））;
xlabel（'\omega'）,ylabel（'幅度谱'）
subplot（3,1,3）;
plot（w,angle（H））;
xlabel（'\omega'）,ylabel（'相位谱'）
legend（'N=21汉明窗'）
N=41时
矩形窗
Wc=0.25*pi;
N=41;
M=（N-1）/2;%位移量
forn=0:
（N-1）
if（n==fix（M））%中间的点单独算
hd（n+1）=Wc/pi;
else
hd（n+1）=sin（Wc*（n-M））/（pi*（n-M））;
end;
end;
win=boxcar（N）;
%%%不同窗函数
h=hd.*win';
[H,w]=freqz（h,1）;
n=0:
1:
N-1;
subplot（3,1,1）;
stem（n,h
升级会员 