哈工大数字信号处理实验报告

H a r b i n I n s t i t u t e o f T e c h n o l o g y 实验 课程名称: 数字信号处理 实验题目:用FFT 作谱分析、用窗函数法设计数字滤波器 院 系: 电子与信息工程学院 班 级: 哈尔滨工业大学 一、实验目的 (1) 进一步加深DFT 算法原理和基本性质的理解。 (2) 熟悉FFT 算法原理和FFT 子程序的应用。 (3) 学习用FFT 对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析误差及其原因,以便在实际中正确应用FFT 。 二、实验步骤 (1) 复习DFT 的定义、性质和用DFT 作谱分析的有关内容。 (2) 复习FFT 算法原理与编程思想,并对照DIT-FFT 运算流图和程序框图。 (3) 编制信号产生程序,产生以下典型信号供谱分析用: 1423()() 1,03()847 403()347 0x n R n n n x n n n n n x n n n =?+≤≤? =-≤≤???-≤≤?? =-≤≤??? 456()cos 4 ()sin 8 ()cos8cos16cos 20x n n x n n x t t t t π π πππ=≤≤=≤≤=++,0n 19,0n 19 (4) 按实验内容要求,上机实验,并写出实验报告。 三、实验原理 1、FFT 产生背景。 离散傅立叶变换从理论上解决了傅立叶变换应用于实际的可能性,但若直接按DFT 公式计算,运算量太大(与N 2成比例)。快速傅立叶变换(FFT )是离散傅立叶变换的快速算法,它大大减少了离散傅立叶变换的运算次量,一般可缩短一、二故数量级,且N 越大改善越明显,从而使DFT 的运算在实际工作中才真正得到广泛应用。 2、FFT 原理。 减少离散傅立叶变换运算次数的方法基于W nk N 的对称性和周期性。 (1)对称性 W W kn N kn N -==*n)-k(N N )(W (2)周期性 W W W nk N k N n N kn N N -+===)()(mod kn N W )( 由此可得 W W W W W k N N k N N N nk N k n N N -=-===+--) 2/(2/)(k)-n(N N 1W 这样可以将DFT 运算中有些项进行合并,利用W nk N 的对称性和周期性,将长序列的DFT 分解为短序列的DFT 。 3、FFT 运算量。 在采用按时间抽取法的情况下,当N=2M 时共有M 级蝶形,每级都有N/2个蝶形运算。计算每个蝶形需要一次复数乘法和两次复数加法,因而每级运算需N/2次复数乘法和N 次复数加法,这样M 级运算总共需要: 复数乘法次数: m F =N/2.M=N/2.log 2N 复数加法次数: a F =NM=Nlog 2N 在一般情况下,复数乘法所需时间比复数加法多,因此以复数乘法为例将DFT 运算量与FFT 运算量进行对比。直接进行DFT 运算复数乘法次数为N 2次,利用FFT 运算复数乘法为N/2.log 2N 次。 计算量之比为: N N N log log N 222 22N FFT DFT ==复数乘法次数复数乘法次数 四、实验结果、源程序及误差分析 1、 实验程序: N1=8; N2=16; x1=[1,1,1,1]; y11=fft(x1,8); y12=fft(x1,16) subplot(2,2,1); stem(0:3,x1); title('函数x1时域图像');box off subplot(2,2,2); stem(0:(N1-1),abs(y11)); title('N=8,x1的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y12)); title('N=16,x1的DFT图像');box off 实验误差分析:理论的DFT函数应为抽样函数,由实验结果可以看出,随着FFT点数的 增加,FFT图像包络函数越趋近于抽样函数。 2、 实验结果: 实验程序: N1=8; N2=16; xa=1:1:4; xb=4:-1:1; x2=[xa,xb]; y21=fft(x2,8); y22=fft(x2,16); subplot(2,2,1); stem(0:7,x2); title('函数x2时域图像'); subplot(2,2,2);box off stem(0:(N1-1),abs(y21)); title('N=8,x2的DFT图像');box off subplot(2,2,4); stem(0:(N2-1),abs(y22)); title('N=16,x2的DFT图像');box off 实验误差分析:理论的DFT函数应为形如抽样函数的平方的函数,由实验结果可以看出,随着FFT点数的增加,FFT图像包络函数接近抽样函数的平方。 3、 实验结果: 实验程序: N1=8; N2=16; xa=4:-1:1; xb=1:1:4; x3=[xa,xb]; y31=fft(x3,8); y32=fft(x3,16); subplot(2,2,1); stem(0:7,x3); title('函数x3时域图像'); subplot(2,2,2);box off stem(0:(N1-1),abs(y31)); title('N=8时,x3的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y32)); title('N=16时,x3的DFT图像');box off 实验误差分析:随着N的增加,曲线接近理论值,具体的分析和2中的信号相同 4、 实验结果: 实验程序: N1=8; N2=16; n=0:1:19 x4=cos(0.25*pi*n); n=0:1:7; x411=cos(0.25*pi*n); n=8:1:15; x412=cos(0.25*pi*n); n=16:1:19; x413=cos(0.25*pi*n); y411=fft(x411,8); y412=fft(x412,8); y413=fft(x413,8); y41=y411+y412+y413; n=0:1:15; x421=cos(0.25*pi*n); n=16:1:19; x422=cos(0.25*pi*n); y421=fft(x421,16); y422=fft(x422,16); y42=y421+y422; subplot(2,2,1); stem(0:19,x4); title('函数x4时域图像');box off subplot(2,2,2); stem(0:(N1-1),abs(y41)); title('N=8,x4的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y42)); title('N=16,x4的DFT图像');box off 实验误差分析:理论的DFT函数为冲激函数,从实验结果来看,有两个明显的峰值可以看到,与理论值接近。 5、 实验结果: 实验程序: n=0:1:19; N1=8; N2=16; x5=sin(0.125*pi*n); n=0:1:7; x511=sin(0.125*pi*n); n=8:1:15; x512=sin(0.125*pi*n); n=16:1:19; x513=sin(0.125*pi*n); y511=fft(x511,8); y512=fft(x512,8); y513=fft(x513,8); y51=y511+y512+y513; n=0:1:15; x521=sin(0.125*pi*n); n=16:1:19; x522=sin(0.125*pi*n); y521=fft(x521,16); y522=fft(x522,16); y52=y521+y522; subplot(2,2,1); stem(0:19,x5); title('函数x5时域图像');box off subplot(2,2,2); stem(0:(N1-1),abs(y51)); title('N=8,x5的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y52)); title('N=16,x5的DFT图像');box off 实验误差分析:理论FFT抽样频谱为单位冲激函数,因此只要满足抽样定理,FFT都可以与理论图形一样的曲线。当N=16时,不满足抽样定理,因此产生了混叠失真,当N=8时,满足抽样定理,因此可得到与理论曲线一致的图形。 6、x4(n)+x5(n) 实验结果: 实验程序: n=0:1:19; N1=8; N2=16; x4=cos(0.25*pi*n); x5=sin(0.125*pi*n); x45=x4+x5; n=0:1:7; x4511=cos(0.25*pi*n)+sin(0.125*pi*n); n=8:1:15; x4512=cos(0.25*pi*n)+sin(0.125*pi*n); n=16:1:19; x4513=cos(0.25*pi*n)+sin(0.125*pi*n); y4511=fft(x4511,8); y4512=fft(x4512,8); y4513=fft(x4513,8); y451=y4511+y4512+y4513; n=0:1:15; x4521=cos(0.25*pi*n)+sin(0.125*pi*n); n=16:1:19; x4522=cos(0.25*pi*n)+sin(0.125*pi*n); y4521=fft(x4521,16); y4522=fft(x4522,16); y452=y4521+y4522; subplot(2,2,1); stem(0:19,x45); title('函数x45时域图像');box off subplot(2,2,2); stem(0:(N1-1),abs(y451)); title('N=8,x45的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y452)); title('N=16,x45的DFT图像');box off 实验误差分析:这个函数的理论FFT抽样频谱为有两个不同频率分量的单位冲激函数,因此只要满足抽样定理,FFT都可以与理论图形一样的曲线。当N=16时,不满足抽样定理,因此产生了混叠失真。 7、x4(n)+j*x5(n) 实验结果: 实验程序: n=0:1:19; N1=8; N2=16; x4=cos(0.25*pi*n); x5=sin(0.125*pi*n); x45=x4+x5; n=0:1:7; x4511=cos(0.25*pi*n)+j*sin(0.125*pi*n); n=8:1:15; x4512=cos(0.25*pi*n)+j*sin(0.125*pi*n); n=16:1:19; x4513=cos(0.25*pi*n)+j*sin(0.125*pi*n); y4511=fft(x4511,8); y4512=fft(x4512,8); y4513=fft(x4513,8); y451=y4511+y4512+y4513; n=0:1:15; x4521=cos(0.25*pi*n)+j*sin(0.125*pi*n); n=16:1:19; x4522=cos(0.25*pi*n)+j*sin(0.125*pi*n); y4521=fft(x4521,16); y4522=fft(x4522,16); y452=y4521+y4522; subplot(2,2,1); stem(0:19,x45); title('函数x45时域图像');box off subplot(2,2,2); stem(0:(N1-1),abs(y451)); title('N=8,x45的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y452)); title('N=16,x45的DFT图像');box off 实验误差分析:理论FFT抽样频谱为两个单位冲激函数,因此只要满足抽样定理,FFT 都可以与理论图形一样的曲线。当N=8时,不满足抽样定理,因此产生了混叠失真,当N=16时,满足抽样定理,因此可得到与理论曲线一致的图形。 8、 实验结果: 实验程序: n=0:1:69; N1=16; N2=32; N3=64; fs=64; x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=0:1:15; x611=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=16:1:31; x612=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=32:1:47; x613=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=48:1:63; x614=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=64:1:69; x615=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); y611=fft(x611,N1); y612=fft(x612,N1); y613=fft(x613,N1); y614=fft(x614,N1); y615=fft(x615,N1); y61=y611+y612+y613+y614+y615; n=0:1:31; x621=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=32:1:63; x622=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=64:1:69; x623=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); y621=fft(x621,N2); y622=fft(x622,N2); y623=fft(x623,N2); y62=y621+y622; n=0:1:63; x631=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); n=64:1:69; x632=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); y631=fft(x631,N3); y632=fft(x632,N3); y63=y631+y632; subplot(2,2,1); stem(0:69,x6); title('函数x6时域图像');box off subplot(2,2,2); stem(0:(N1-1),abs(y61)); title('N=16,x6的DFT图像');box off subplot(2,2,3); stem(0:(N2-1),abs(y62)); title('N=32,x6的DFT图像');box off subplot(2,2,4); stem(0:(N3-1),abs(y63)); title('N=64,x6的DFT图像');box off 实验误差分析:理论FFT抽样频谱为两个单位冲激函数,因此只要满足抽样定理,FFT 都可以与理论图形一样的曲线。当N=16时,不满足抽样定理,因此产生了混叠失真,当N=32,64时,满足抽样定理,因此可得到与理论曲线一致的图形。 五、思考题 (1) 在N=8时,x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢? 答:从实验结果看出,N=8时,x2(n)和x3(n)的幅频特性相同。而N=16时,两者不同。原因是N=8时,x3(n)只是以长度8为周期,将其延拓称周期序列然后加以移位,最后截取主值区间的序列值。因此x3(n)是x2(n)进行循环位移后的结果。移位后只影响DFT相频特性而不影响幅频特性。而当N=16时,两个序列不满足循环位移的性质,因此幅频特性发生变化。 (2) 如果周期信号的周期预先不知道,如何用FFT进行谱分析? 答:若该信号是模拟信号,用FFT进行频谱分析时,首先对信号进行采样,使之变成离散信号,然后用FFT对离散信号进行频谱分析。由采样定理,采样频率fs应当大于信号最高频率的2倍。如果信号的周期预先不知道,可以选择一个较大的采样频率,用试验法逐步增大采样频率,观察频谱特性。 一、实验目的 (1)熟悉矩形窗、汉宁窗、海明窗和布莱克曼窗。 (2) 掌握用上述窗函数法设计FIR 数字滤波器的原理和方法。 (3) 熟悉线性相位FIR 数字滤波器特性。 (4) 了解各种窗函数对滤波特性的影响。 二、实验原理与方法 如果所希望的滤波器的理想频率响应函数为Hd(e j ω),则其对应的单位脉冲响应为 1 ()()2j j n d d h n H e e d π ωωπ ω π - = ? 用窗函数w(n)将hd(n)截断, 并进行加权处理, 得到: ()()()d h n h n n ω= h(n)作为实际设计的FIR 数字滤波器的单位脉冲响应序列, 其频率响应函数H(ej ω)为 1 ()()N j j n n H e h n e ωω--==∑ 如果要求线性相位特性, 则h(n)还必须满足: ()(1)h n h N n =±-- 根据上式中的正、 负号和长度N 的奇偶性又将线性相位FIR 滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如, 要设计线性相位低通特性, 可选择h(n)=h(N-1-n)一类, 而不能选h(n)=-h(N-1-n)一类。 三、实验内容与步骤 (1)复习用窗函数法设计FIR 数字滤波器一节内容,阅读本实验原理,掌握设计步骤。 (2)编写程序,其中幅度特性要求用dB 示。 备注:不许用 freqz() 这个函数 设 ()[()]()()() ()k I H k DFT h n H k H k jH k H k ==+= 画图时,用20log|H(k)|打印幅度特性。第k 点对应的频率为 k N w k π 2= 为使曲线包络更接近H(ejw)的幅度特性曲线,DFT 变换区间要选大些。例如窗口长度N=33 时,可通过在h(n)末尾补零的方法,使长度变为64,再进行64点DFT ,则可得到更精确的幅度衰减特性曲线。 四、上机实验内容 1、用升余弦窗设计一线性相位低通FIR 数字滤波器,截止频率Wc= 4 Π rad 。窗口N=15,33。要求在两种窗口长度情况下,分别求出h(n),打印出相应的幅频特性和相频曲线。观察3db 带宽和20db 带宽,窗口N 对滤波特性的影响。 设计低通FIR 数字滤波器时,一般以理想低通滤波特性为逼近函数 即 ,()0 1211()()22sin ()() c c j a c j d c j j n j a j n d d c e H e N a h n H e e e e d n a n a ωω πωωωωωπω ωωωωπω ππ ωπ----- ?≤?=? <