第2学年《数字信号处理实验讲义》要点Word格式文档下载.docx
《第2学年《数字信号处理实验讲义》要点Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《第2学年《数字信号处理实验讲义》要点Word格式文档下载.docx（24页珍藏版）》请在冰豆网上搜索。
在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换（DFT）。
这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x（n）的长度为N时，它的DFT和反变换定义分别为：
有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。
而快速傅立叶变换
并不是与
不同的另外一种变换，而是为减少
计算次数的一种快速有效的算法。
这种快速算法，主要是利用了
下面两个特性使长序列的
分解为更小点数的
所实现的。
2、用FFT计算线性卷积
用FFT可以实现两个序列的圆周卷积。
在一定的条件下，可以使圆周卷积等于线性卷积。
一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度N≥N1＋N2-1，对于长度不足N的两个序列，分别将他们补零延长到N。
当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。
有两种方法：
重叠相加法、重叠保留法。
3、利用FFT求解离散傅立叶变换（IFFT）
逆变换IFFT为
求X（k）的逆FFT可以分为以下3个步骤：
（1）取X（k）的共扼，得到
；
（2）求
的FFT，得N
（3）取
的共扼，并除以N，即得到x（n）。
采用这种方法可以直接用FFT程序来计算逆FFT。
有关IFFT的具体应用，与FFT一一对应，在此不再赘述。
4、MATLAB相关函数、程序
MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有fft、Ifft、fft2、ifft2,fftn、ifftn和fftshift、Ifftshift等。
当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。
所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。
fft和ifft函数一维快速正傅里叶变换和逆傅里叶变换
调用方式:
Y=fft（X）
参数说明
·
如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换；
如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；
如果X是（N
维数组，则是对第一个非单元素的维进行离散傅立叶变换；
Y＝fft（X,N）
N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。
Y＝fft（X,[],dim）或Y＝fft（X,N,dim）
在参数dim指定的维上进行离散傅立叶变换；
当X为矩阵时，dim用来指定变换的实施方向：
dim=1,表明变换按列进行；
dim=2表明变换按行进行。
函数Ifft的参数应用与函数Fft完全相同。
应用说明
【实例1】fft的应用
X=[2128];
Y=fft（X）
运行结果
Y＝
13.00000+7.0000i-5.00000-7.0000
【实例2】fft（X,N,dim）的应用
A=[2578;
1405;
3851;
9127];
Z=fft（A,[],1）
【实例3】fft在信号分析中的应用
使用频率分析方法从受噪声污染的信号x（t）中鉴别出有用的信号。
程序：
t=0:
0.001:
1;
%采样周期为0.001s,即采样频率为1000Hz;
%产生受噪声污染的正县正弦波信号；
x=sin（2*pi*100*t）+sin（2*pi*200*t）+rand（size（t））;
subplot（2,1,1）
plot（x（1:
50））;
%画出时域内的信号；
Y=fft（x,512）;
%对X进行512点的傅立叶变换；
f=1000*（0:
256）/512;
%设置频率轴（横轴）坐标，1000为采样频率；
subplot（2,1,2）
plot（f,Y（1:
257））;
%画出频域内的信号
图1-1时域信号和频域信号的比较
由上面的图1-1可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。
但是通过对x（t）作傅立叶变换，把时域信号变换到频域进行分析，可以明显看出信号中100Hz和200Hz的两个频率分量。
fft2和ifft2函数调用方式
1）Y=fft2（X）
如果X是向量，则此傅立叶变换即变成一维傅立叶变换fft；
如果X是矩阵，则是计算该矩阵的二维快速傅立叶变换；
数据二维傅立叶变换fft2（X）相当于fft（fft（X）
，即先对X的列做一维傅立叶变换，然后再对变换结果的行做一维傅立叶变换。
2）Y＝fft2（X,M,N）
通过对X进行补零或截断，使得成为（M*N）的矩阵。
函数Ifft2的参数应用与函数Fft2完全相同。
【实例6-4】fft2、ifft2的应用
A=[25789;
13750;
26149;
81526];
Y=fft2（A）;
B=ifft2（Y）;
运行结果：
四、实验内容
1、分别对下列序列进行频谱分析：
编制DFT程序及FFT程序，并比较DFT程序与FFT程序的运行时间。
（1）实指数序列
（2）复指数序列
（3）周期为N的正弦序列
（4）周期为N的余弦序列
（5）矩形序列
从以上五种输入序列中选择3种作为实验内容进行实验。
其中（3）、（4）只能选其中一个。
2、序列
，试绘制
及它的离散傅里叶变换
图。
五、思考题
对一个有限长序列进行DFT等价于将该序列周期延拓后进行DFS展开，因为DFS也只是取其中一个周期来计算，所以FFT在一定条件下也可以用以分析周期信号序列。
如果实正弦信号
用16点FFT来做DFS运算，得到的频谱是信号本身的真实谱吗？
为什么？
六、实验报告要求
1、简述实验目的及原理。
2、按实验步骤附上实验信号序列和幅频特性曲线，分别分析所得到的图形，说明参数改变对时域和频域的影响。
3、总结实验中的主要结论。
4、简要回答思考题。
实验二IIR数字滤波器的设计
1、掌握IIR数字滤波器的设计原理、设计方法和设计步骤；
2、能根据给定的滤波器指标进行滤波器设计；
3、掌握数字巴特沃斯滤波器的设计原理和步骤；
4、加深对冲激响应不变法和双线性变换法设计IIR数字滤波器的了解，掌握MATLAB函数实现冲激响应变换的方法。
设计IIR滤波器时，首先根据模拟滤波器的指标设计出相应的模拟滤波器Ha（s），然后将设计好的模拟滤波器Ha（s）转换成满足给定技术指标的数字滤波器H（z）。
1、模拟原型低通滤波器Ha（s）的设计
模拟低通滤波器的幅频特性如图2-1所示。
图2-1模拟低通滤波器的幅频特性
模拟低通滤波器的技术指标有通带截止频率
、通带内的最大衰减系数
、阻带截止频率
和阻带内的最小衰减系数
。
，
已知
、
和
，可由止式求出滤波器的阶数N。
求出的N可能有小数部分一般取大于等于N的最小整数。
模拟滤波器的技术指标给定后，需要根据这组指标设计模拟系统函数
，使其逼近理想滤波器特性，一般是根据幅度平方函数
来逼近的。
典型的模拟滤波器有巴特沃斯（Butterworth）滤波器、切比雪夫（Chebyshev）滤波器。
2、实验用MATLAB语言工具函数简介
在MATLAB的数字信号处理工具箱中提供了相应的设计函数，常用的有：
1）Butterworth模拟/数字滤波器设计
调用格式1：
[N,Wn]=buttord（Wp,Ws,Rp,Rs，’s’）
输入参数：
Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减；
输出参数：
N符合要求的滤波器最小阶数，Wn为Butterworth滤波器固有频率（3dB）。
调用格式2：
[b,a]=butter（n,wn,'
ftype'
’s’）
[b,a]=butter（n,wn,'
）
说明：
N和Wn分别为滤波器的阶数和3dB截止频率。
利用此函数可以获得滤波器系统函数的分子多项式（b）和分母多项式（a）的系数
选项中加入‘S’用于设计各种模拟Butterworth滤波器；
不加设计各种数字Butterworth滤波器
Ftype为缺省，设计低通滤波器
Ftype＝hign，设计高通滤波器
Ftype＝stop，设计带阻滤波器
调用格式3：
[z,p,k]=buttap（N）:
设计一个N阶的归一化的巴特沃斯原型低通模拟滤波器，返回滤波器的零点、极点和增益，此时z为空。
例2-1设计一个5阶Butterworth数字高通滤波器，阻带截止频率为250Hz。
设采样频率为1000Hz。
解：
源代码如下：
[b,a]=butter（5,250/500,'
high'
[z,p,k]=butter（5,250/500,'
freqz（b,a,512,1000）
例2-2：
设计一个10阶的带通巴特沃斯数字滤波器，带通频率为100Hz到200Hz，采样频率为1000Hz，绘出该滤波器的幅频于相频特性，以及其冲击响应图
clearall;
N=10;
Wn=[100200]/500;
[b,a]=butter（N,Wn,’bandpass’）;
freqz（b,a,128,1000）
figure
（2）
[y,t]=impz（b,a,101）;
stem（t,y）
例2-3设计一个巴特沃斯模拟低通滤波器，要求通带截止频率
，通带最大衰减
，阻带截止频率
，阻带最小衰减
要求绘出滤波器的幅频特性曲线。
参考程序如下：
fp=5000;
wp=2*pi*fp;
fs=12000;
ws=2*pi*fs;
rp=2;
rs=30;
[N,Wn]=buttord（wp,ws,rp,rs,'
s'
）;
[b,a]=butter（N,Wn,'
freqs（b,a）
运行结果如图2-2所示。
图2-2巴特沃斯模拟低通滤波器的幅频特性和相频特性
如果将freqs（b,a）改为：
H=20*log10（abs（freqs（b,a,w）））;
运行结果如图2-3所示。
图2-3巴特沃斯模拟低通滤波器幅频特性的分贝表示
由图2-3可知该设计结果满足在通带最大衰减
的要求。
2）chebyshevI、chebyshevII型模拟/数字滤波器设计
调用格式：
[b,a]=cheby1（n,Rp,wn,'
）
[b,a]=cheby2（n,Rs,wn,'
例2-4设计一个7阶chebyshevII型数字低通滤波器，截止频率为3000Hz，Rs＝30dB。
源程序如下：
[b,a]=cheby2（7,30,300/500'
）；
freqz（b,a,512,1000）
3、模拟域的频率变换法
1、lp2lp低通到低通模拟滤波器变换。
[bt,at]=lp2lp（b,a,w0），将系统函数表示的截止频率为1rad/s的模拟低通滤波器原型变换为截止频率为w0的低通滤波器。
2、lp2hp低通到高通模拟滤波器变换。
[bt,at]=lp2hp（b,a,w0），将系统函数表示的截止频率为1rad/s的模拟低通滤波器原型变换为截止频率为w0的高通滤波器。
3、lp2bp低通到带通模拟滤波器变换。
[bt,at]=lp2bp（b,a,w0,bw），将系统函数表示的截止频率为1rad/s的模拟低通滤波器原型变换为中心频率为w0、带宽为bw的带通滤波器。
如果被设计的滤波器低端截止频率为w1，高端截止频率为w2，则
W0=sqrt（w1*w2），bw=w2-w1
4、lp2bs低通到带阻模拟滤波器变换。
[bt,at]=lp2bs（b,a,w0,bw），将系统函数表示的截止频率为1rad/s的模拟低通滤波器原型变换为中心频率为w0、带宽为bw的带阻滤波器。
W0=sqrt（w1*w2），bw=w2-w1。
4、脉冲响应不变法和双线性变换法
1）冲激响应不变法
一般来说，在要求时域冲激响应能模仿模拟滤波器的场合，一般使用冲激响应不变法。
冲激响应不变法一个重要特点是频率坐标的变换是线性的，因此如果模拟滤波器的频响带限于折叠频率的话，则通过变换后滤波器的频响应可不失真的反映原响应与频率的关系。
例2-5设计一个中心频率为500HZ，带宽为600Hz的数字带通滤波器，采样频率为1000Hz。
[z,p,k]=buttap（3）;
[b,a]=zp2tf（z,p,k）;
[bt,at]=lp2bp（b,a,500*2*pi,600*2*pi）;
[bz,az]=impinvar（bt,at,1000）;
%将模拟滤波器变换成数字滤波器
freqz（bz,az,512,'
whole'
1000）
2）双线性变换法
双线性变换——将s域映射到z域。
与冲激响应不变法比较，双线性变换的主要优点是靠频率的非线性关系得到S平面与Z平面的单值一一对应关系，整个值对应于单位圆一周。
所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。
[numd,dend]=bilinear（num,den,Fs），将模拟域系统函数转换为数字域的系统函数，Fs为采样频率。
例2-6基于Butterworth模拟滤波器原型，用双线性转换设计数字滤波器，其中参数指标为：
通带截止频率：
通带波动值：
阻带截止频率：
阻带波动值：
首先确定滤波器的阶数N，同时根据
确定
＝0.5。
接着使用bilinear进行双线性转换，最后绘制在频域上的各种图像，其源代码如下：
％数字滤波器指标
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
％将数字滤波器指标反转变化为模拟滤波器的参数
T=1;
fs=1/T;
omegap=（2/T）*tan（wp/T）;
omegas=（2/T）*tan（ws/T）;
ep=sqrt（10^（Rp/10）-1）;
Ripple=sqrt（1/（1+ep*ep））;
Attn=1/（10^（As/20））;
％butterworth原型模拟滤波器的设计
[cs,ds]=afd_butt（omegap,omegas,Rp,As）;
％双线性变换
[b,a]=bilinear（cs,ds,T,fs）;
％频域图像的绘制
freqz（b,a）;
4、实验内容及步骤
1、采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求：
wp=0.25*pi，rp=1db，ws=0.4*pi，as=15db，滤波器采样频率Fs=100hz。
2、利用MATLAB编程设计一个Butterworth数字带通滤波器，指标要求如下：
通带边缘频率：
ωp1=0.45rad，ωp2=0.65rad，通带峰值起伏:
αp≤1dB;
阻带边缘频率：
ωs1=0.3rad，ωs2=0.8rad，最小阻带衰减：
αs≥40dB。
分别用脉冲响应不变法和双线性变换法进行IIR数字滤波器的设计。
Ts=0.125ms
1.双线性变换法和冲激响应不变法相比较，有哪些优点和缺点？
实验三FIR数字滤波器的设计
1、进一步理解FIR滤波器的线性相位特性，熟悉FIR滤波器的幅频特性、相频特性和零点分布情况。
2、掌握用窗函数法设计FIR数字滤波器的原理及方法，了解各种窗函数对滤波器性能的影响。
1、标准型FIR数字滤波器的MATLAB实现
在MATLAB的数字信号处理工具箱中提供了函数fir1。
fir1是采用经典窗函数法设计线性相位FIR数字滤波器，且具有标准低通、带通、高通和带阻等类型。
语法格式：
B=fir1（n,Wn）
B=fir1（n,Wn,’ftype’）
B=fir1（n,Wn,window）
B=fir1（n,Wn,’ftype’,window）
其中：
n为FIR滤波器的阶数，对于高通、带阻滤波器n取偶数。
Wn是截止频率，其取值在0～1之间，对于带通、带阻滤波器，Wn=[W1,W2],且W1<
W2。
‘ftype’为滤波器类型，缺省时为低通或带通滤波器，为’high’时是高通滤波器,为’stop’为带阻滤波器。
window为窗函数，列向量，其长度为n+1，缺省时，自动取hamming窗。
输出参数B为FIR滤波器系数向量h（n），即单位响应，长度为n+1。
特别强调
的长度与滤波器阶数间的关系：
FIR滤波器的系统函数可表示为：
，
的长度为
，而滤波器的阶数为
阶。
2、用窗函数法设计FIR数字滤波器
常用的窗函数主要用于数字FIR滤波器的设计中，可根据实际情况选择合适的窗函数。
注意，在以下的公式中，由于MATLAB的矢量下标规定从1而不是从0开始，因此这些表达式可能与前面章节中的有关信号的起点从0开始稍有不同。
Matlab中提供了用于窗函数法设计FIR数字滤波器的函数，其调用格式如下：
1）矩形窗——w=boxcar（N）：
产生一长度为N的矩形窗。
时域表达式：
2）三角窗——w=triang（N）：
产生一长度为N的三角窗。
当n为奇数时，
当n为偶数时，
3）巴特利特窗——w=bartlett（N）：
产生一长度为N的巴特利特窗。
（3.10）
（3.11）
巴特利特窗与三角窗很相似，但巴特利特窗在第1和第N各采样点数上都是0，而三角窗不是。
当n为奇数时，triang（N-2）等于bartlett（N）。
4）汉宁窗（升余弦窗）——w=hanning（N）：
产生一长度为N的汉宁窗（升余弦窗）。
汉宁窗时域表达式：
（3.12）
w=hamming（N）：
产生一长度为N的哈明窗（改进的升余弦窗）。
哈明窗时域表达式：
（3.13）
下面举例说明如何运用上面的函数进行FIR数字滤波器的设计。
例3-1分别用矩形窗和哈明窗设计FIR低通滤波器，设窗宽
，截止频率
，要求绘出两种窗函数设计的滤波器幅频曲线。
N=11;
h1=fir1（N-1,0.2,boxcar（N））;
h2=fir1（N-1,0.2,hamming（N））;
w=0:
0.01:
pi;
H1=freqz（h1,1,w）;
H1db=20*log10（abs（H1）/max（abs（H1）））;
H2=freqz（h2,1,w）;
H2db=20*log10（abs（H2）/max（abs（H1）））;
subplot（1,2,1）
plot（w/pi,abs（H1）,'
-.'
w/pi,abs（H2））
legend（'
矩形窗'
'
哈明窗'
xlabel（'
w/pi'
ylabel（'
幅频响应'
subplot（1,2,2）
pl
升级会员 