基于Matlab设计FIR滤波器

2万 · 2014-3-4 19:51

一、摘要　　除了采用编程的方法实现滤波器之外，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的FDATool

MATLAB集成了一套功能强大的滤波器设计工具FDATool（Filter Design & Analysis Tool），可以完成多种滤波器的设计、分析和性能评估。

单击MATLAB主窗口下方的“Start”按钮，如图B.1所示，选择菜单“ToolBox” →“Filter Design” →“Filter Design & Analysis Tool（FDATool）”命令，打开FDATool，如图B.2所示。

　　　　　　　　　　　　　　　　　　　　图B.1 FDATool的启动

　　　　　　　　　　　　　　　　　　　　图B.2 FDATool的主界面

另外，在MATLAB主命令窗口内键入“fdatool”，同样可打开FDATool程序界面。

3.3 选择Design Filter

FDATool界面左下侧排列了一组工具按钮，其功能分别如下所述：

● 创建多速率滤波器（Create a Multirate Filter）

● 滤波器转换（TransForm Filter）

● 设置量化参数（Set Quantization Parameters）

● 实现模型（Realize Model）

● 零极点编辑器（Pole-zero Editor）

● 导入滤波器（Import Filter）

● 设计滤波器（Design Filter）

选择其中的按钮，进入设计滤波器界面，进行下列选择，如图B.3所示。

　　　　　　　　　　　　　　　　　　图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）。

　　　　　　　　　　　　图B.4 FIR滤波器幅频响应

在图的左侧列出了当前滤波器的相关信息：

● 滤波器类型为Direct Form FIR（直接I型FIR滤波器）

● 滤波器阶数为15

选择菜单“Analysis”→“Phase Response”，启动相频响应分析，如图B.5所示。由该图可以看到设计的FIR滤波器在通带内其相位响应为线性的，即该滤波器是一个线性相位的滤波器。

　　　　　　　　　　　　　　图B.5 滤波器相频响应

图B.6显示了滤波器幅频特性与相频特性的比较，这可以通过菜单“Analysis”→“Magnitude and Phase Response”来启动分析。

　　　　　　　　　　　　　　图B.6 滤波器幅频和相频响应

选择菜单“Analysis”→“Group Delay Response”，启动群时延分析。

FDATool还提供了以下几种分析工具：

● 群时延响应分析。

● 冲激响应分析（Impulse Response），如图B.7所示。

● 阶跃响应分析（Step Response），如图B.8所示。

● 零极点图分析（Pole/Zero Plot），如图B.9所示。

　　　　　　　　　　　　　　图B.7 冲激响应

　　　　　　　　　　　　　　图B.8 阶跃响应

　　　　　　　　　　　　　　图B.9 零极点图

求出的FIR滤波器的系数可以通过选择菜单“Analysis”→“Filter Coefficients”来观察。如图B.10所示，图中列出了FDATool计算的15阶直接I型FIR滤波器的部分系数。

　　　　　　　　　　　　　　图B.10 滤波器系数

3.5 量化

可以看到，FDATool计算出的值是一个有符号的小数，如果建立的FIR滤波器模型需要一个整数作为滤波器系数，就必须进行量化，并对得到的系数进行归一化。为此，单击FDATool左下侧的工具按钮进行量化参数设置。量化参数有三种方式：双精度、单精度和定点。在使用定点量化前，必须确保MATLAB中已经安装定点工具箱并有相应的授权。

3.6 导出滤波器系数

为导出设计好的滤波器系数，选择FDATool菜单的“File”→“Export”命令，打开Export（导出）对话框，如图B.11所示。

图B.11 滤波器系数Export对话框

在该窗口中，选择导出到工作区（Workplace）。这时滤波器系数就存入到一个一维变量Num中了。不过这时Num中的元素是以小数形式出现的：

Num=

Columns 1 through 9

-0.0369 0.0109 0.0558 0.0054 -0.0873 -0.0484 0.1805 0.4133 0.4133

Columns 10 through 16

0.1805 -0.0484 -0.0873 0.0054 0.0558 0.0109 -0.0369

由此，可以得到低通滤波器的系数。

资料来自网络仅供大家联系使用

2万 · 2014-3-4 19:52

往往到最后一步滤波器系数就存入到一个一维变量Num大家可能就不知道到哪儿去找到系数
很简答就是在Command Window中输入Num 回车就显示啦

2万 · 2014-3-4 19:52

接下来给大家一个历程虽然写的是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[n] = sin(nΩ1)/n/pi = sin(0.64pi*n)/n/pi
窗函数为：
w[n] = 0.42 - 0.5cos(2pi*n/24) + 0.8cos(4pi*n/24)
则滤波器脉冲响应为：
h[n] = h1*w[n] |n|<=12
h[n] = 0             |n|>12
3.3  滤波器的差分方程根据滤波器的脉冲响应计算出h[n]，然后将脉冲响应值移位为因果序列
这里计算h[n]的值，采用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.018

