简述采用窗函数法设计FIR数字滤波器的设计步骤及主要公式。

将模拟频率转化为数字频率，设取样时间为T（要满足抽样定理）

Ωp=2πfpT Ωs=2πfsT

过渡带宽度△Ω=Ωp-Ωs

阻带衰减已经超过74db，要选用Kaiser窗了，Kaiser的参数可变，要根据公式确定滤波器的参数

一般都选用Ⅰ型线性相位滤波器即滤波器阶数M为偶数，程序如下：

wp=;ws=;Ap=1;As=100;

Rp=1-10^(-005Ap);Rs=10^(-005As);

f=[fp fs];

a=[0 1];

dev=[Rp Rs];

[M,wc,beta,ftype]=kaiserord(f,a,dev);

M=mod(M,2)+M;

h=fir1(M,wc,ftype,kaiser(M+1,beta));

omega=linspace(0,pi,512);

mag=freqz(h,[1],omega);

plot(omega/pi,20log10(abs(mag)));

grid;

omega1=linspace(0,wp,512);

h1=freqz(h,[1],omega1);

omega2=linspace(ws,pi,512);

h2=freqz(h,[1],omega2);

fprintf('Ap=%4f\n',-20log10(min(abs(h1))));

fprintf('As=%4f\n',-20log10(max(abs(h2))));

运行程序可以得到滤波器的通阻带衰减，画出频率响应，若同阻带衰减不满足要求还可以使用滤波器的优化，一般使用的等波纹FIR进行优化

MATLAB 信号处理常用函数 

一、 波形产生 

函数名  功能   

sawtooth  产生锯齿波或三角波  

Sinc  产生sinc或函数sin(pit)/(pit) 

Square  产生方波  

Diric  产生Dirichlet或周期sinc函数

二、 滤波器分析和实现 

函数名  功能   

Abs  求绝对值（幅值）  

Freqs  模拟滤波器频率响应 

Angle  求相角  

Freqspace  频率响应中的频率间隔 

Conv  求卷积  

Freqz  数字滤波器频率响应 

Fftfilt  重叠相加法FFT滤波器实现  

Grpdelay  平均滤波器延迟（群延迟） 

Filter  直接滤波器实现  

Impz  数字滤波器的冲激响应 

Filtfilt  零相位数字滤波  

Zplane  离散系统零极点图 

Filtie  Filter 函数初始条件选择

三、 线性系统变换 

函数名  功能  

Convmtx  卷积矩阵  

Ss2tf    变系统状态空间形式为传递函数形式 

Ploy2rc  从多项式系数中计算反射系数 

Ss2zp   变系统状态空间形式为零极点增益形式 

Rc2ploy  从反射系数中计算多项式系数 

Tf2ss   变系统传递函数形式为状态空间形式 

Residuez   Z变换部分分式展开或留数计算 

Tf2zp     变系统传递函数形式为零极点增益形式 

Sos2ss    变系统二阶分割形式为状态空间形式 

Zp2sos    变系统零极点形式为二阶分割形式 

Sos2zp    变系统二阶分割形式为零极点增益形式 

Zp2tf     变系统零极点增益形式为传递函数形式 

Ss2sos    变系统状态空间形式为二阶分割形式

   

四、 IIR滤波器设计 

Besself   Bessel（贝塞尔）模拟滤波器设计 

Cheby2   Chebyshev（切比雪夫）II型模拟滤波器设计 

Butter    Butterworth（巴特沃思）模拟滤波器设计 

Ellip    椭圆模拟滤波器设计 

Cheby1  Chebyshev（切比雪夫）I 型模拟滤波器设计 

Yulewalk  递归数字滤波器设计

五、 IIR滤波器阶选择 

Buttord    Butterworth（巴特沃思）滤波器阶的选择 

Cheb2ord   Chebyshev（切比雪夫）II型滤波器阶的选择 

Ehebord   Chebyshev（切比雪夫）I 型滤波器阶的选择 

Clipord  椭圆滤波器设计阶的选择 模拟原型滤波器设计 

Besselap  Bessel模拟低通滤波器原型  

Cheb2ap   Chebyshev（切比雪夫）II型低通滤波器原型 

Buttap    Butterworth（巴特沃思）模拟低通滤波器原型 

Ellipap  椭圆模拟低通滤波器原型 

Cheb1ap  Chebyshev（切比雪夫）I 型低通滤波器原型

    

六、 频率变换 

Lp2bp  低通到带通模拟滤波器转换  

