数字信号处理课程设计实验.docx
《数字信号处理课程设计实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计实验.docx（24页珍藏版）》请在冰点文库上搜索。
数字信号处理课程设计实验
华北电力大学实验报告
实验环境
MATLAB6.5
实验名称
实验一：
FFT的应用
实验目的
1、熟悉MATLAB在数字信号处理中的应用。
2、掌握利用FFT计算序列线性卷积的基本原理及编程实现。
3、掌握对连续信号进行采样的基本原理和方法，并利用FFT对信号进行频谱分析。
实验原理
1、离散时间信号的表示：
在数字信号处理中，所有信号都是离散时间信号——序列，所有的序列都可以表示为如下形式：
x（n）={…,x（-1）,x（0）,x
（1）,…}
用matlab的stem（n,x）即可实现简单的表示。
若要表示具有特定采样频率的信号，需要定义时间轴向量。
2、快速傅里叶变换：
采用时间抽取基2FFT算法。
MATLAB中提供了fft和ifft函数来分别计算DFT和IDFT。
fft和ifft函数是用机器语言，而不是用MATLAB指令写成的，因此它的执行速度很快。
fft函数的用法：
y=fft（x）;%计算x的快速离散傅里叶变换y
y=fft（x,N）;%计算x的N点FFT。
当x的长度大于N时，截断x；否则补零
ifft函数的用法：
y=ifft（x）;%计算x的快速离散傅里叶反变换y
y=ifft（x,N）;%计算x的N点IFFT
3、利用FFT计算线性卷积：
步骤如下：
1）将x1（n）和x2（n）都延长到N点，N=N1+N2-1
2）计算x1（n）的N点DFT，即：
X1（k）=FFT[x1（n）]
3）计算x2（n）的N点DFT,即：
X2（k）=FFT[x2（n）]
4）计算Y（k）=X1（k）X2（k）
5）计算的反变换,即y（n）=IFFT[X1（k）X2（k）]
4、利用FFT对信号进行谱分析：
包括振幅谱，相位谱和功率谱。
参数：
fs>=2f0
在matlab中，提供了计算模值的函数abs和计算相角的函数angle，其用法如下：
Abs（x）;
Angle（x）;
实验内容
1.对于两个序列：
x（n）=nR16（n）,h（n）=R8（n）
（1）在同一图形窗口中绘出两序列的时域图形。
（2）利用FFT编程计算两序列的线性卷积，绘出的时域图形。
设计方案：
首先定义x序列，为离散序列，再通过stem函数将x依次与y序列对应，即可表示两序列。
然后利用FFT计算俩序列的线性卷积，方法与步骤如实验原理所述，注意为使计算快速，可延长时把两序列都延长到最接近的2的整数幂的值，通过函数N=2^nextpow2（N1+N2-1）实现，然后两序列补零，求傅里叶变换，相乘，求反变换，依次即可实现。
源程序：
%两序列时域波形
N1=16;
x1=0:
N1-1;
y1=x1;
N2=8;
x2=0:
N2-1;
y2=ones（1,N2）;
subplot（2,1,1）;
stem（x1,y1,'b'）;holdon;stem（x2,y2,'r'）;
title（'两序列时域波形'）;
%线性卷积时域波形
N=2^nextpow2（N1+N2-1）;
y3=[y1（1:
N1）zeros（1,N-N1）];
y4=[y2（1:
N2）zeros（1,N-N2）];
Y1=fft（y3,N）;
Y2=fft（y4,N）;
Y=Y1.*Y2;
y=ifft（Y,N）;
n=0:
N-1;
subplot（2,1,2）;
stem（n,y）;
title（'线性卷积时域波形'）;
2．利用FFT对信号进行谱分析
对于连续信号xa（t）=cos（2πf1t）+5cos（2πf2t）+cos（2πf3t），其中f1=6.5kHz,f2=7kHz,f3=9kHz,以采样频率fs=32kHz对其进行采样，
（1）对xa（t）信号采集16点样本，分别作16点和补零到256点的FFT，并分别绘出对应的幅频特性曲线。
（2）对xa（t）信号采集256点样本，分别作256点和512点的FFT，并分别绘出对应的幅频特性曲线。
（3）比较
（1）和
（2）中的结果，分析采样点数和傅里叶变换点数对FFT的影响，说明高密度频谱和高分辨率频谱的特点与区别。
设计方案：
若要表示具有特定采样频率的信号，先需要根据采样频率定义时间轴向量，然后定义函数，依次表示出16个采样点，补零到256个采样点，256个采样点和补零至512个采样点时的时域波形和频谱，图形表示时直接使用stem（）函数即可，补零用x2=[x（1:
N1）zeros（1,N2-N1）]这样的形式，快速傅里叶变换为fft（）函数，幅值用abs（）函数直接求解。
源程序：
f1=6500;f2=7000;f3=9000;fs=32000;
t=0:
1/fs:
1;
x=cos（2*f1*pi*t）+5*cos（2*f2*pi*t）+cos（2*f3*pi*t）;
N1=16;N2=256;N3=512;
n=N3-1;
%16个采样点的时域波形，频谱
n1=0:
N1-1;
x1=x（1:
N1）;
subplot（2,4,1）;stem（n1,x1）;title（'16个采样点的时域信号x1'）;
Y1=fft（x1,N1）;
Y1=abs（Y1）;
k1=0:
N1-1;
f1=fs/N1*k1;
subplot（2,4,5）;stem（f1,Y1）;title（'信号x1的频谱'）;
%补零到256个采样点的时域波形，频谱
n2=0:
N1-1;
x2=[x（1:
N1）zeros（1,N2-N1）];
m1=0:
255;
subplot（2,4,2）;stem（m1,x2）;title（'补零到256点的时域信号x2'）;
Y2=fft（x2,N2）;
Y2=abs（Y2）;
k2=0:
N2-1;
f2=fs/N2*k2;
subplot（2,4,6）;
stem（f2,Y2）;
title（'信号x2的频谱'）;
%256个采样点的时域波形，频谱
n3=0:
N2-1;
x3=x（1:
N2）;
subplot（2,4,3）;stem（n3,x3）;title（'256个采样点的时域信号x3'）;
Y3=fft（x,N2）;
Y3=abs（Y3）;
k3=0:
N2-1;
f3=fs/N2*k3;
subplot（2,4,7）;stem（f3,Y3）;title（'信号x3的频谱'）;
%补零到512个采样点的时域波形，频谱
n4=0:
N2-1;
x4=[x（1:
N2）zeros（1,N3-N2）];
m2=0:
511;
subplot（2,4,4）;stem（m2,x4）;title（'补零到512点的时域信号x4'）;
Y4=fft（x4,N3）;
Y4=abs（Y4）;
k4=0:
N3-1;
f4=fs/N3*k4;
subplot（2,4,8）;stem（f4,Y4）;title（'信号x4的频谱'）;
实验结果及分析
1.对于两个序列：
x（n）=nR16（n）,h（n）=R8（n）
（1）在同一图形窗口中绘出两序列的时域图形。
（2）利用FFT编程计算两序列的线性卷积，绘出的时域图形。
实验结果：
实验分析：
两信号分别为x（n）=nR16（n）,h（n）=R8（n），为时间离散信号，两个信号同时在第一个图形中表示出来，结果与信号表达式相同，结果正确。
两信号卷积x（n）*h（n）为所有h（n-k）与x（k）的乘积的累加和。
所以在n=15时会出现卷积的最大值。
X（n）和h（n）的线性卷积的结果y（n）是一个有限序列，其非零值长度为N1+N2-1,最终matlab实现的结果为22个非零值，实验结果正确。
2．利用FFT对信号进行谱分析
实验结果：
实验分析：
程序运行结果如上图所示，可以得到如下结论：
当取样点N=16时，从频谱图上几乎无法看出信号的任何频率信息。
将16点中的信号补零到N=256点时，频谱的谱线相当密，但从中仍然很难看出信号的频率成分，称之为高密度频谱。
它只是在信号X1的频谱基础上增加采样密度，但不增加分辨率，无法提取有用的频谱成分。
当对序列取足够的采样点时，如256个采样点，可以从幅频特性中清晰地看出信号的频率成分，称之为高分辨率频谱。
频谱应为6.5kHz，7kHz，9kHz。
另一组为23kHz,25kHz,25.5kHz。
（如6.5k对称的为-6.5k采样频率为32k,32-6.5为25.5，故第二组为含25.5kHz的频谱，其余同理即可得，如上面第二三个图所示）
问题：
高密度频谱与高分辨率频谱有哪些特点和区别：
答：
高密度频谱与高分辨率频谱的谱线都非常密集，采样频率都非常高。
但高密度频谱是当时间域长度不变时，频域内对它的频谱进行提高采样频率，细化当前分辨率下的频谱，故对分辨率没有改善。
而高分辨率频谱才是我们所希望得到的频谱结果。
实验名称
设计性实验一：
IIR数字滤波器的设计
实验目的
1、本实验为设计性实验。
2、掌握用双线性变换法设计IIR数字滤波器的基本原理和设计方法。
3、掌握用双线性变换法设计IIR数字Butterworth滤波器的原理和设计方法。
实验原理
1、IIR数字滤波器设计的原理：
IIR数字滤波器的系统函数H（z）是z-1的有理分式，设计IIR数字滤波器，就是要确定系统函数H（z）的阶数N以及分子分母多项式的系数ak和bk，使其频率特性满足指定的要求。
用双线性变化法设计IIR数字滤波器的步骤如下：
（1）确定数字滤波器的通带截止频率wp、阻带截止频率ws、通带最大衰减Rp和阻带最小衰减Rs；
（2）对截止频率wp和ws进行频率预畸变；
（3）计算模拟滤波器的阶数N和频率Wc，进而求得s平面上的极点，从而构成模拟滤波器的传递函数Ha（s），使得满足给定的技术指标；
（4）利用双线性变换公式将模拟滤波器的传递函数Ha（s）变换成数字滤波器的系统函数H（z）；
（5）若需设计其他形式的滤波器，则根据频率变换公式将数字低通滤波器原型变换到相应形式的滤波器即可。
模拟域频率与数字域频率的关系为
Wc=（2/T）tan（w/2）
这是一种非线性的关系。
这种非线性关系使得模拟滤波器和数字滤波器的频率响应与对应频率的关系上发生了畸变，也造成了相位的非线性变化。
为保证各边界频率点为预先指定的频率，在确定模拟低通滤波器的系统函数之前，必须按上式进行频率预畸变。
2、Matlab中应用到的函数：
在MATLAB的数字信号处理工具箱中，提供了用双线性变换法设计数字Butterworth低通滤波器的相关函数。
（1）buttord：
Butterworth滤波器阶数选择函数
[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）;
其中，Wp为通带截止角频率，Ws为阻带截止角频率，Rp为通带最大衰减，Rs为阻带最小衰减；N是符合要求的滤波器最小阶数，Wn是Butterworth滤波器固有角频率（3dB）。
Wp、Ws、Wn均归一化在之间，单位为πrad/s；'s'表示用于模拟滤波器，去掉则用于数字滤波器。
（2）buttap：
Butterworth模拟低通滤波器的建立函数
[Z,P,K]=buttap（N）;
给出N阶Butterworth模拟滤波器的零点向量Z、极点向量P和增益K。
产生的滤波器在左半平面的单位圆附近有N个极点，没有零点。
（3）zp2tf：
零极点增益模型到传递函数模型的转换函数
[Bap,Aap]=zp2tf（Z,P,K）;
其中，Z、P、K分别为零极点增益模型的零点向量、极点向量和增益；Bap、Aap分别为传递函数分子和分母的多项式系数向量。
（4）lp2lp：
从低通原型向低通的转换函数
[b,a]=lp2lp（Bap,Aap,Wn）;
把截止角频率为1rad/s的模拟低通原型滤波器转换成截止角频率为Wn的模拟低通滤波器。
（5）bilinear：
双线性变换函数
[bz,az]=bilinear（b,a,Fs）;
用双线性变换法把模拟低通滤波器转换为数字低通滤波器。
其中，Fs是采样频率，bz和az分别是传递函数的分子和分母的多项式系数向量。
（6）freqz：
数字滤波器的频率响应函数
[H,W]=freqz（B,A,N）;
返回数字滤波器均匀分布在上的N点频率向量W和与之对应的N点频率响应向量H。
A和B分别是滤波器系统函数分子和分母的多项式系数向量；N最好选用2的整数幂，以便使用FFT快速运算，N的缺省值为512。
[H,F]=freqz（B,A,N,Fs）
对在[0,Fs/2]上等间隔采样N点，采样点频率及相应的频率响应值分别记录在F和H中。
实验内容
IIR数字滤波器的设计
用双线性变换法设计一个IIR数字Butterworth低通滤波器。
技术指标为：
通带截止频率fp=1kHz，阻带截止频率fs=1.5kHz，通带衰减Rp≤1dB，阻带衰减Rs≥40dB，采样频率Fs=10kHz。
绘出滤波器的幅频特性曲线和相频特性曲线，判断设计是否符合要求。
设计方案：
按照双线性变换法设计IIR滤波器的步骤，首先指标转换，w=2pi*f，因为后面用到buttord（）函数，所以指标转换时归一化处理。
对wp,ws进行预畸变，将数字域频率转化成模拟域频率，利用buttord函数求滤波器阶数，通过函数由模拟滤波器的零点向量，极点向量，增益求模拟低通的传递函数。
再通过双线性变换函数bilinear（）将模拟低通转换为数字低通。
通过freqz函数求数字滤波器的频率响应函数即可。
源程序：
clearall;
fp=1000;fs=1500;Fs=10000;
Rp=1.5;Rs=25;
%指标转换
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Fs=Fs/Fs;
%预畸变
Wp=2*tan（wp/2）;
Ws=2*tan（ws/2）;
%由阶数，零点向量，极点向量，增益实现模拟低通的传递函数
[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）;
[z,p,k]=buttap（N）;
[Bap,Aap]=zp2tf（z,p,k）;
%由低通原型转向低通
[b,a]=lp2lp（Bap,Aap,Wn）;
%双线性变换实现模拟低通变为数字低通
[bz,az]=bilinear（b,a,Fs）;
%频率响应
[h,w]=freqz（bz,az,512,Fs）;
%输出幅频
h1=abs（h）;
subplot（2,1,1）;
h2=20*log10（（h1+eps）/max（h1））;
plot（w*10000,h2）;
title（'幅频特性'）;
gridon;
%输出相频
subplot（2,1,2）;
plot（w*10000,angle（h））;
title（'相频特性'）;
gridon;
实验结果及分析
实验结果：
实验分析：
上图为，程序运行后滤波器的幅度特性和相频特性。
从图中可以看出：
幅度响应曲线在通带截止频率（f=1000hz）即在通带内有更大部分的幅度接近于0，可以从图中看出，通带截止频率为1000Hz，阻带截止频率1500Hz为下降3dB是的频率值，纵坐标换成dB时可清楚看到。
幅频特性图中可以清晰的看到数字滤波器的通带、过渡带、阻带三部分。
IIR数字滤波器不能做到严格的线性相位。
所以结果符合设计要求。
问题：
冲激响应不变法和双线性变换法比较：
答：
相同点：
首先设计模拟滤波器，再将模拟滤波器转换为数字滤波器。
不同点：
冲激响应不变法的幅频响应有失真，频率之间呈线性。
双线性变换法的幅频响应无混迭，频率之间有失真。
实验名称
设计性实验二：
FIR数字滤波器的设计
实验目的
1、本实验为设计性实验。
2、掌握用窗函数法设计IIR数字滤波器的基本原理和设计方法。
3、掌握用窗函数法设计线性相位FIR数字低通滤波器的编程实现。
实验原理
1、设计原理：
FIR数字滤波器的独特优点是容易得到严格的线性相位，所以FIR总是最稳定的。
FIR数字滤波器主要采用非递归结构，而且可以借助FFT计算。
因为所要设计的是FTR数字滤波器，其单位冲激响应h（n）必然是有限长的，所以要用有限长的h（n）来逼近无限长的hd（n）。
最有效的方法昰截断hd（n），即用有限长的窗函数w（n）来截取hd（n），表示h（n）=hd（n）w（n）,这种方法称为窗函数法。
用窗函数法设计FIR数字滤波器的步骤如下：
（1）给出所要设计的FIR数字滤波器的技术指标，如通带截止频率wp、阻带截止频率ws、通带衰减Rp和阻带衰减As；
（2）根据允许的过渡带宽度及阻带衰减，初步选择窗函数和N值。
（3）若选用理想低通逼近，求出理想低通的冲激响应hd（n）。
理想低通的截止频率选择为wc=（wp+ws）/2
（4）将hd（n）与窗函数相乘得FIR数字滤波器的冲激响应h（n）
（5）计算FIR数字滤波器的频率响应，并验证是否达到所要求的指标。
2、matlab中用到的函数：
在MATLAB的数字信号处理工具箱中提供了用窗函数法设计线性相位FIR数字滤波器的函数fir1，它具有标准低通、带通、高通和带阻等类型。
fir1：
滤波器设计函数
B=fir1（N,Wn,'ftype',window）;
其中，N为FIR滤波器的阶数，Wn为滤波器截止频率，取值范围【0,1】rad/s，对于带通、带阻滤波器，Wn=[W1,W2]，且W1B为FIR数字滤波器的系数向量，长度为N+1。
MATLAB中提供了很多窗函数，常用的有rectwin、bartlett、hanning、hamming、blackman、kaiser。
实验内容
FIR数字滤波器的设计
用窗函数法设计一个线性相位FIR数字低通滤波器。
技术指标为：
通带截止角频率ωp=0.2π，阻带截止角频率ωs=0.3π，通带衰减Rp≤1dB，阻带衰减Rs≥40dB。
绘出滤波器的幅频特性曲线和相频特性曲线，判断设计是否符合要求。
根据相同的滤波器要求，选用不同的窗函数进行设计，比较各种窗函数对FIR数字滤波器频率特性的影响。
设计方案：
Wp为通带截止频率，Ws为阻带截止频率，W为过渡带带宽。
阶数N通过函数N=ceil（6.6*pi/W）+1求，也可以用主瓣过渡带宽度来求，如汉明窗为8pi/N。
首先产生理想低通滤波器的单位冲激响应hd，然后单位冲激响应加海明窗，直接用hamming函数实现。
理想低通与海明窗相乘得到截断后实际单位脉冲响应。
通过freqz函数计算实际滤波器的幅度响应和相位响应即可。
源程序：
1、海明窗
Wp=0.2*pi;
Ws=0.3*pi;
W=Ws-Wp;
N=ceil（6.6*pi/W）+1;
%理想低通滤波器的单位冲激响应
n=0:
1:
N-1;
Wc=（Ws+Wp）/2;
n1=（N-1）/2;
n=0:
1:
N-1;
m=n-n1+eps;
hd=sin（Wc*m）./（pi*m）;
subplot（2,3,1）;stem（n,hd）;title（'理想单位脉冲响应hd（n）'）;
gridon;
%加海明窗
w_ham=（hamming（N））';
subplot（2,3,2）;stem（n,w_ham）;title（'海明窗'）;
gridon;
%实际的单位脉冲响应
h=hd.*w_ham;
subplot（2,3,4）;stem（n,h）;title（'实际单位脉冲响应h（n）'）;
gridon;
%计算实际滤波器的幅度响应
[H,w]=freqz（h）;
mag=abs（H）;
db=20*log10（mag）;
subplot（2,3,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;
gridon;
%计算实际滤波器的相位响应
pha=angle（H）;
subplot（2,3,6）;plot（w/pi,pha）;title（'相位响应'）;
gridon;
2、凯泽窗
Wp=0.2*pi;
Ws=0.3*pi;
W=Ws-Wp;
Rs=60;
N=ceil（（Rs-7.95）/（14.36*W/（2*pi））+1）+1;
%理想低通滤波器的单位冲激响应
n=0:
1:
N-1;
beta=0.1102*（Rs-8.7）;%线性相位斜率
Wc=（Ws+Wp）/2;
n1=（N-1）/2;
n=0:
1:
N-1;
m=n-n1+eps;
hd=sin（Wc*m）./（pi*m）;
subplot（2,3,1）;stem（n,hd）;title（'理想单位脉冲响应hd（n）'）;
gridon;
%加凯泽
w_kai=（kaiser（N,beta））';
subplot（2,3,2）;stem（n,w_kai）;title（'凯泽窗'）;
gridon;
%实际的单位脉冲响应
h=hd.*w_kai;
subplot（2,3,4）;stem（n,h）;title（'实际单位脉冲响应h（n）'）;
gridon;
%计算实际滤波器的幅度响应
[H,w]=freqz（h）;
mag=abs（H）;
db=20*log10（mag）;
subplot（2,3,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;
gridon;
%计算实际滤波器的相位响应
pha=angle（H）;
subplot（2,3,6）;plot（w/pi,pha）;title（'相位响应'）;
gridon;
实验结果及分析
实验结果：
实验分析：
根据实验的技术指标，观察幅度响应（dB）的图形，可以看出通带截止角频率wp为0.2pi，阻带截止角频率ws为0.3pi，由图中可以观察到过渡带在0.2pi和0.3pi之间。
通带衰减Rp非常小，阻带衰减Rs可以看出>40dB，与实验的技术指标要求相同。
观察相频特性，可以看出FIR滤波器有更严格的线性相位。
通过比较如上两个窗函数—海明窗和凯泽窗，可以得到如下结论：
（1）凯泽窗的最小阻带衰减要比海明窗小，
（2）得到相同的主瓣过渡带凯泽窗的阶数（窗函数长度）比海明窗小，（3）旁瓣峰值幅度凯泽窗要比海明窗小更理想。
问题：
窗函数的长度N对FIR滤波器的性能有什么影响？
答：
窗函数的宽度N越大，窗函数的频谱的主瓣越窄，因而过渡带也越窄。
问题：
吉布斯效应是怎么产生的，增加窗函数长度能否消除？
为什么？
答：
吉布斯效应是由窗函数的频谱的旁瓣造成的。
在一般情况下，对窗函数的要求是：
（1）尽量减少窗函数频谱的旁瓣高度，使能量集中在主瓣，减少通带/阻带中的波纹。
（2）主瓣的宽度尽量窄，获得较陡的过渡带。
窗函数的长度N不能改变主瓣和旁瓣的相对比例，增大N可以减小波纹幅度，但效果并不是很好。
问题：
不同的窗函数对FIR滤波器的性能有什么影响？
答：
不同的窗函数对滤波器旁瓣峰值幅度、主瓣过渡带幅度、最小阻带衰减等方面均有不同的特性。
实验
升级会员 