Columns 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[n] = 0.001x[n-2] - 0.002x[n-3] - 0.002x[n-4] + 0.01x[n-5] - 0.009x[n-6] - 0.018[n-7]

　　　　+ 0.049x[n-8] -0.02x[n-9] - 0.11x[n-10] + 0.28x[n-11] + 0.64x[n-12] + 0.28x[n-13]

　　　　- 0.11[n-14] - 0.02x[n-15] + 0.049x[n-16] - 0.018x[n-17] - 0.009x[n-18] + 0.1x[n-19]

　　　　- 0.002x[n-20] - 0.002x[n-21] + 0.001x[n-22]

3.4 DSP实现

根据滤波器系数，编写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 SAMPLEF  10000
#define PI 3.1415926

float InputWave();
float FIR();

float fHn[FIRNUMBER]={ 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[FIRNUMBER]={ 0.0 };
float fInput,fOutput;
float fSignal1,fSignal2;
float fStepSignal1,fStepSignal2;
float f2PI;
int i;
float fIn[256],fOut[256];
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[nIn]=fInput;
  nIn++; nIn%=256;
  fOutput=FIR();
  fOut[nOut]=fOutput;
  nOut++;
  if ( nOut>=256 )
  {
nOut=0;  /* 请在此句上设置软件断点 */
  }
}
}

float InputWave()
{
for ( i=FIRNUMBER-1;i>0;i-- )
  fXn=fXn[i-1];
fXn[0]=sin(fSignal1)+cos(fSignal2)/6.0;
fSignal1+=fStepSignal1;
if ( fSignal1>=f2PI ) fSignal1-=f2PI;
fSignal2+=fStepSignal2;
if ( fSignal2>=f2PI ) fSignal2-=f2PI;
return(fXn[0]);
}

float FIR()
{
float fSum;
fSum=0;
for ( i=0;i<FIRNUMBER;i++ )
{
  fSum+=(fXn*fHn);
}
return(fSum);
}
3.5  仿真结果

　　左上角的波形，为混叠有高频噪声的输入波形；右上角的是该波形的频谱。
　　左下角的波形，为经过FIR低通滤波器之后的输出波形；右下角的是该波形的频谱。
由实验结果可知，FIR滤波器能起到很好的滤波效果。

2万 · 2014-3-4 19:55

一、摘要　　前面一篇**介绍了通过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和β值。

第二步：通过程序设计滤波器。

程序如下：

b = fir1(29,0.4,kaiser(30,4.55));

[h1,w1]=freqz(b,1);

plot(w1/pi,20*log10(abs(h1)));

axis([0,1,-80,10]);

grid;

xlabel('归一化频率/p') ;

ylabel('幅度/dB') ;

波形如下：

方法二：

采用[n,Wn,beta,ftype] = kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器。

　　这里的函数kaiserord(f,a,dev)或者kaiserord(f,a,dev,fs)：

　　f为对应的频率，fs为采样频率；当f用数字频率表示时，fs则不需要写。

　　a=[1 0]为由f指定的各个频带上的幅值向量，一般只有0和1表示；a和f长度关系为（2*a的长度）- 2=（f的长度）

　　devs=[0.05 10^(-2.5)]用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等，计算公式：

阻带衰减误差=αs，通带衰减误差=αp，可有滤波器指标中的3、4条得到。

　　fs缺省为2Hz。

程序如下：

fcuts = [0.3 0.5]; %归一化频率omega/pi，这里指通带截止频率、阻带起始频率

mags = [1 0];

devs = [0.05 10^(-2.5)];

[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs); %计算出凯塞窗N，beta的值

hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

freqz(hh);

波形如下：

　　实际中，一般调用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=[800 1000];

mval=[1 0];

dev=[0.0559 0.01];

fs=4000;

[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);

b=remez(N,fpts,mag,wt);

[h,w]=freqz(b,1,256);

plot(w*2000/pi,20*log10(abs(h)));

grid;

xlabel('频率/Hz') ;

ylabel('幅度/dB');

波形如下：

例3 利用MATLAB编程设计一个数字带通滤波器，指标要求如下：通带边缘频率：Ωp1=0.45pi，Ωp2=0.65pi，通带峰值起伏：δ1<=1[dB]。阻带边缘频率：Ωs1=0.3pi，Ωs2=0.8pi，最小阻带衰减：δ2>=40[dB] 。方法一：窗函数法

程序如下：