Lp2bs  低通到带阻模拟滤波器变换 

Lp2hp  低通到高通模拟滤波器变换  

Lp2lp  低通到低通模拟滤波器转换

七、 滤波器离散化 

Blinear  双线性变换  

Impinvar  冲激响应不变法

八、 FIR滤波器设计 

Fir1 基于窗函数的 FIR 滤波器设计—标准响应 

Intfilt  内插FIR滤波器设计 

Fir2   基于窗函数的 FIR 滤波器设计—任意响应 

Remez  Firls  最小二乘FIR滤波器设计  

Remezord  Parks-McCellan 最优 FIR 滤波器 j阶估计

九、 窗函数 

Boxcar  矩形窗  

Hanning  Hanning（汉宁）窗 

Triang  三角窗  

Blackman  Blackman（布莱克曼）窗 

Bartlett  Bartlett（巴特得特）窗  

Chebwin  Chebyshev（切比雪夫）窗 

Hamming  Hamming（汉明）窗  

Kaiser  Kaiser（凯泽）窗

十、 变换 

Ctz  线性调频Z变换  

Fft  一维快速傅里叶变换 

Dct  离散余弦变换  

Ifft  一维快速傅里叶逆变换 

Idct  逆离散余弦变换  

Fftshift  重新排列 fft的输出 

Dftmtx  离散傅里叶变换矩阵  

Hilbert  Hilbert（希尔伯特）变换

 

十一、 统计信号处理 

Cov  协方差矩阵  

Psd  信号功率谱密度（PSD）估计 

Xcov  互协方差函数估计  

Tfe  从输入输出中估计传递函数 

Corrcoef  相关系数矩阵  

Periodogram  采用周期图法估计功率谱密度 

Xcoor  互相关系数估计  

Pwelch  采用 Welch方法估计功率谱密度 

Cohere  相关函数平方幅值估计  

Rand  生成均匀分布的随机数 

Csd  互谱密度估计  

Randn  生成正态分布的随机数

十二、 自适应滤波器部分 

Adaptfiltlms  最小均方(LMS)自适应算法  

Adaptfiltrls  递推最小二乘(RLS)自适应算法 

Adaptfiltnlms 归一化最小均方(NLMS)自适应算法

十三、 时频分析与小波变换部分 

Spectrogram  短时傅里叶变换  

Idwt  单级离散一维小波逆变换 

Waveinfo  介绍小波工具箱中所有小波的信息 

Wavedec  多级离散一维小波分解 

Cwt  连续一维小波变换  

Appcoef  一维小波变换近似系数 

Dwt  单级离散一维小波变换  

Detcoef  一维小波变换细节系数

十四、 二维信号处理 

Conv2  二维卷积  

Xcorr2  二维互相关参数 

Fft2  二维快读傅里叶变换  

Dwt2  单级离散二维小波变换 

Ifft2  二维逆快速傅里叶变换  

Idwt2  单级离散二维小波逆变换 

Filter2  二维数字滤波器  

Waverec2  多级离散二维小波分解

1FIR数字滤波器原理

假设理想低通滤波器的截止频率为ωc＝2πfc，且具有线性相位，群延时为a，即频率响应：

航空重力勘探理论方法及应用

表示成幅度函数和相位函数形式：

航空重力勘探理论方法及应用

则幅度函数：

航空重力勘探理论方法及应用

在通带范围｜ω｜ ≤ωc（截止频率）内Hd（ejω）的幅度为1，相位为－ωα；对应的时间域（或空间域）滤波函数为：

航空重力勘探理论方法及应用

有限脉冲响应FIR（Finite Impulse Response）数字滤波器要求用有限长的单位冲击响应h（n）来逼近无限长的理想滤波器的单位冲击响应hd（n），最常用和有效的方法就是用一个有限长（长度为N）的“窗函数”序列w（n）来截取hd（n）的主要成分（陈玉东，2005）：

航空重力勘探理论方法及应用

实际上是用有限长的h（n）去逼近hd（n），通过这种方式得到的频率响应H（ejω）近似于理想频率响应Hd（ejω）（在频率域内采用均方差最小准则逼近）。按照线性相位滤波器的约束要求，h（n）必须是偶对称的，其对称中心应为它长度的一半：h（n）=h（N-1-n），而且 ；所以同时要求窗函数w（n）也必须是关于中心偶对称：w（n）=w（N-1-n）。

2几种常见窗函数

（1）矩形窗

长度为N的矩形窗函数为：

