基于Matlab设计FIR滤波器
一、摘要 除了采用编程的方法实现滤波器之外，Matlab中自带工具箱FDATool也能很方便快速的实现滤波器的设计。另外FPGA、DSP等实现数字滤波器算法时，经常要用到滤波器系数，采用FDATool工具箱也能快速的得到滤波器系数。二、实验平台 Matlab7.1
三、实现步骤3.1 滤波器指标若需要设计一个16阶的FIR滤波器（h(0)=0），给定的参数如下：(1) 低通滤波器(2) 采样频率Fs为48kHz，滤波器Fc为10.8kHz(3) 输入序列位宽为9位（最高位为符号位）在此利用MATLAB来完成FIR滤波器系数的确定。3.2 打开MATLAB的FDAToolMATLAB集成了一套功能强大的滤波器设计工具FDATool（Filter Design & Analysis Tool），可以完成多种滤波器的设计、分析和性能评估。单击MATLAB主窗口下方的“Start”按钮，如图B.1所示，选择菜单“ToolBox” →“Filter Design” →“Filter Design & Analysis Tool（FDATool）”命令，打开FDATool，如图B.2所示。
http://pic002.cnblogs.com/images/2011/314385/2011112210041151.jpg
图B.1 FDATool的启动
http://pic002.cnblogs.com/images/2011/314385/2011112210081837.jpg
图B.2 FDATool的主界面
另外，在MATLAB主命令窗口内键入“fdatool”，同样可打开FDATool程序界面。3.3 选择Design FilterFDATool界面左下侧排列了一组工具按钮，其功能分别如下所述：● 创建多速率滤波器（Create a Multirate Filter）● 滤波器转换（TransForm Filter）● 设置量化参数（Set Quantization Parameters）● 实现模型（Realize Model）● 零极点编辑器（Pole-zero Editor）● 导入滤波器（Import Filter）● 设计滤波器（Design Filter）选择其中的按钮，进入设计滤波器界面，进行下列选择，如图B.3所示。http://pic002.cnblogs.com/images/2011/314385/2011112210105662.jpg
图B.3 FDATool设计FIR滤波器
● 滤波器类型（Filer Type）为低通（Low Pass）● 设计方法（Design Method）为FIR，采用窗函数法（Window）● 滤波器阶数（Filter order）定制为15● 窗口类型为Kaiser，Beta为0.5● Fs为48kHz，Fc为10.8kHz最后单击Design Filter图标，让MATLAB计算FIR滤波器系数并作相关分析。其系统函数H(z)可用下式来表示：H(z)=显然上式可以写成:H(z)=即可以看成是一个15阶的FIR滤波器的输出结果经过了一个单位延时单元，所以在FDATool中,把它看成15阶FIR滤波器来计算参数。因此，设置滤波器阶数时，要比要求的小1。3.4 滤波器分析计算完FIR滤波器系数以后，往往需要对设计好的FIR滤波器进行相关的性能分析，以便了解该滤波器是否满足设计要求。分析操作步骤如下：选择FDATool的菜单“Analysis”→“Magnitude Response”，启动幅频响应分析如图B.4所示，x轴为频率，y轴为幅度值（单位为dB）。http://pic002.cnblogs.com/images/2011/314385/2011112210220424.jpg
图B.4 FIR滤波器幅频响应在图的左侧列出了当前滤波器的相关信息：● 滤波器类型为Direct Form FIR（直接I型FIR滤波器）● 滤波器阶数为15选择菜单“Analysis”→“Phase Response”，启动相频响应分析，如图B.5所示。由该图可以看到设计的FIR滤波器在通带内其相位响应为线性的，即该滤波器是一个线性相位的滤波器。http://pic002.cnblogs.com/images/2011/314385/2011112210232724.jpg
图B.5 滤波器相频响应图B.6显示了滤波器幅频特性与相频特性的比较，这可以通过菜单“Analysis”→“Magnitude and Phase Response”来启动分析。http://pic002.cnblogs.com/images/2011/314385/2011112210250033.jpg
图B.6 滤波器幅频和相频响应选择菜单“Analysis”→“Group Delay Response”，启动群时延分析。FDATool还提供了以下几种分析工具：● 群时延响应分析。● 冲激响应分析（Impulse Response），如图B.7所示。● 阶跃响应分析（Step Response），如图B.8所示。● 零极点图分析（Pole/Zero Plot），如图B.9所示。http://pic002.cnblogs.com/images/2011/314385/2011112210255589.jpg
图B.7 冲激响应http://pic002.cnblogs.com/images/2011/314385/2011112210264654.jpg
图B.8 阶跃响应http://pic002.cnblogs.com/images/2011/314385/2011112210274118.jpg
图B.9 零极点图求出的FIR滤波器的系数可以通过选择菜单“Analysis”→“Filter Coefficients”来观察。如图B.10所示，图中列出了FDATool计算的15阶直接I型FIR滤波器的部分系数。http://pic002.cnblogs.com/images/2011/314385/2011112210284492.jpg
图B.10 滤波器系数3.5 量化可以看到，FDATool计算出的值是一个有符号的小数，如果建立的FIR滤波器模型需要一个整数作为滤波器系数，就必须进行量化，并对得到的系数进行归一化。为此，单击FDATool左下侧的工具按钮进行量化参数设置。量化参数有三种方式：双精度、单精度和定点。在使用定点量化前，必须确保MATLAB中已经安装定点工具箱并有相应的授权。3.6 导出滤波器系数为导出设计好的滤波器系数，选择FDATool菜单的“File”→“Export”命令，打开Export（导出）对话框，如图B.11所示。http://pic002.cnblogs.com/images/2011/314385/2011112210294755.jpg
图B.11 滤波器系数Export对话框在该窗口中，选择导出到工作区（Workplace）。这时滤波器系数就存入到一个一维变量Num中了。不过这时Num中的元素是以小数形式出现的：Num=Columns 1 through 9-0.03690.01090.05580.0054-0.0873-0.04840.18050.41330.4133Columns 10 through 160.1805-0.0484-0.08730.00540.05580.0109-0.0369 由此，可以得到低通滤波器的系数。
资料来自网络 仅供大家联系使用
往往到最后一步滤波器系数就存入到一个一维变量Num大家可能就不知道到哪儿去找到系数
很简答就是在Command Window中输入Num回车就显示啦 接下来给大家一个历程 虽然写的是c28的由于不使用外设所以基于任何一款DSP都可以运行
一、摘要 采用DSP做FIR算法
二、实验平台 Matlab7.1 + CCS3.1
三、实验内容 根据要求设计低通FIR滤波器。
要求：通带边缘频率10KHz，阻带边缘频率22KHz，阻带衰减75dB，采样频率50KHz。
四、实验步骤3.1参数计算窗函数选定：阻带衰减75dB，选择blackman窗
截止频率：2pi*（10+（22-10）/2）/50 = 0.64pi
窗函数长度：blackman窗的过渡带宽为5.98，单位为2pi/N，而要设计的低通滤波器的过渡带宽为2pi*12/50=0.48pi，二者相等，得N=24.9，取25。
3.2滤波器的脉冲响应理想低通滤波器脉冲响应：
h1 = sin(nΩ1)/n/pi = sin(0.64pi*n)/n/pi
窗函数为：
w = 0.42 - 0.5cos(2pi*n/24) + 0.8cos(4pi*n/24)
则滤波器脉冲响应为：
h = h1*w |n|<=12
h = 0 |n|>12
3.3滤波器的差分方程根据滤波器的脉冲响应计算出h，然后将脉冲响应值移位为因果序列
这里计算h的值，采用Matlab计算，有2种方法。
一种是用程序，代码如下：
Window=blackman(25);
b=fir1(24,0.64,Window);
freqz(b,1)
系数如下：
h1 =Columns 1 through 8 0.000 0.000 0.001 -0.002 -0.002 0.010 -0.009 -0.018Columns 9 through 16 0.049 -0.020 -0.110 0.280 0.640 0.280 -0.110 -0.020
Columns 17 through 24 0.049 -0.018 -0.009 0.010 -0.002 -0.002 0.001 0.000
Columns 25 0.000
另外一种是通过FDATool工具箱设置参数之后生成滤波器(这里设置滤波器的阶数时，设置为24)，之后得到滤波器的系数。系数值同上。最后，得到滤波器的差分方程：y = 0.001x - 0.002x - 0.002x + 0.01x - 0.009x - 0.018 + 0.049x -0.02x - 0.11x + 0.28x + 0.64x + 0.28x - 0.11 - 0.02x + 0.049x - 0.018x - 0.009x + 0.1x - 0.002x - 0.002x + 0.001x3.4DSP实现根据滤波器系数，编写DSP实现的程序，下面是采用CCS软件仿真的形式，验证FIR滤波器的效果。
//#include "DSP281x_Device.h" // DSP281x Headerfile Include File
//#include "DSP281x_Examples.h" // DSP281x Examples Include File
#include "f2812a.h"
#include"math.h"
#define FIRNUMBER 25
#define SIGNAL1F 1000
#define SIGNAL2F 4500
#define SAMPLEF10000
#define PI 3.1415926
float InputWave();
float FIR();
float fHn={ 0.0,0.0,0.001,-0.002,-0.002,0.01,-0.009,
-0.018,0.049,-0.02,0.11,0.28,0.64,0.28,
-0.11,-0.02,0.049,-0.018,-0.009,0.01,
-0.002,-0.002,0.001,0.0,0.0
};
float fXn={ 0.0 };
float fInput,fOutput;
float fSignal1,fSignal2;
float fStepSignal1,fStepSignal2;
float f2PI;
int i;
float fIn,fOut;
int nIn,nOut;
main(void)
{
nIn=0; nOut=0;
f2PI=2*PI;
fSignal1=0.0;
fSignal2=PI*0.1;
fStepSignal1=2*PI/30;
fStepSignal2=2*PI*1.4;
while ( 1 )
{
fInput=InputWave();
fIn=fInput;
nIn++; nIn%=256;
fOutput=FIR();
fOut=fOutput;
nOut++;
if ( nOut>=256 )
{
nOut=0;/* 请在此句上设置软件断点 */
}
}
}
float InputWave()
{
for ( i=FIRNUMBER-1;i>0;i-- )
fXn=fXn;
fXn=sin(fSignal1)+cos(fSignal2)/6.0;
fSignal1+=fStepSignal1;
if ( fSignal1>=f2PI ) fSignal1-=f2PI;
fSignal2+=fStepSignal2;
if ( fSignal2>=f2PI ) fSignal2-=f2PI;
return(fXn);
}
float FIR()
{
float fSum;
fSum=0;
for ( i=0;i<FIRNUMBER;i++ )
{
fSum+=(fXn*fHn);
}
return(fSum);
}
3.5仿真结果
http://pic002.cnblogs.com/images/2011/314385/2011112416103655.jpg
左上角的波形，为混叠有高频噪声的输入波形；右上角的是该波形的频谱。
左下角的波形，为经过FIR低通滤波器之后的输出波形；右下角的是该波形的频谱。
由实验结果可知，FIR滤波器能起到很好的滤波效果。 一、摘要 前面一篇**介绍了通过FDATool工具箱实现滤波器的设计，这里通过几个例子说明采用Matlab语言设计FIR滤波器的过程。
二、实验平台 Matlab7.1
三、实验原理 以低通滤波器为例，其常用的设计指标有：
[*]通带边缘频率fp(数字频率为Ωp)
[*]阻带边缘频率fst (数字频率为Ωst)
[*]通带内最大纹波衰减δp=-20log10(1-αp)，单位为 dB
[*]阻带最小衰减αs=-20log10(αs)，单位为 dB
[*]阻带起伏αs
[*]通带峰值起伏αp
其中，以1、2、3、4条最为常用。5、6条在程序中估算滤波器阶数等参数时会用到。
数字频率 = 模拟频率/采样频率
四、实例分析例1用凯塞窗设计一FIR低通滤波器，通带边界频率Ωp=0.3pi，阻带边界频率 Ωs=0.5pi，阻带衰减δs不小于50dB。方法一：手动计算滤波器阶数N和β值，之后在通过程序设计出滤波器。第一步：通过过渡带宽度和阻带衰减，计算滤波器的阶数B和β值。http://pic002.cnblogs.com/images/2011/314385/2011112316182982.jpg
第二步：通过程序设计滤波器。程序如下：b = fir1(29,0.4,kaiser(30,4.55));=freqz(b,1);plot(w1/pi,20*log10(abs(h1))); axis(); grid;xlabel('归一化频率/p') ;ylabel('幅度/dB') ;
波形如下：http://pic002.cnblogs.com/images/2011/314385/2011112316205032.jpg
方法二：采用 = kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器。 这里的函数kaiserord(f,a,dev)或者kaiserord(f,a,dev,fs)： f为对应的频率，fs为采样频率；当f用数字频率表示时，fs则不需要写。 a=为由f指定的各个频带上的幅值向量，一般只有0和1表示；a和f长度关系为（2*a的长度）- 2=（f的长度） devs=用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等，计算公式：阻带衰减误差=αs，通带衰减误差=αp，可有滤波器指标中的3、4条得到。 fs缺省为2Hz。
程序如下：fcuts = ; %归一化频率omega/pi，这里指通带截止频率、阻带起始频率mags = ;devs = ; = kaiserord(fcuts,mags,devs);%计算出凯塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); freqz(hh); 波形如下：http://pic002.cnblogs.com/images/2011/314385/2011112316234228.jpg
实际中，一般调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W(ω)，调用函数remez可进行等波纹滤波器的设计，直接求出滤波器系数。函数remezord中的数组fedge为通带和阻带边界频率，数组mval是两个边界处的幅值，而数组dev是通带和阻带的波动，fs是采样频率单位为Hz。 例2利用雷米兹交替算法设计等波纹滤波器，设计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动 阻带最小衰减At=40dB，采样频率fs=4000Hz。
解：在MATLAB中可以用remezord 和remez两个函数设计 程序如下： fedge=;
mval=;
dev=;
fs=4000;
=remezord(fedge,mval,dev,fs);
b=remez(N,fpts,mag,wt);
=freqz(b,1,256);
plot(w*2000/pi,20*log10(abs(h)));
grid;
xlabel('频率/Hz') ;
ylabel('幅度/dB'); 波形如下：http://pic002.cnblogs.com/images/2011/314385/2011112320210539.jpg
例3 利用MATLAB编程设计一个数字带通滤波器，指标要求如下：通带边缘频率：Ωp1=0.45pi，Ωp2=0.65pi，通带峰值起伏：δ1<=1。阻带边缘频率：Ωs1=0.3pi，Ωs2=0.8pi，最小阻带衰减：δ2>=40 。方法一：窗函数法 程序如下：=kaiserord(,,);%用kaiserord函数估计出滤波器阶数n和beta参数h1=fir1(n,wn,ftype,kaiser(n+1,bta),'noscale');=freqz(h1,1,256);figure(1)subplot(2,1,1)plot(w1/pi,20*log10(abs(hh1)))gridxlabel('归一化频率w');ylabel('幅度/db');subplot(2,1,2)plot(w1/pi,angle(hh1))gridxlabel('归一化频率w');ylabel('相位/rad'); 波形如下： http://pic002.cnblogs.com/images/2011/314385/2011112320424359.jpg
滤波器系数为：
h1 = Columns 1 through 8 0.0041 0.0055 -0.0091 -0.0018 -0.0056 -0.0000 0.0391 -0.0152 Columns 9 through 16 -0.0381 0.0077 -0.0293 0.0940 0.0907 -0.2630 -0.0517 0.3500 Columns 17 through 24 -0.0517 -0.2630 0.0907 0.0940 -0.0293 0.0077 -0.0381 -0.0152 Columns 25 through 31 0.0391 -0.0000 -0.0056 -0.0018 -0.0091 0.0055 0.0041
如果直接用freqz(h1,1,256)，得幅频特性和相频特性曲线：http://pic002.cnblogs.com/images/2011/314385/2011112320474023.jpg 方法二：等波纹法设计
程序如下：
=remezord(,,);%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w，使remez函数设计出的滤波器满足f、a及dev指定的性能要求。
h2=remez(n,fpts,mag,wt);%设计出等波纹滤波器
=freqz(h2,1,256);
figure(2)
subplot(2,1,1)
plot(w2/pi,20*log10(abs(hh2)))
grid
xlabel('归一化频率w');ylabel('幅度/db');
subplot(2,1,2)
plot(w2/pi,angle(hh2))
grid
xlabel('归一化频率w');ylabel('相位/rad');
h2
波形如下：
http://pic002.cnblogs.com/images/2011/314385/2011112320550677.jpg
滤波器系数如下：
h2 = Columns 1 through 9 -0.0013 0.0092 -0.0255 -0.0642 0.1177 0.0922 -0.2466 -0.0466 0.3116 Columns 10 through 17 -0.0466 -0.2466 0.0922 0.1177 -0.0642 -0.0255 0.0092 -0.0013
如果直接用freqz(h2,1,256);得幅频特性和相频特性曲线：http://pic002.cnblogs.com/images/2011/314385/2011112320564685.jpg
方法三：采用FDATool工具 这种方法需要事先计算出滤波器的阶数，bate值，然后设置相应参数，最后生成滤波器。
设置界面如下图所示：
http://pic002.cnblogs.com/images/2011/314385/2011112409551425.jpg
将上述圈圈的区域设置好之后，生成滤波器，最后通过analysis菜单可以观察生成的滤波器的各种特性曲线和滤波器系数。这里的滤波器系数跟方法一的一样。
波形如下：
http://pic002.cnblogs.com/images/2011/314385/2011112410002856.jpg
五、结果分析5.1滤波器设计总结 FIR滤波器实现一般采用窗函数法和等纹波设计法。窗函数法还包含两个分支，一种是用公式先手动算出N值和其他对应得窗函数参数值，再代入窗函数和fir1实现，一种是用函数*rord估算出N和相应参数再用fir1实现。不过要注意*rord会低估或高估阶次n，可能会使滤波器达不到指定的性能，这时应稍微增加或降低阶次。如果截止频率在0或Nyquist频率附近，或者设定的dev值较大，则得不到正确结果。 滤波器实现形式及特点：由于一般的滤波器在利用窗函数是其通带波纹和阻带波纹不同（一般为第一个阻带波纹最大）因此，在满足第一个阻带衰减旁瓣时，比其频率高的旁瓣，它们的衰减都大大超出要求。而根据阻带衰减与项数的近似关系N = P(δ2)*fs/TW，可得当阻带衰减越大，所需项数越多。5.2窗函数法和等波纹设计的不同之处 窗函数设计是通过最小平方积分办法来设计的，即该滤波器的误差为：
http://pic002.cnblogs.com/images/2011/314385/2011112410164455.jpg
即要求最小方法来设计滤波器，这样的滤波器更忠实于理想滤波器（即滤波系数更接近于理想滤波器）。证明如下：http://pic002.cnblogs.com/images/2011/314385/2011112410184087.jpg
因此，幅度频谱差值越小，实际滤波器就越接近理想滤波器。 而等波纹滤波器是通过最大加权误差最小化来实现，其误差为： http://pic002.cnblogs.com/images/2011/314385/2011112410194334.jpg
要求该误差最小来实现滤波器，得出来的滤波系数较窗函数设计相差较远。以下通过对例3中的h1及h2作比较。%sigsum是用来对数组各元素进行求和function y=sigsum(n1,n2,n,x);y=0; for i=n1+1-min(n):n2+1-min(n) y=y+x(i);end n=0.001:30.001;h=2*cos(0.55*pi*(n-15)).*sin(0.175*pi*(n-15))./(pi*(n-15));delta1=h-h1;n=0.001:16.001;h=2*cos(0.55*pi*(n-15)).*sin(0.175*pi*(n-15))./(pi*(n-15));delta2=h-h2;y1=sigsum(0,30,,(abs(delta1).^2))/31;y2=sigsum(0,16,,(abs(delta2).^2))/17; 结果如下：y1 = 1.9099e-004 y2 = 0.0278 由此得到用窗函数实现的滤波系数比用等波纹滤波器系数的每一项更接近于理想滤波器（y1为用窗函数实现的与理想滤波器的差值，y2为用等波纹滤波器实现的与理想滤波器的差值）； http://pic002.cnblogs.com/images/2011/314385/2011112410260532.jpg
对比二者的幅度频谱可知，等波纹滤波器阻带边缘比用窗函数实现的更平滑（理想滤波器为垂直下降的）。 从设计的角度考虑，由于窗函数设计法都是通过已有的窗函数对理想滤波器的改造，因此，可以用手算的办法方便的设计滤波器。而等波纹滤波器，其实现是通过大量的迭代运算来实现，这样的方法一般只能通过软件来设计。 项数的问题由于等波纹滤波器能较平均的分布误差，因此对于相同的阻带衰减，其所需的滤波系数比窗函数的要少。5.3 几点说明1.相频特性曲线形状不同说明
http://pic002.cnblogs.com/images/2011/314385/2011112410325263.jpg
上面第一个图是用角度为单位画出来的，下面的图是用rad单位画出来的。从图形可以观察到在0.3到0.8数字频率间两个图都是严格的线性相位，至于下面的图为什么在这个区间会有跳变是因为rad的区间只有-pi——pi，当相位由-pi继续增加时只能跳到pi而不能大于pi，而角度表示则可以连续增大。
2.调用firl或者reme函数时，用scale（缺省方式）对滤波器进行归一化，即滤波器通带中心频率处的响应幅值为0db。用noscale不对滤波器归一化。
好贴 顶顶顶
huangzj121 发表于 2014-3-6 08:26 static/image/common/back.gif
IIR的怎么搞啊 指导一下 补充资料
信号处理——Matlab使用指南
顶 请问：ADC采样的16位AD编码(0~0xffff)无符号数据，如何代入fir滤波函数里面进行计算？ 都是浮点的。 gongjian32 发表于 2014-3-6 14:34 static/image/common/back.gif
请问：ADC采样的16位AD编码(0~0xffff)无符号数据，如何代入fir滤波函数里面进行计算？ 都是浮点的。 ...
直接送进去算不行吗 ？
必须要送进去浮点吗 那就除个1.0
好贴！ 写的太好了 颓园花落 发表于 2014-6-8 18:52 static/image/common/back.gif
写的太好了
资料来自收集分享 供大家学习
谢谢支持 给力～顶 写的太好了写的太好了 滤波器设计 好资料 楼主，我想从fdatool导出十进制整数系数，量化完了还是小数呀？具体怎么设置呢？谢谢～
页:
[1]
2