[n,wn,bta,ftype]=kaiserord([0.3 0.45 0.65 0.8],[0 1 0],[0.01 0.1087 0.01]);%用kaiserord函数估计出滤波器阶数n和beta参数

h1=fir1(n,wn,ftype,kaiser(n+1,bta),'noscale');

[hh1,w1]=freqz(h1,1,256);

figure(1)

subplot(2,1,1)

plot(w1/pi,20*log10(abs(hh1)))

grid

xlabel('归一化频率w');ylabel('幅度/db');

subplot(2,1,2)

plot(w1/pi,angle(hh1))

grid

xlabel('归一化频率w');ylabel('相位/rad');

波形如下：

滤波器系数为：

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)，得幅频特性和相频特性曲线：

方法二：等波纹法设计
程序如下：
[n,fpts,mag,wt]=remezord([0.3 0.45 0.65 0.8],[0 1 0],[0.01 0.1087 0.01]);%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w，使remez函数设计出的滤波器满足f、a及dev指定的性能要求。
h2=remez(n,fpts,mag,wt);%设计出等波纹滤波器
[hh2,w2]=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

波形如下：

滤波器系数如下：

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);得幅频特性和相频特性曲线：

方法三：采用FDATool工具 这种方法需要事先计算出滤波器的阶数，bate值，然后设置相应参数，最后生成滤波器。
设置界面如下图所示：

　　将上述圈圈的区域设置好之后，生成滤波器，最后通过analysis菜单可以观察生成的滤波器的各种特性曲线和滤波器系数。这里的滤波器系数跟方法一的一样。

波形如下：

五、结果分析5.1 滤波器设计总结

　　FIR滤波器实现一般采用窗函数法和等纹波设计法。窗函数法还包含两个分支，一种是用公式先手动算出N值和其他对应得窗函数参数值，再代入窗函数和fir1实现，一种是用函数*rord估算出N和相应参数再用fir1实现。不过要注意*rord会低估或高估阶次n，可能会使滤波器达不到指定的性能，这时应稍微增加或降低阶次。如果截止频率在0或Nyquist频率附近，或者设定的dev值较大，则得不到正确结果。

　　滤波器实现形式及特点：由于一般的滤波器在利用窗函数是其通带波纹和阻带波纹不同（一般为第一个阻带波纹最大）因此，在满足第一个阻带衰减旁瓣时，比其频率高的旁瓣，它们的衰减都大大超出要求。而根据阻带衰减与项数的近似关系N = P(δ2)*fs/TW，可得当阻带衰减越大，所需项数越多。

5.2 窗函数法和等波纹设计的不同之处　　窗函数设计是通过最小平方积分办法来设计的，即该滤波器的误差为：
　　

　　即要求最小方法来设计滤波器，这样的滤波器更忠实于理想滤波器（即滤波系数更接近于理想滤波器）。

证明如下：

因此，幅度频谱差值越小，实际滤波器就越接近理想滤波器。

　　而等波纹滤波器是通过最大加权误差最小化来实现，其误差为：

要求该误差最小来实现滤波器，得出来的滤波系数较窗函数设计相差较远。

以下通过对例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,[0:30],(abs(delta1).^2))/31;

y2=sigsum(0,16,[0:16],(abs(delta2).^2))/17;

结果如下：

y1 =

1.9099e-004

y2 =

0.0278

　　由此得到用窗函数实现的滤波系数比用等波纹滤波器系数的每一项更接近于理想滤波器（y1为用窗函数实现的与理想滤波器的差值，y2为用等波纹滤波器实现的与理想滤波器的差值）；

　　对比二者的幅度频谱可知，等波纹滤波器阻带边缘比用窗函数实现的更平滑（理想滤波器为垂直下降的）。

　　从设计的角度考虑，由于窗函数设计法都是通过已有的窗函数对理想滤波器的改造，因此，可以用手算的办法方便的设计滤波器。

而等波纹滤波器，其实现是通过大量的迭代运算来实现，这样的方法一般只能通过软件来设计。

　　项数的问题由于等波纹滤波器能较平均的分布误差，因此对于相同的阻带衰减，其所需的滤波系数比窗函数的要少。

5.3 几点说明1.相频特性曲线形状不同说明

　　上面第一个图是用角度为单位画出来的，下面的图是用rad单位画出来的。从图形可以观察到在0.3到0.8数字频率间两个图都是严格的线性相位，至于下面的图为什么在这个区间会有跳变是因为rad的区间只有-pi——pi，当相位由-pi继续增加时只能跳到pi而不能大于pi，而角度表示则可以连续增大。

2.调用firl或者reme函数时，用scale（缺省方式）对滤波器进行归一化，即滤波器通带中心频率处的响应幅值为0db。用noscale不对滤波器归一化。