航空重力勘探理论方法及应用

（2）三角形窗（Bartlett）

长度为N的三角形窗函数为：

航空重力勘探理论方法及应用

（3）汉宁窗（Hanning）

长度为N的汉宁窗函数为：

航空重力勘探理论方法及应用

（4）海明窗（Hamming）

为使得旁瓣更小，可将汉宁窗改进成海明窗，长度为N的海明窗函数为：

航空重力勘探理论方法及应用

（5）布拉克曼窗（Blackman）

为进一步有效抑制旁瓣，可以再加上余弦的二次谐波分量，得到长度为N的布拉克曼窗函数为：

航空重力勘探理论方法及应用

（6）凯泽窗（Kaiser）

长度为N的凯泽窗函数为：

航空重力勘探理论方法及应用

其中I0（x）为第一类变形零阶贝塞尔函数，α＝（N-1）/2。β是一个可自由选择的参数，它可同时调整窗函数谱主瓣宽度与旁瓣幅值；β越大，则窗函数w（n）变化越快、变得越窄，频谱旁瓣就越小，但主瓣宽度相应增加。一般选择4＜β＜9，相当于窗函数频谱旁瓣幅度与主瓣幅度的比值由31%变到0047%。β＝0时相当于矩形窗（陈玉东，2005）。

3窗函数FIR滤波器

式（7-4-3）至式（7-4-8）窗函数都满足关于中心偶对称的线性相位滤波器的约束要求，结合式（7-4-1）至式（7-4-2）可以得到相应窗函数的FIR低通数字滤波器函数（郭志宏，罗锋，等，2007）:

航空重力勘探理论方法及应用

航空重力勘探理论方法及应用

用该滤波器窗口对时间域（或空间域）长度为M的数据序列逐点进行窗口滑动卷积求和计算（实际处理时窗口中点作为输出计算点，则一边损失半个滤波窗口数据），就可获得FIR滤波后的数据（郭志宏，段树岭，等，2009）：

航空重力勘探理论方法及应用

h（n）为滤波器系数，x（n）、y（n）分别为输入、输出数据序列。

4窗函数FIR滤波试验

（1）GT-1A型航空重力数据

图7-4-1至图7-4-2分别为GT-1A型航空重力系统获得的一条原始未滤波、100 s和60 s滤波自由空间重力异常测线数据，其中飞机的飞行速度约60m/s，剖面图横轴为测线基准点号，基准点间距约30 m。图7-4-1中GT-1A型系统航空原始未滤波自由空间重力测线数据的高频干扰非常之严重，噪声幅度在－5 000×10-5m·s-2至5000×10-5m·s-2的大范围内变化，而幅度通常只有（10-3～10-4）m·s-2的由密度和构造变化等地质因素引起的重力异常信号（图7-4-2）则完全淹没在高频干扰中。图7-4-2中GT-1A型系统航空100 s、60 s滤波自由空间重力测线数据是采用GT-1A型航空重力系统自带软件模块由图7-4-1的航空原始未滤波自由空间重力测线数据获得的滤波数据，滤波后高频干扰已基本消除，油气和矿产地球物理勘查所需的重力异常则较好的显现出来。

图7-4-1 GT-1A型航空重力系统原始未滤波自由空间重力异常

（2）几种窗函数FIR滤波试验

根据式（7-4-1）至式（7-4-10），我们研制了窗函数法FIR数字滤波计算软件，用各种窗函数FIR滤波器对图7-4-1的GT-1A航空原始未滤波自由空间重力测线数据分别进行了截止波长为100 s、60 s长度（按v=60m/s的航速计算，截止波长A。分别为6km、36km，按fc=v/λc计算的截止频率分别为001 Hz、00167 Hz）的低通滤波试验计算，试验结果见图7-4-3至图7-4-8。为了图形对比方便，各剖面图中仍然保留了测线边部两端的半个滤波窗口数据，这些数据由于存在边部效应，因而是不准确的，实际应用时应该去掉。从试验结果图可以看到，矩形窗和三角窗FIR滤波后异常整体形状虽然也与图7-4-2类似，但其上叠加了高频扰动，尤其是矩形窗FIR滤波结果，这就是通常所说的“吉布斯”振荡效应（陈玉东，2005）。如果在图7-4-3至图7-4-4的基础上，采用空间域非线性曲率滤波方法（郭志宏，刘浩军，等，2003），用中国国土资源航空物探遥感中心的“空中探针”系统（刘浩军，薛典军，等，2003）中的滤波软件进一步处理，则可获得消除扰动后接近图7-4-2效果的异常数据。从汉宁窗、海明窗、布拉克曼窗以及凯泽窗FIR滤波试验结果看到，通过选择合适的窗口长度、截至波长等滤波参数，基本都获得了令人满意的效果。

