Harbin Institute of Technology
课程名称：实验题目：利用院 系：
班 级：
姓 名：
学 号：
指导教师：
实验时间：实验报告 数字信号处理I FFT 对信号进行频谱分析窗函数法设计FIR 数字滤波器 电子与信息工程学院 冀振元 2016年11月
哈尔滨工业大学
实验一 利用FFT 对信号进行频谱分析
1. 实验目的
(1)进一步加深DFT 算法原理和基本性质的理解(因为FFT 只是DFT 的一种快速算法，所以FFT 的运算结果必然满足DFT 的基本性质) 。
(2)熟悉FFT 算法原理和FFT 子程序的应用。
(3)学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT 。
2. 实验步骤
(1)复习DFT 的定义、性质和用DFT 作谱分析的有关内容。
(2)复习FFT 算法原理与编程思想，并对照DIT-FFT 运算流图和程序框图。
(3)编制信号产生子程序，产生以下典型信号供谱分析用：
x 1(n ) =R 4(n )
⎧n +1, 0≤n ≤3⎪x 2(n ) =⎨8-n 4≤n ≤7
⎪⎩0
⎧4-n 0≤n ≤3⎪x 3(n ) =⎨n -34≤n ≤7
⎪0⎩
x 4(n ) =cos
x 5(n ) =sin π4n n ，0≤n ≤19，0≤n ≤19π
8
x 6(t ) =cos8πt +cos16πt +cos 20πt
(4)编写主程序。
(5)按实验内容要求，上机实验，并写出实验报告。
3. 实验原理
用FFT 对信号作频分析是学习数字信号处理的重要内容，经常需要进行分析的信号是模拟信号的时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D 和分析误差。频谱分辨率直接和FFT 的变换区间N 有关，因为FFT 能够实现的频率分辨率是2π/N，因此要求2π/N小于等于D 。可以根据此式选择FFT 的变换区间N 。误差主要来自于用FFT 作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N 较大时，离散谱的包络才能逼近连续谱，因此N 要适当选择大一些。
4. 实验结果
(1) x1序列及其8点、16点FFT 谱分析
由图可知，采样频率越高，得到的频谱的分辨率就越高，当N 趋近于无穷大时频谱包络接近与理论函数。
(2) x2序列及其8点、16点FFT 谱分析
同样地，由图可知，采样频率越高，得到的频谱的分辨率就越高，当N 趋近于无穷大时频谱包络接近与理论函数。
(3) x3序列及其8点、16点FFT 谱分析
注意到，x3的8点FFT 与x2的8点FFT 的幅度谱相同，而二者的16点FFT 幅度谱则不同。
(4) x4序列及其8点、16点FFT 谱分析
由图可知，当N=8和N=16时均满足采样定理，故其FFT 可以得到与理论相符的频谱图。
(5) x5序列及其8点、16点FFT 谱分析
与x4不同，x5的FFT 在N=8时不满足抽样定理，因而产生了混叠；N=16时，满足抽样定理，可得到与理论曲线一致的图形。
(6) x6序列及其16点、32点、64点FFT 谱分析
由图可知，x6的FFT 在N=16时不满足抽样定理，因而产生了混叠；而在N=32,64时，满足抽样定理，得到了与理论曲线一致的图形。
(7) x7=x4+x5序列及其8点、16点FFT 谱分析
易知x7序列的理论FFT 抽样频谱应有两个不同频率分量的单位冲激函数，故只要满足抽样定理，便可以得到与理论图形一样的曲线。由图形可知，当N=8时，不满足抽样定理，产生了混叠；当N=16时，满足抽样定理，得到了与理论曲线一致的图形。
(8) x8=x4+jx5序列及其8点、16点FFT 谱分析
同样地，由图形可知，当N=8时，不满足抽样定理，产生了混叠；当N=16时，满足抽样定理，得到了与理论曲线一致的图形。
5. 误差分析
实验中得到的频谱图与理论基本相符。N 较小时，有可能会不满足采样定理，从而造成频谱混叠失真，产生误差；而当N 较大时，可以比较完整地反映原信号的真实频谱，误差相对较小。
6. 参数选取
由上述实验结果可知，N 要大于信号的长度M ，这样才能由抽样序列不失真地恢复出原信号，故尽可能地增大N 可以有效地反映原信号的信息。
7. 思考题
(1)在N =8时，x 2(n ) 和x 3(n ) 的幅频特性会相同吗? 为什么? N =16呢? N=8时二者的幅频特性相同，N=16时二者的幅频特性不同。 原因：x3可由x2移位得到，而二者不为0的区间长度恰为8，故在N=8时，经周期延拓、移位、截取主值区间，二者的幅频特性相同；而在N=16时，经过了补零，二者便有了区别，幅频特性也就不同了。还注意到，二者的相频特性都不同，这是因为循环位移会影响其傅里叶变换后的相频特性，故不同。
(2)如果周期信号的周期预先不知道， 如何用FFT 进行谱分析? 应先适当截取M 点进行FFT ，再将截取的长度扩大1倍重新截取，比较二者结果，若二者的差别满足分析误差要求，就可以近似得到该信号的频谱，如不满足误差要求则继续将截取的长度加倍进行FFT ，直到结果满足误差要求为止。
8. 实验结论
对离散序列作DFT ，采样时要正确确定采样周期，当周期N 大于离散序列的长度M ，即N ≥M 时，才能得到正确的频谱。当N ＜M 时，
会发生频谱混叠现象，产生失真。
FFT 是DFT 的一种快速算法，其运算次数可大幅度减少，N 越大时减少越明显，因此在实际中常利用FFT 进行计算。
附：实验一全部matlab 代码
%% x1序列的产生与谱分析 clc;clear;close all ;
x1=[1 1 1 1 zeros(1,4)]; %x1序列的产生 subplot(3,2,1:2); stem(0:7,x1); grid on ;
title('x1序列' );
xlabel('n' ); x1_8=fftshift(fft(x1,8)); hx1_8=abs(x1_8);
ax1_8=angle(x1_8); subplot(3,2,3); stem(-4:3,hx1_8); grid on ;
title('N=8时x1的幅度谱' ); xlabel('K' );
ylabel('|Hw|'); subplot(3,2,4); stem(-4:3,ax1_8); grid on ;
title('N=8时x1的相位谱' );
xlabel('K' );ylabel(' 相角' ); x1_16=fftshift(fft(x1,16)); hx1_16=abs(x1_16);
ax1_16=angle(x1_16); subplot(3,2,5); stem(-8:7,hx1_16); grid on ;
title('N=16时x1的幅度谱' );
xlabel('K' );ylabel('|Hw|'); subplot(3,2,6); stem(-8:7,ax1_16); grid on ;
title('N=16时x1的相位谱' );
xlabel('K' );ylabel(' 相角' );
%% x2序列的产生与谱分析 clc;clear;close all ;
x2=[1 2 3 4 4 3 2 1]; subplot(3,2,1:2); stem(0:7,x2);
%画出x1序列 %将x1进行8点的FFT 变换并得到其幅度和相角%画出x1的8点FFT 的幅度谱 %画出x1的16点FFT 的相位谱 %将x1进行16点的FFT 变换并得到其幅度和相角 %画出x1的16点FFT 的幅度谱 %画出x1的16点FFT 的相位谱 %x2序列的产生
title('x2序列' );
xlabel('n' ); %画出x2序列 x2_8=fftshift(fft(x2,8)); hx2_8=abs(x2_8);
ax2_8=angle(x2_8); %将x2进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3); stem(-4:3,hx2_8); grid on ;
title('N=8时x2的幅度谱' );
xlabel('K' );ylabel('|Hw|'); subplot(3,2,4); stem(-4:3,ax2_8); grid on ;
title('N=8时x2的相位谱' );
xlabel('K' );ylabel(' 相角' ); x2_16=fftshift(fft(x2,16)); hx2_16=abs(x2_16);
ax2_16=angle(x2_16); subplot(3,2,5); stem(-8:7,hx2_16); grid on ;
title('N=16时x2的幅度谱' );
xlabel('K' );ylabel('|Hw|'); subplot(3,2,6); stem(-8:7,ax2_16); grid on ;
title('N=16时x2的相位谱' );
xlabel('K' );ylabel(' 相角' );
%% x3序列的产生与谱分析 clc;clear;close all ;
x3=[4 3 2 1 1 2 3 4]; subplot(3,2,1:2); stem(0:7,x3); grid on ;
title('x3序列' );
xlabel('n' ); x3_8=fftshift(fft(x3,8)); hx3_8=abs(x3_8);
ax3_8=angle(x3_8); subplot(3,2,3); stem(-4:3,hx3_8);
%画出x2的8点FFT 的幅度谱 %画出x2的16点FFT 的相位谱 %将x2进行16点的FFT 变换并得到其幅度和相角%画出x2的16点FFT 的幅度谱 %画出x2的16点FFT 的相位谱 %x3序列的产生 %画出x3序列 %将x3进行8点的FFT 变换并得到其幅度和相角
title('N=8时x3的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x3的8点FFT 的幅度谱 subplot(3,2,4); stem(-4:3,ax3_8); grid on ;
title('N=8时x3的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x3的16点FFT 的相位谱 x3_16=fftshift(fft(x3,16)); hx3_16=abs(x3_16);
ax3_16=angle(x3_16); %将x3进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5); stem(-8:7,hx3_16);
grid on ;title('N=16时x3的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x3的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax3_16); grid on ;title('N=16时x3的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x3的16点FFT 的相位谱
%% x4序列的产生与谱分析 clc;clear;close all ;
x4=cos(pi/4*(0:1:19)); %x4序列的产生 subplot(3,2,1:2);stem(0:19,x4); grid on ;title('x4序列' );
xlabel('n' ); %画出x4序列 x4_8=fftshift(fft(x4,8));
hx4_8=abs(x4_8); ax4_8=angle(x4_8); %将x4进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx4_8); grid on ;title('N=8时x4的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x4的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax4_8); grid on ;title('N=8时x4的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x4的16点FFT 的相位谱 x4_16=fftshift(fft(x4,16));
hx4_16=abs(x4_16); ax4_16=angle(x4_16); %将x4进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx4_16); grid on ;title('N=16时x4的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x4的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax4_16); grid on ;title('N=16时x4的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x4的16点FFT 的相位谱
%% x5序列的产生与谱分析 clc;clear;close all ;
x5=sin(pi/8*(0:1:19)); %x5序列的产生 subplot(3,2,1:2);stem(0:19,x5); grid on ;title('x5序列' );
xlabel('n' ); %画出x5序列 x5_8=fftshift(fft(x5,8));
hx5_8=abs(x5_8); ax5_8=angle(x5_8); %将x5进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx5_8); grid on ;title('N=8时x5的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x5的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax5_8); grid on ;title('N=8时x5的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x5的16点FFT 的相位谱 x5_16=fftshift(fft(x5,16));
hx5_16=abs(x5_16); ax5_16=angle(x5_16); %将x5进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx5_16); grid on ;title('N=16时x5的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x5的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax5_16); grid on ;title('N=16时x5的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x5的16点FFT 的相位谱
%% x6序列的产生与谱分析 clc;clear;close all ; fs=64; N=64;
t=(0:1:(N-1))/fs;
x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t); %x6序列的产生 subplot(4,2,1:2);stem(0:63,x6); grid on ;title('x6序列' );
xlabel('n' ); %画出x6序列 x6_16=fftshift(fft(x6,16));
hx6_16=abs(x6_16); ax6_16=angle(x6_16); %将x6进行16点的FFT 变换并得到其幅度和相角 subplot(4,2,3);stem(-8:7,hx6_16); grid on ;title('N=16时x6的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x6的16点FFT 的幅度谱 subplot(4,2,4);stem(-8:7,ax6_16); grid on ;title('N=16时x6的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x6的16点FFT 的相位谱 x6_32=fftshift(fft(x6,32));
hx6_32=abs(x6_32); ax6_32=angle(x6_32); %将x6进行32点的FFT 变换并得到其幅度和相角 subplot(4,2,5);stem(-16:15,hx6_32);
grid on ;title('N=32时x6的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x6的32点FFT 的幅度谱 subplot(4,2,6);stem(-16:15,ax6_32); grid on ;title('N=32时x6的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x6的32点FFT 的相位谱 x6_64=fftshift(fft(x6,64));
hx6_64=abs(x6_64); ax6_64=angle(x6_64); %将x6进行64点的FFT 变换并得到其幅度和相角 subplot(4,2,7);stem(-32:31,hx6_64); grid on ;title('N=64时x6的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x6的64点FFT 的幅度谱 subplot(4,2,8);stem(-32:31,ax6_64); grid on ;title('N=64时x6的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x6的64点FFT 的相位谱
%% x7=x4+x5序列的产生与谱分析 clc;clear;close all ; x4=cos(pi/4*(0:1:19)); x5=sin(pi/8*(0:1:19));
x7=x4+x5; %x7=x4+x5序列的产生 subplot(3,2,1:2);stem(0:19,x7); grid on ;title('x7序列' );
xlabel('n' ); %画出x7序列 x7_8=fftshift(fft(x7,8));
hx7_8=abs(x7_8); ax7_8=angle(x7_8); %将x7进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx7_8); grid on ;title('N=8时x7的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x7的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax7_8); grid on ;title('N=8时x7的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x7的16点FFT 的相位谱 x7_16=fftshift(fft(x7,16));
hx7_16=abs(x7_16); ax7_16=angle(x7_16); %将x7进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx7_16); grid on ;title('N=16时x7的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x7的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax7_16); grid on ;title('N=16时x7的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x7的16点FFT 的相位谱
%% x8=x4+jx5序列的产生与谱分析 clc;clear;close all ; x4=cos(pi/4*(0:1:19));
x5=sin(pi/8*(0:1:19));
x8=x4+1i*x5; %x8=x4+jx5序列的产生 subplot(3,2,1:2);stem(0:19,x8); grid on ;title('x8序列' );
xlabel('n' ); %画出x8序列 x8_8=fftshift(fft(x8,8));
hx8_8=abs(x8_8); ax8_8=angle(x8_8); %将x8进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx8_8); grid on ;title('N=8时x8的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x8的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax8_8); grid on ;title('N=8时x8的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x8的16点FFT 的相位谱 x8_16=fftshift(fft(x8,16));
hx8_16=abs(x8_16); ax8_16=angle(x8_16); %将x8进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx8_16); grid on ;title('N=16时x8的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x8的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax8_16); grid on ;title('N=16时x8的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x8的16点FFT 的相位谱
实验二 窗函数法设计FIR 数字滤波器
1. 实验目的
(1) 熟悉矩形窗、汉宁窗、海明窗等常用窗函数。
(2) 掌握用上述窗函数法设计FIR 数字滤波器的原理和方法。 (3) 熟悉线性相位FIR 数字滤波器特性。 (4) 了解各种窗函数对滤波特性的影响。
2. 实验原理与方法
如果所希望的滤波器的理想频率响应函数为Hd (e jω) ，则其对应的单位冲激响应为
1
h d (n ) =
2π
⎰π
-
π
H d (e j ω) e j ωn d ω
用窗函数ω(n ) 将h d (n ) 截断，得到:
h (n ) =h d (n ) ω(n )
频率响应函数H (e jω) 为
h (n ) 就作为实际设计的FIR 数字滤波器的单位冲激响应序列， 其
H (e j ω) =∑h (n ) e -j ωn
n =0
N -1
如果要求线性相位特性， 则h (n ) 还必须满足：
h (n ) =±h (N -1-n )
根据上式中的正、负号和长度N 取为奇数或偶数又将线性相位FIR 滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。
第一类：h(n)=h(N-1-n)且N 为奇数，此类情况适合低通、高通、带通、带阻；
第二类：h(n)=h(N-1-n)且N 为偶数，此类情况适合设计低通、带通滤波器；
第三类：h(n)=-h(N-1-n)且N 为奇数，此类情况只适合设计带通滤波器；
第四类：h(n)=-h(N-1-n)且N 为偶数，此类情况适合设计高通、带通滤波器。
3. 实验内容及步骤
(1)复习用窗函数法设计FIR 数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。 (2) 编写程序。
①编写能产生三种窗函数的子程序。
②编写主程序。其中幅度特性要求用分贝（dB ）表示。
设
H (k ) =DFT[h (n )]
H (k ) =H k (k ) +jH I (k )
H (k ) =
画图时，用20log|H(k)|画出幅度特性。
(1) 15阶矩形窗、汉宁窗和海明窗
(2)33阶矩形窗、汉宁窗和海明窗
(1)窗函数法特点
窗口法设计的主要优点是简单，使用方便。窗口函数大多有封闭的公式可循，性能、参数都已有表格、资料可供参考， 计算程序简便， 所以很实用，缺点是通带和阻带的截止频率不易控制。 (2)窗口长度N 和窗函数类型对滤波特性的影响
增加窗的长度，可以减少窗的主瓣宽度，从而减少过渡带的带宽，
但增加N 并不能减少带内波动以及加大阻带。这两个指标只能从窗函数的形状上找解决方法。
6. 综合比较
对选用以上三种窗口设计出的滤波器的结果进行比较可以看出，矩形窗的阻带衰减最小，其过渡带宽也最窄；汉宁窗和海明窗衰减比矩形窗大，其过度带宽比矩形窗大；而海明窗又比汉宁窗阻带衰减大，过渡带也大；从矩形窗、汉宁窗到海明窗，阻带衰减逐渐加大，过渡带宽也逐渐加大。在选定窗口的情况下N 值越大，过渡带越窄，所需的计算量和存储量也加大，因此实际设计的时候应根据需要选用不同的窗口。
7. 实验结论
选择窗函数时要注意：
(1)要从保持最大信息和消除旁瓣的综合效果出发来考虑问题，要使窗
函数频谱中的主瓣宽度尽量窄，能量尽可能集中在主瓣内，从而在谱分析时获得较高的频率分辨力，在数字滤波器设计中获得较小的过渡带。
(2)窗谱的旁瓣高度应尽量小而且随频率尽快衰减，以减小谱估计时泄漏失真。在设计数字滤波器时减小通带的波动，提高阻带的衰减。但主瓣既窄，旁瓣又小衰减又快的窗函数是不容易找到的，比如矩形窗旁瓣很大，但其主瓣宽度是最窄的，因此，在数据处理时通常需要做综合考虑取其折中。
(3)在应用窗函数时，除了要考虑窗谱本身的特性外，还应当充分考虑被分析信号的特点以及具体的处理要求。要考虑信号中的信息量的分布，增强信号中所需要的信息部分，压制信号中不需要的信息部分，以感兴趣的有效信息与窗函数作用后的综台效果为最好来选用窗函数，使得处理结果有足够的频谱检测能力和频谱幅值估计精度。
(4)信号的加窗处理，重要的问题是在于根据信号的性质和研究目的来选用窗函数。如对于频率分辨率要求很高，谱估计幅值精度要求不很高的信号，处理时可选用主瓣较窄而便于分辨的矩形窗。
(5)对于要处理的信号，不知选用哪一种时窗函数好时，往往多用几种窗函数进行处理，然后比较用几种窗处理的结果和试验验证的综合考虑来决定选用什么样的窗函数。
8. 思考题
(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何
用窗函数法设计线性相位低通滤波器? 写出设计步骤。
a) 求出对应的数字频率。
通带截止频率ωp=
阻带截止频率ωst=Ωpfs Ωst
fs
b ）求ℎd(n)。设Hd(ejω) 为理想线性相位滤波器
1sin [ωc(n−τ)] n≠τ1π
−jωτjωnℎd=∫eedω={π(n−τ) ω2π−πc n=τπ
其中τ为线性相位所必须的位移，且满足τ=N−12.
c ）求窗函数。由阻带衰减δ2确定窗形状，由过度带宽确定N 。 d ）求h(n)。由窗函数表达式ω(n)确定FIR 滤波器的h(n)。
h (n ) =ℎd(n)∙ω(n)
e ）由h(n)求Hd(ejω) ，检验各项指标是否满足要求。如不满足要求，则要改变N ，或改变窗形状，然后重新计算。
(2) 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应ℎd(n).
①假设理想带通滤波器的频率响应：
-j ωt 0⎧e H d (e j ω) =⎨⎩00
②对H d (e j ω) 做IDTFT 得单位冲激响应h d (n ) ：
1{sin [(n −t0) ω2]−sin [(n−t0)ω1]} n≠t01π(n−t0) jωjωℎd(n) =∫H(e) ∙edω= 2π−πd 1{π(ω2−ω1) n=t0π
③根据指标要求如通带最大衰减、阻带最小衰减可选择窗函数的类型和点数N 。
h (n ) =h d (n ) ⋅ω(n )
④对h (n ) 做DTFT 可得H j ω
d (e ) ，验证指标是否满足设计要求，若不满
足重新选择窗口类型和窗口长度，直至满足设计指标要求。
附：实验二全部matlab 代码
clc;clear;close all ;
N1=15;N2=33; %分别表示15和33阶 wc=pi/4;
a1=(N1-1)/2;a2=(N2-1)/2;
n1=0:N1-1;n2=0:N2-1;
n1=n1';n2=n2';
Hd1=sin(wc*(n1-a1))./(pi*(n1-a1));Hd1(8)=wc/pi;
Hd2=sin(wc*(n2-a2))./(pi*(n2-a2));Hd2(17)=wc/pi;
w1_boxcar=boxcar(N1);
w2_boxcar=boxcar(N2);
w1_hanning=hanning(N1);
w2_hanning=hanning(N2);
w1_hamming=hamming(N1);
w2_hamming=hamming(N2);
h1_boxcar=Hd1.*w1_boxcar; %最终设计的滤波器系数h(n)=Hd(n)*w(n) h1_hanning=Hd1.*w1_hanning;
h1_hamming=Hd1.*w1_hamming;
h2_boxcar=Hd2.*w2_boxcar;
h2_hanning=Hd2.*w2_hanning;
h2_hamming=Hd2.*w2_hamming;
FFT_N=1024; %做FFT 查看幅频和相频曲线(FFT点数应比滤波器阶数要多，使频谱细化) h1_boxcar_fft=fftshift(fft(h1_boxcar,FFT_N));
h1_boxcar_fft_hw=20*log10(abs(h1_boxcar_fft));
h1_boxcar_fft_fai=angle(h1_boxcar_fft);
h1_hanning_fft=fftshift(fft(h1_hanning,FFT_N));
h1_hanning_fft_hw=20*log10(abs(h1_hanning_fft));
h1_hanning_fft_fai=angle(h1_hanning_fft);
h1_hamming_fft=fftshift(fft(h1_hamming,FFT_N));
h1_hamming_fft_hw=20*log10(abs(h1_hamming_fft));
h1_hamming_fft_fai=angle(h1_hamming_fft);
h2_boxcar_fft=fftshift(fft(h2_boxcar,FFT_N));
h2_boxcar_fft_hw=20*log10(abs(h2_boxcar_fft));
h2_boxcar_fft_fai=angle(h2_boxcar_fft);
h2_hanning_fft=fftshift(fft(h2_hanning,FFT_N));
h2_hanning_fft_hw=20*log10(abs(h2_hanning_fft));
h2_hanning_fft_fai=angle(h2_hanning_fft);
h2_hamming_fft=fftshift(fft(h2_hamming,FFT_N));
h2_hamming_fft_hw=20*log10(abs(h2_hamming_fft));
h2_hamming_fft_fai=angle(h2_hamming_fft);
xias_w=(-512:511)/512;
figure; %画出15阶滤波器 subplot(231);plot(xias_w,h1_boxcar_fft_hw);
grid on ;title('15阶矩形窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(234);plot(xias_w,h1_boxcar_fft_fai);
grid on ;title('15阶矩形窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(232);plot(xias_w,h1_hanning_fft_hw);
grid on ;title('15阶汉宁窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(235);plot(xias_w,h1_hanning_fft_fai);
grid on ;title('15阶汉宁窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(233);plot(xias_w,h1_hamming_fft_hw);
grid on ;title('15阶海明窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(236);plot(xias_w,h1_hamming_fft_fai);
grid on ;title('15阶海明窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
figure; subplot(231);plot(xias_w,h2_boxcar_fft_hw);
grid on ;title('33阶矩形窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(234);plot(xias_w,h2_boxcar_fft_fai);
grid on ;title('33阶矩形窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(232);plot(xias_w,h2_hanning_fft_hw);
grid on ;title('33阶汉宁窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(235);plot(xias_w,h2_hanning_fft_fai);
grid on ;title('33阶汉宁窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(233);plot(xias_w,h2_hamming_fft_hw);
grid on ;title('33阶海明窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(236);plot(xias_w,h2_hamming_fft_fai);
grid on ;title('33阶海明窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
%画出33阶滤波器
Harbin Institute of Technology
课程名称：实验题目：利用院 系：
班 级：
姓 名：
学 号：
指导教师：
实验时间：实验报告 数字信号处理I FFT 对信号进行频谱分析窗函数法设计FIR 数字滤波器 电子与信息工程学院 冀振元 2016年11月
哈尔滨工业大学
实验一 利用FFT 对信号进行频谱分析
1. 实验目的
(1)进一步加深DFT 算法原理和基本性质的理解(因为FFT 只是DFT 的一种快速算法，所以FFT 的运算结果必然满足DFT 的基本性质) 。
(2)熟悉FFT 算法原理和FFT 子程序的应用。
(3)学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT 。
2. 实验步骤
(1)复习DFT 的定义、性质和用DFT 作谱分析的有关内容。
(2)复习FFT 算法原理与编程思想，并对照DIT-FFT 运算流图和程序框图。
(3)编制信号产生子程序，产生以下典型信号供谱分析用：
x 1(n ) =R 4(n )
⎧n +1, 0≤n ≤3⎪x 2(n ) =⎨8-n 4≤n ≤7
⎪⎩0
⎧4-n 0≤n ≤3⎪x 3(n ) =⎨n -34≤n ≤7
⎪0⎩
x 4(n ) =cos
x 5(n ) =sin π4n n ，0≤n ≤19，0≤n ≤19π
8
x 6(t ) =cos8πt +cos16πt +cos 20πt
(4)编写主程序。
(5)按实验内容要求，上机实验，并写出实验报告。
3. 实验原理
用FFT 对信号作频分析是学习数字信号处理的重要内容，经常需要进行分析的信号是模拟信号的时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D 和分析误差。频谱分辨率直接和FFT 的变换区间N 有关，因为FFT 能够实现的频率分辨率是2π/N，因此要求2π/N小于等于D 。可以根据此式选择FFT 的变换区间N 。误差主要来自于用FFT 作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N 较大时，离散谱的包络才能逼近连续谱，因此N 要适当选择大一些。
4. 实验结果
(1) x1序列及其8点、16点FFT 谱分析
由图可知，采样频率越高，得到的频谱的分辨率就越高，当N 趋近于无穷大时频谱包络接近与理论函数。
(2) x2序列及其8点、16点FFT 谱分析
同样地，由图可知，采样频率越高，得到的频谱的分辨率就越高，当N 趋近于无穷大时频谱包络接近与理论函数。
(3) x3序列及其8点、16点FFT 谱分析
注意到，x3的8点FFT 与x2的8点FFT 的幅度谱相同，而二者的16点FFT 幅度谱则不同。
(4) x4序列及其8点、16点FFT 谱分析
由图可知，当N=8和N=16时均满足采样定理，故其FFT 可以得到与理论相符的频谱图。
(5) x5序列及其8点、16点FFT 谱分析
与x4不同，x5的FFT 在N=8时不满足抽样定理，因而产生了混叠；N=16时，满足抽样定理，可得到与理论曲线一致的图形。
(6) x6序列及其16点、32点、64点FFT 谱分析
由图可知，x6的FFT 在N=16时不满足抽样定理，因而产生了混叠；而在N=32,64时，满足抽样定理，得到了与理论曲线一致的图形。
(7) x7=x4+x5序列及其8点、16点FFT 谱分析
易知x7序列的理论FFT 抽样频谱应有两个不同频率分量的单位冲激函数，故只要满足抽样定理，便可以得到与理论图形一样的曲线。由图形可知，当N=8时，不满足抽样定理，产生了混叠；当N=16时，满足抽样定理，得到了与理论曲线一致的图形。
(8) x8=x4+jx5序列及其8点、16点FFT 谱分析
同样地，由图形可知，当N=8时，不满足抽样定理，产生了混叠；当N=16时，满足抽样定理，得到了与理论曲线一致的图形。
5. 误差分析
实验中得到的频谱图与理论基本相符。N 较小时，有可能会不满足采样定理，从而造成频谱混叠失真，产生误差；而当N 较大时，可以比较完整地反映原信号的真实频谱，误差相对较小。
6. 参数选取
由上述实验结果可知，N 要大于信号的长度M ，这样才能由抽样序列不失真地恢复出原信号，故尽可能地增大N 可以有效地反映原信号的信息。
7. 思考题
(1)在N =8时，x 2(n ) 和x 3(n ) 的幅频特性会相同吗? 为什么? N =16呢? N=8时二者的幅频特性相同，N=16时二者的幅频特性不同。 原因：x3可由x2移位得到，而二者不为0的区间长度恰为8，故在N=8时，经周期延拓、移位、截取主值区间，二者的幅频特性相同；而在N=16时，经过了补零，二者便有了区别，幅频特性也就不同了。还注意到，二者的相频特性都不同，这是因为循环位移会影响其傅里叶变换后的相频特性，故不同。
(2)如果周期信号的周期预先不知道， 如何用FFT 进行谱分析? 应先适当截取M 点进行FFT ，再将截取的长度扩大1倍重新截取，比较二者结果，若二者的差别满足分析误差要求，就可以近似得到该信号的频谱，如不满足误差要求则继续将截取的长度加倍进行FFT ，直到结果满足误差要求为止。
8. 实验结论
对离散序列作DFT ，采样时要正确确定采样周期，当周期N 大于离散序列的长度M ，即N ≥M 时，才能得到正确的频谱。当N ＜M 时，
会发生频谱混叠现象，产生失真。
FFT 是DFT 的一种快速算法，其运算次数可大幅度减少，N 越大时减少越明显，因此在实际中常利用FFT 进行计算。
附：实验一全部matlab 代码
%% x1序列的产生与谱分析 clc;clear;close all ;
x1=[1 1 1 1 zeros(1,4)]; %x1序列的产生 subplot(3,2,1:2); stem(0:7,x1); grid on ;
title('x1序列' );
xlabel('n' ); x1_8=fftshift(fft(x1,8)); hx1_8=abs(x1_8);
ax1_8=angle(x1_8); subplot(3,2,3); stem(-4:3,hx1_8); grid on ;
title('N=8时x1的幅度谱' ); xlabel('K' );
ylabel('|Hw|'); subplot(3,2,4); stem(-4:3,ax1_8); grid on ;
title('N=8时x1的相位谱' );
xlabel('K' );ylabel(' 相角' ); x1_16=fftshift(fft(x1,16)); hx1_16=abs(x1_16);
ax1_16=angle(x1_16); subplot(3,2,5); stem(-8:7,hx1_16); grid on ;
title('N=16时x1的幅度谱' );
xlabel('K' );ylabel('|Hw|'); subplot(3,2,6); stem(-8:7,ax1_16); grid on ;
title('N=16时x1的相位谱' );
xlabel('K' );ylabel(' 相角' );
%% x2序列的产生与谱分析 clc;clear;close all ;
x2=[1 2 3 4 4 3 2 1]; subplot(3,2,1:2); stem(0:7,x2);
%画出x1序列 %将x1进行8点的FFT 变换并得到其幅度和相角%画出x1的8点FFT 的幅度谱 %画出x1的16点FFT 的相位谱 %将x1进行16点的FFT 变换并得到其幅度和相角 %画出x1的16点FFT 的幅度谱 %画出x1的16点FFT 的相位谱 %x2序列的产生
title('x2序列' );
xlabel('n' ); %画出x2序列 x2_8=fftshift(fft(x2,8)); hx2_8=abs(x2_8);
ax2_8=angle(x2_8); %将x2进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3); stem(-4:3,hx2_8); grid on ;
title('N=8时x2的幅度谱' );
xlabel('K' );ylabel('|Hw|'); subplot(3,2,4); stem(-4:3,ax2_8); grid on ;
title('N=8时x2的相位谱' );
xlabel('K' );ylabel(' 相角' ); x2_16=fftshift(fft(x2,16)); hx2_16=abs(x2_16);
ax2_16=angle(x2_16); subplot(3,2,5); stem(-8:7,hx2_16); grid on ;
title('N=16时x2的幅度谱' );
xlabel('K' );ylabel('|Hw|'); subplot(3,2,6); stem(-8:7,ax2_16); grid on ;
title('N=16时x2的相位谱' );
xlabel('K' );ylabel(' 相角' );
%% x3序列的产生与谱分析 clc;clear;close all ;
x3=[4 3 2 1 1 2 3 4]; subplot(3,2,1:2); stem(0:7,x3); grid on ;
title('x3序列' );
xlabel('n' ); x3_8=fftshift(fft(x3,8)); hx3_8=abs(x3_8);
ax3_8=angle(x3_8); subplot(3,2,3); stem(-4:3,hx3_8);
%画出x2的8点FFT 的幅度谱 %画出x2的16点FFT 的相位谱 %将x2进行16点的FFT 变换并得到其幅度和相角%画出x2的16点FFT 的幅度谱 %画出x2的16点FFT 的相位谱 %x3序列的产生 %画出x3序列 %将x3进行8点的FFT 变换并得到其幅度和相角
title('N=8时x3的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x3的8点FFT 的幅度谱 subplot(3,2,4); stem(-4:3,ax3_8); grid on ;
title('N=8时x3的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x3的16点FFT 的相位谱 x3_16=fftshift(fft(x3,16)); hx3_16=abs(x3_16);
ax3_16=angle(x3_16); %将x3进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5); stem(-8:7,hx3_16);
grid on ;title('N=16时x3的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x3的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax3_16); grid on ;title('N=16时x3的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x3的16点FFT 的相位谱
%% x4序列的产生与谱分析 clc;clear;close all ;
x4=cos(pi/4*(0:1:19)); %x4序列的产生 subplot(3,2,1:2);stem(0:19,x4); grid on ;title('x4序列' );
xlabel('n' ); %画出x4序列 x4_8=fftshift(fft(x4,8));
hx4_8=abs(x4_8); ax4_8=angle(x4_8); %将x4进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx4_8); grid on ;title('N=8时x4的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x4的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax4_8); grid on ;title('N=8时x4的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x4的16点FFT 的相位谱 x4_16=fftshift(fft(x4,16));
hx4_16=abs(x4_16); ax4_16=angle(x4_16); %将x4进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx4_16); grid on ;title('N=16时x4的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x4的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax4_16); grid on ;title('N=16时x4的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x4的16点FFT 的相位谱
%% x5序列的产生与谱分析 clc;clear;close all ;
x5=sin(pi/8*(0:1:19)); %x5序列的产生 subplot(3,2,1:2);stem(0:19,x5); grid on ;title('x5序列' );
xlabel('n' ); %画出x5序列 x5_8=fftshift(fft(x5,8));
hx5_8=abs(x5_8); ax5_8=angle(x5_8); %将x5进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx5_8); grid on ;title('N=8时x5的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x5的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax5_8); grid on ;title('N=8时x5的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x5的16点FFT 的相位谱 x5_16=fftshift(fft(x5,16));
hx5_16=abs(x5_16); ax5_16=angle(x5_16); %将x5进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx5_16); grid on ;title('N=16时x5的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x5的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax5_16); grid on ;title('N=16时x5的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x5的16点FFT 的相位谱
%% x6序列的产生与谱分析 clc;clear;close all ; fs=64; N=64;
t=(0:1:(N-1))/fs;
x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t); %x6序列的产生 subplot(4,2,1:2);stem(0:63,x6); grid on ;title('x6序列' );
xlabel('n' ); %画出x6序列 x6_16=fftshift(fft(x6,16));
hx6_16=abs(x6_16); ax6_16=angle(x6_16); %将x6进行16点的FFT 变换并得到其幅度和相角 subplot(4,2,3);stem(-8:7,hx6_16); grid on ;title('N=16时x6的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x6的16点FFT 的幅度谱 subplot(4,2,4);stem(-8:7,ax6_16); grid on ;title('N=16时x6的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x6的16点FFT 的相位谱 x6_32=fftshift(fft(x6,32));
hx6_32=abs(x6_32); ax6_32=angle(x6_32); %将x6进行32点的FFT 变换并得到其幅度和相角 subplot(4,2,5);stem(-16:15,hx6_32);
grid on ;title('N=32时x6的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x6的32点FFT 的幅度谱 subplot(4,2,6);stem(-16:15,ax6_32); grid on ;title('N=32时x6的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x6的32点FFT 的相位谱 x6_64=fftshift(fft(x6,64));
hx6_64=abs(x6_64); ax6_64=angle(x6_64); %将x6进行64点的FFT 变换并得到其幅度和相角 subplot(4,2,7);stem(-32:31,hx6_64); grid on ;title('N=64时x6的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x6的64点FFT 的幅度谱 subplot(4,2,8);stem(-32:31,ax6_64); grid on ;title('N=64时x6的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x6的64点FFT 的相位谱
%% x7=x4+x5序列的产生与谱分析 clc;clear;close all ; x4=cos(pi/4*(0:1:19)); x5=sin(pi/8*(0:1:19));
x7=x4+x5; %x7=x4+x5序列的产生 subplot(3,2,1:2);stem(0:19,x7); grid on ;title('x7序列' );
xlabel('n' ); %画出x7序列 x7_8=fftshift(fft(x7,8));
hx7_8=abs(x7_8); ax7_8=angle(x7_8); %将x7进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx7_8); grid on ;title('N=8时x7的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x7的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax7_8); grid on ;title('N=8时x7的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x7的16点FFT 的相位谱 x7_16=fftshift(fft(x7,16));
hx7_16=abs(x7_16); ax7_16=angle(x7_16); %将x7进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx7_16); grid on ;title('N=16时x7的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x7的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax7_16); grid on ;title('N=16时x7的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x7的16点FFT 的相位谱
%% x8=x4+jx5序列的产生与谱分析 clc;clear;close all ; x4=cos(pi/4*(0:1:19));
x5=sin(pi/8*(0:1:19));
x8=x4+1i*x5; %x8=x4+jx5序列的产生 subplot(3,2,1:2);stem(0:19,x8); grid on ;title('x8序列' );
xlabel('n' ); %画出x8序列 x8_8=fftshift(fft(x8,8));
hx8_8=abs(x8_8); ax8_8=angle(x8_8); %将x8进行8点的FFT 变换并得到其幅度和相角 subplot(3,2,3);stem(-4:3,hx8_8); grid on ;title('N=8时x8的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x8的8点FFT 的幅度谱 subplot(3,2,4);stem(-4:3,ax8_8); grid on ;title('N=8时x8的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x8的16点FFT 的相位谱 x8_16=fftshift(fft(x8,16));
hx8_16=abs(x8_16); ax8_16=angle(x8_16); %将x8进行16点的FFT 变换并得到其幅度和相角 subplot(3,2,5);stem(-8:7,hx8_16); grid on ;title('N=16时x8的幅度谱' );
xlabel('K' );ylabel('|Hw|'); %画出x8的16点FFT 的幅度谱 subplot(3,2,6);stem(-8:7,ax8_16); grid on ;title('N=16时x8的相位谱' );
xlabel('K' );ylabel(' 相角' ); %画出x8的16点FFT 的相位谱
实验二 窗函数法设计FIR 数字滤波器
1. 实验目的
(1) 熟悉矩形窗、汉宁窗、海明窗等常用窗函数。
(2) 掌握用上述窗函数法设计FIR 数字滤波器的原理和方法。 (3) 熟悉线性相位FIR 数字滤波器特性。 (4) 了解各种窗函数对滤波特性的影响。
2. 实验原理与方法
如果所希望的滤波器的理想频率响应函数为Hd (e jω) ，则其对应的单位冲激响应为
1
h d (n ) =
2π
⎰π
-
π
H d (e j ω) e j ωn d ω
用窗函数ω(n ) 将h d (n ) 截断，得到:
h (n ) =h d (n ) ω(n )
频率响应函数H (e jω) 为
h (n ) 就作为实际设计的FIR 数字滤波器的单位冲激响应序列， 其
H (e j ω) =∑h (n ) e -j ωn
n =0
N -1
如果要求线性相位特性， 则h (n ) 还必须满足：
h (n ) =±h (N -1-n )
根据上式中的正、负号和长度N 取为奇数或偶数又将线性相位FIR 滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。
第一类：h(n)=h(N-1-n)且N 为奇数，此类情况适合低通、高通、带通、带阻；
第二类：h(n)=h(N-1-n)且N 为偶数，此类情况适合设计低通、带通滤波器；
第三类：h(n)=-h(N-1-n)且N 为奇数，此类情况只适合设计带通滤波器；
第四类：h(n)=-h(N-1-n)且N 为偶数，此类情况适合设计高通、带通滤波器。
3. 实验内容及步骤
(1)复习用窗函数法设计FIR 数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。 (2) 编写程序。
①编写能产生三种窗函数的子程序。
②编写主程序。其中幅度特性要求用分贝（dB ）表示。
设
H (k ) =DFT[h (n )]
H (k ) =H k (k ) +jH I (k )
H (k ) =
画图时，用20log|H(k)|画出幅度特性。
(1) 15阶矩形窗、汉宁窗和海明窗
(2)33阶矩形窗、汉宁窗和海明窗
(1)窗函数法特点
窗口法设计的主要优点是简单，使用方便。窗口函数大多有封闭的公式可循，性能、参数都已有表格、资料可供参考， 计算程序简便， 所以很实用，缺点是通带和阻带的截止频率不易控制。 (2)窗口长度N 和窗函数类型对滤波特性的影响
增加窗的长度，可以减少窗的主瓣宽度，从而减少过渡带的带宽，
但增加N 并不能减少带内波动以及加大阻带。这两个指标只能从窗函数的形状上找解决方法。
6. 综合比较
对选用以上三种窗口设计出的滤波器的结果进行比较可以看出，矩形窗的阻带衰减最小，其过渡带宽也最窄；汉宁窗和海明窗衰减比矩形窗大，其过度带宽比矩形窗大；而海明窗又比汉宁窗阻带衰减大，过渡带也大；从矩形窗、汉宁窗到海明窗，阻带衰减逐渐加大，过渡带宽也逐渐加大。在选定窗口的情况下N 值越大，过渡带越窄，所需的计算量和存储量也加大，因此实际设计的时候应根据需要选用不同的窗口。
7. 实验结论
选择窗函数时要注意：
(1)要从保持最大信息和消除旁瓣的综合效果出发来考虑问题，要使窗
函数频谱中的主瓣宽度尽量窄，能量尽可能集中在主瓣内，从而在谱分析时获得较高的频率分辨力，在数字滤波器设计中获得较小的过渡带。
(2)窗谱的旁瓣高度应尽量小而且随频率尽快衰减，以减小谱估计时泄漏失真。在设计数字滤波器时减小通带的波动，提高阻带的衰减。但主瓣既窄，旁瓣又小衰减又快的窗函数是不容易找到的，比如矩形窗旁瓣很大，但其主瓣宽度是最窄的，因此，在数据处理时通常需要做综合考虑取其折中。
(3)在应用窗函数时，除了要考虑窗谱本身的特性外，还应当充分考虑被分析信号的特点以及具体的处理要求。要考虑信号中的信息量的分布，增强信号中所需要的信息部分，压制信号中不需要的信息部分，以感兴趣的有效信息与窗函数作用后的综台效果为最好来选用窗函数，使得处理结果有足够的频谱检测能力和频谱幅值估计精度。
(4)信号的加窗处理，重要的问题是在于根据信号的性质和研究目的来选用窗函数。如对于频率分辨率要求很高，谱估计幅值精度要求不很高的信号，处理时可选用主瓣较窄而便于分辨的矩形窗。
(5)对于要处理的信号，不知选用哪一种时窗函数好时，往往多用几种窗函数进行处理，然后比较用几种窗处理的结果和试验验证的综合考虑来决定选用什么样的窗函数。
8. 思考题
(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何
用窗函数法设计线性相位低通滤波器? 写出设计步骤。
a) 求出对应的数字频率。
通带截止频率ωp=
阻带截止频率ωst=Ωpfs Ωst
fs
b ）求ℎd(n)。设Hd(ejω) 为理想线性相位滤波器
1sin [ωc(n−τ)] n≠τ1π
−jωτjωnℎd=∫eedω={π(n−τ) ω2π−πc n=τπ
其中τ为线性相位所必须的位移，且满足τ=N−12.
c ）求窗函数。由阻带衰减δ2确定窗形状，由过度带宽确定N 。 d ）求h(n)。由窗函数表达式ω(n)确定FIR 滤波器的h(n)。
h (n ) =ℎd(n)∙ω(n)
e ）由h(n)求Hd(ejω) ，检验各项指标是否满足要求。如不满足要求，则要改变N ，或改变窗形状，然后重新计算。
(2) 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应ℎd(n).
①假设理想带通滤波器的频率响应：
-j ωt 0⎧e H d (e j ω) =⎨⎩00
②对H d (e j ω) 做IDTFT 得单位冲激响应h d (n ) ：
1{sin [(n −t0) ω2]−sin [(n−t0)ω1]} n≠t01π(n−t0) jωjωℎd(n) =∫H(e) ∙edω= 2π−πd 1{π(ω2−ω1) n=t0π
③根据指标要求如通带最大衰减、阻带最小衰减可选择窗函数的类型和点数N 。
h (n ) =h d (n ) ⋅ω(n )
④对h (n ) 做DTFT 可得H j ω
d (e ) ，验证指标是否满足设计要求，若不满
足重新选择窗口类型和窗口长度，直至满足设计指标要求。
附：实验二全部matlab 代码
clc;clear;close all ;
N1=15;N2=33; %分别表示15和33阶 wc=pi/4;
a1=(N1-1)/2;a2=(N2-1)/2;
n1=0:N1-1;n2=0:N2-1;
n1=n1';n2=n2';
Hd1=sin(wc*(n1-a1))./(pi*(n1-a1));Hd1(8)=wc/pi;
Hd2=sin(wc*(n2-a2))./(pi*(n2-a2));Hd2(17)=wc/pi;
w1_boxcar=boxcar(N1);
w2_boxcar=boxcar(N2);
w1_hanning=hanning(N1);
w2_hanning=hanning(N2);
w1_hamming=hamming(N1);
w2_hamming=hamming(N2);
h1_boxcar=Hd1.*w1_boxcar; %最终设计的滤波器系数h(n)=Hd(n)*w(n) h1_hanning=Hd1.*w1_hanning;
h1_hamming=Hd1.*w1_hamming;
h2_boxcar=Hd2.*w2_boxcar;
h2_hanning=Hd2.*w2_hanning;
h2_hamming=Hd2.*w2_hamming;
FFT_N=1024; %做FFT 查看幅频和相频曲线(FFT点数应比滤波器阶数要多，使频谱细化) h1_boxcar_fft=fftshift(fft(h1_boxcar,FFT_N));
h1_boxcar_fft_hw=20*log10(abs(h1_boxcar_fft));
h1_boxcar_fft_fai=angle(h1_boxcar_fft);
h1_hanning_fft=fftshift(fft(h1_hanning,FFT_N));
h1_hanning_fft_hw=20*log10(abs(h1_hanning_fft));
h1_hanning_fft_fai=angle(h1_hanning_fft);
h1_hamming_fft=fftshift(fft(h1_hamming,FFT_N));
h1_hamming_fft_hw=20*log10(abs(h1_hamming_fft));
h1_hamming_fft_fai=angle(h1_hamming_fft);
h2_boxcar_fft=fftshift(fft(h2_boxcar,FFT_N));
h2_boxcar_fft_hw=20*log10(abs(h2_boxcar_fft));
h2_boxcar_fft_fai=angle(h2_boxcar_fft);
h2_hanning_fft=fftshift(fft(h2_hanning,FFT_N));
h2_hanning_fft_hw=20*log10(abs(h2_hanning_fft));
h2_hanning_fft_fai=angle(h2_hanning_fft);
h2_hamming_fft=fftshift(fft(h2_hamming,FFT_N));
h2_hamming_fft_hw=20*log10(abs(h2_hamming_fft));
h2_hamming_fft_fai=angle(h2_hamming_fft);
xias_w=(-512:511)/512;
figure; %画出15阶滤波器 subplot(231);plot(xias_w,h1_boxcar_fft_hw);
grid on ;title('15阶矩形窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(234);plot(xias_w,h1_boxcar_fft_fai);
grid on ;title('15阶矩形窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(232);plot(xias_w,h1_hanning_fft_hw);
grid on ;title('15阶汉宁窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(235);plot(xias_w,h1_hanning_fft_fai);
grid on ;title('15阶汉宁窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(233);plot(xias_w,h1_hamming_fft_hw);
grid on ;title('15阶海明窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(236);plot(xias_w,h1_hamming_fft_fai);
grid on ;title('15阶海明窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
figure; subplot(231);plot(xias_w,h2_boxcar_fft_hw);
grid on ;title('33阶矩形窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(234);plot(xias_w,h2_boxcar_fft_fai);
grid on ;title('33阶矩形窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(232);plot(xias_w,h2_hanning_fft_hw);
grid on ;title('33阶汉宁窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(235);plot(xias_w,h2_hanning_fft_fai);
grid on ;title('33阶汉宁窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(233);plot(xias_w,h2_hamming_fft_hw);
grid on ;title('33阶海明窗幅频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
subplot(236);plot(xias_w,h2_hamming_fft_fai);
grid on ;title('33阶海明窗相频曲线' );
xlabel(' 归一化频率(w/(fs/2))');ylabel(' 幅度衰减(dB)');
%画出33阶滤波器