图7-4-2 GT-1A型航空重力系统100s、60s滤波自由空间重力异常

图7-4-3 矩形窗FIR低通滤波截止波长100s、60s航空自由空间重力异常

表7-4-1为图7-4-3至图7-4-8所示的各种窗函数FIR低通滤波截止波长100 s、60 s长度航空自由空间重力异常与图7-4-2所示的GT-1A型航空重力系统100 s、60 s滤波自由空间重力异常（作为标准）的比较，通过两者之差值的统计结果来衡量吻合程度。从统计表中可以看到，除了矩形窗、三角窗外，其他几种窗函数FIR低通滤波结果的差异值都在±1×10-5m·s-2以内，均方差值则多数为03×10-5m·s-2左右，可见吻合程度还是比较好的。

图7-4-4 三角窗FIR低通滤波截止波长100s、60s航空自由空间重力异常

5结论

1）通过选择合适的窗形、窗口长度、滤波参数，窗函数法FIR低通数字滤波器可以在航空重力数据的滤波处理中发挥应有的作用。

图7-4-5 汉宁窗FIR低通滤波截止波长100s、60s航空自由空间重力异常

图7-4-6 海明窗FIR低通滤波截止波长100s、60s航空自由空间重力异常

图7-4-7 布拉克曼窗FIR低通滤波截止波长100s、60s航空自由空间重力异常

图7-4-8 凯泽窗（β＝6）FIR低通滤波截止波长100s、60s航空自由空间重力异常

2）为了获得与GT-1A型航空重力系统100 s、60 s低通滤波（60m/s航速）对应的自由空间重力测线数据，所选择汉宁、海明、布拉克曼、凯泽窗的长度通常为400点（2 Hz采样率），FIR低通滤波对应的截止频率分别为001 Hz、00167 Hz。

3）窗函数法不但可以设计FIR低通滤波器，还可设计FIR高通、带通、带阻滤波器等。通常一个高通滤波器相当于一个全通滤波器减去一个低通滤波器；一个带通滤波器相当于两个低通滤波器相减；而一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器。

表7-4-1 窗函数FIR滤波试验结果与GT-1A系统滤波结果的差值统计

4）除了窗函数法FIR低通滤波器，其他诸如等波纹法FIR低通滤波器、无限脉冲响应IIR低通滤波、Kalman滤波等方法（周坚鑫，刘浩军，等，2001；陈玉东，2005）均可用于航空重力数据的低通数字滤波处理中。

无限脉冲响应滤波器是数位滤波器的一种，简称IIR数位滤波器(infinite impulse response filter)。由于无限脉冲响应滤波器中存在反馈回路，因此对于脉冲输入信号的响应是无限延续的。

有限脉冲响应滤波器是数字滤波器的一种，简称FIR数字滤波器（finite impulse response filter）。这类滤波器对于脉冲输入信号的响应最终趋向于0，因此是有限的，而得名。它是相对于无限脉冲响应滤波器(IIR)而言。

有限脉冲响应滤波器（FIR filter）的优点：

1 脉冲响应（impulse response）为有限长：造成当输入数位讯号为有限长的时候，输出数位讯号也为有限长。

2 比无限脉冲响应滤波器（IIR filter）较容易最佳化（optimize）。

3 线性相位（linear phase）：造成h(n)\,是偶对称（even）或奇对称（odd）且有限长。

4 一定是稳定的（stable）：因为Z转换（Z transform）后所有的极点（pole）都在单位圆内。

有限脉冲响应滤波器（FIR filter）的缺点：

设计方式较无限脉冲响应滤波器（IIR filter）不容易。

无限脉冲响应滤波器(IIR filter)的优点：

较容易设计以及实现。

无限脉冲响应滤波器(IIR filter)的缺点：

1 脉冲响应(impulse response)为无限长：造成当输入数位讯号为有限长的时候，输出数位讯号会变成无限长。

2 比有限脉冲响应滤波器(FIR filter)较不易最佳化(optimize)。

3 不一定是稳定的(stable)：因为Z转换(Z transform)后所有的极点(pole)不一定都在单位圆内。

reference：w开头的被baidu屏蔽的某网站