实时FIR滤波器的实现

1、赣南师院物理与电子信息学院课程设计报告姓名： 陈如意 邹丽蓉 贺玲玲班级：09电信本学号：090802003090802071090802053时间：2012年4月25日设计题目基于DSP的实时FIR滤波器的实现设计要求熟悉AD50/McBSP的初始化设置及其使用方法，以及实时R滤波器的实现了解DSP/BIOS程序设计。熟悉DSPLIB中FIR函数的调用及参数设置。掌握FIR滤波器的原理与设计方法；学会在CCS平台编写C程序实现FIR滤波。设计过程1.设计原理1.1 FIR滤波器设计原理在数字信号处理应用中往往需要设计线性相位的滤波器，FIR滤波器在保 证幅度特性满足技术要求的同时，很容易做到

2、严格的线性相位特性。FIR滤波 器不断地对输入样本x(n)延时后，再作乘法累加算法，将滤波结果y(n) 输出，因此，FIR实际上是一种累加运算。在数字滤波器中，FIR滤波器的最主要的特点是没有反馈回路，故不存在 不稳定的为题，同时，可以在幅度特性随意设置的同时，保证精确的线性相位。 稳定和线性相位特性是FIR滤波器的突出优点。另外，它还有以下特点：设计 方式是线性的；硬件容易实现；滤波器过渡过程具有有限区间；相对IIR滤波 器而言，阶次较高，其延迟也要比同样性能的IIR滤波器大得多。FIR数字滤波器系统的传递函数为：H ( z )=里=习顷z - n XCo通过反z变换，数字滤波器的差分方程为

3、：y (n)=祝 h x(n - k)k =0由上式可以得出如下图所示的直接型结构，这种结构又可以成为卷积型结 构。将转置理论应用于图1.1可以得到转置直接型结构。将式中的系统函数H(z)分解成若干一阶和二阶多项式的连乘积：H (z) = H H ik (z )EI H 2 k (z) k=1k=1则由此式可以构成如图1.2所示的级联型结构。其中H 1k(z) = O0D +%)z-1为一 阶节，H2k (z) = a02 + ak) + a2k为二阶节。每个一阶节、二阶节可用图1.2级 联型结构实现。当M1=M2时，即得到图1.3所示的具体结构。这种结构的每 一节都便于控制零点，在需要控制传

4、输零点时可以采用。但是它所需要的系数 a比直接型的h(n)多。在对滤波器计算时间没有特殊要求的时候可以采用这 种形式。若需要严格考虑滤波器的计算时间则需要同时考虑它们的优点及缺点 来设计。这在算法设计的时候要使用软件编辑环境来计算运行的时间问题。通 常FIR的计算时间都较长。很多时候我们需要牺牲时间来获得想要得到的滤波区功能。图1.1 FIR滤波器直接型结构图H H/r)e 为卜-4 取0)h w t - 图1.2级联型结构图,j)村腴h 卜 q广|图1.3级联型具体结构FIR滤波器实质上就是一个分节的延迟线，把每一节的输出用滤波器系数 进行加权累加，便得到滤波器的输出结果，它总是稳定并且可实

5、现的。在一些 工程实际应用(如：图像处理、数据调制解调)中，往往对相位要求较高。FIR 滤波器可以实现严格的线性相位，从而得到了广泛应用。它的差分方程数学表达式为：y(n) = N 1 b x(n - k)k=0式中，N是FIR滤波器的抽头系数，x(n)表示在n时刻输入的信号样 值，h(n)表示滤波器的第n级抽头系数。横截型FIR滤波器的结构如图1.4 所示：X(n)z-1z -1z -1ifUi fifh(0)h(1)h(2)h(N-2)h(N-1)y(n),勉图1.4 FIR滤波器的横截型结构2设计方案2.1方案一：窗函数法窗函数法也称为傅里叶级数法。理想的数字滤波器频率特性H(eJW)是

6、无 法实现的，FIR的设计就是要寻找一个可以得到的频率特性H(ejw) = tlh(n)e-网来逼近H(ejw)，这相当于用一个可实现的单位脉冲响应 n=0h (n)去逼近一个理想单位脉冲响应气(n)。七(n)可由理想频率特性H(ejw)通过傅氏反变换得到，h (n) = J H (ejw)d。一般来说，这样得到的理想单 d2兀d兀位脉冲响应序列七(n)是个无限长序列，因而是非因果的。设有一个截止频率 为的理想线性相位低通，延时为，其频率特性是：H ew)=ej%0 小 Hd I 0 兀tC得到:h (n)=，低(n-T sd兀(n-t )这是一个以n=t为中心偶对称的无限长非因果序列，要想用

7、一个有限长 的因果序列去逼近它，最简单的方法是截取从0到N-1的一段来表示它，及 h(n) = hd (n) (0 n N -1);其他 N： h(n) = 0。同时，为了保证线性相位，还要满足偶对称h(n) = h(N 1 n)。这就像透过一个窗口观看到的一段气(n)，因此h(n)就表示成七(n)和一个窗口函数的 乘积，这样对h(n)的求解就变为h(n)=气(n)*%,这里的W就被称为窗函数， 既然一个频域上的标准的矩形窗口对应于时域是一个无限长的序列，那么在时 域上截取一段必然会造成频域的矩形窗口的失真。相应的，截取出的信号也相 应失真，从而需要改变原来窗口的形状来修正经过时域截取后的窗口

8、失真。2.2方案二：频率抽样法窗函数法是从时域出发，用窗函数截取理想的h (n)得到h(n)，以此有限 d长的h(n)近似气(n)，这样得到的频率响应 H(ejw)逼近于理想的频率响应间隔抽样% (ejw)% (ejw)。频率抽样法是从频率出发，将给定的理想频率响应% (ejw)加以等=Hd (k)皿gN然后以此H (k )作为FIR滤波器的频率响应抽样值H (k )，再根据DFT (离散傅 d氏变换)定义由频域这N个抽样值来唯一定义一个有限长序列h(n)，同样也可以算出FIR滤波器的系统抽样函数H(z)及频率响应H(ejw)，可以推出频率响应H(ejw)是频率抽样值H(k)与线性相位因子e-

9、j低n一 1)/2及如下内插函数S(o ,k)的线一 sin 竺性组合。S(,k)=篇e-jN 一 2N , 冗k、sin(: 丁)2 N所以，在各频率取样点上，实际滤波器的频率响应是严格地和所要求的滤波器的频 率响应一致的，逼近误差为零，但在抽样点之间的频率响应是各取样点的内插函数的延 伸叠加而成，有一定的逼近误差，误差大小取决于频率响应曲线的圆滑程度和抽样点的 密度。为了减少误差就要增加抽样点数目即增大采样频率，抽样点之间的理想频率特性 变化越陡，则逼近误差越大，在理想频率特性的不连续点附近会产生肩峰和纹波。频率 抽样法的优点是可以在频域设计，适于利用最优化方法，而且这种方法特别适用于窄带

10、 选频滤波器，但频率抽样法的抽样频率只能是2n /N的整数倍或2n /N的整数倍加上n /N。不能保证截止频率。的准确取值，要实现精确的就必须取N足够大，相应的计 算量也很大。此外，它的阻带最大衰减一般，也只有30-50dB左右，很难满足频域特性 要求较高的场合。2.3方案三:Chebyshev逼近法窗函数法和频率采样法设计的滤波器的频率特性都是在不同意义上对所 给理想频率特性气(ejw)的逼近。由数值逼近理论可知，对某个函数f(x)的 逼近一般有以下三种方法：插值法(Interpolating Way)最小平方逼近法(Least Square Approaching Way)一致逼近法(C

11、onsistent Approaching Way)切比雪夫最佳一致逼近的基本思想是，对于给定区间a,b上的连续函数 f(x)，在所有n次多项式的集合中，寻找一个多项式p(x)，使它在a,b上对 f(x)的偏差和其它一切属于平的多项式p(xMWf(x)的偏差相比是最小的，即： nmax| p( x) - f (x) = minmax(p( x) - f (x)切比雪夫逼近理论，解决了 p (x)的存在性、唯一性和如何构造等问题。但 切比雪夫逼近法的效率初始极值频率点的估计，且通带和阻带内波纹数较多， 这是Chebyshev方法的两个主要缺点。FIR窗函数法的基本步骤：把气(ejw)经傅里叶变换

12、，得hd (n)；对匕(n)自然截短到所需的长度；将截短后的hd (n)右移M个采样间隔，得h(n)；将h(n)乘以合适的窗函数得到所要滤波器的冲击响应，窗函数以n=M对 称。利用所求得的单位抽样响应，即可用硬件构成滤波器的转移函数H(z)， 也可利用h(n)在计算机上用软件实现滤波。3.2典型窗函数介绍数据窗在FIR滤波器的窗函数设计中起着重要的作用，它的性能的好坏直 接影响着滤波器的过渡带宽和衰减的大小。对窗函数总的要求，是希望他的频 谱中的主瓣尽量窄，旁瓣幅度尽量小，尤其是第一旁瓣。此外，窗函数还应该 满足下列要求以便可以定量地比较各函数的性能。W(n)应是非负的实偶函 数，为了使滤波器

13、获得较大的主旁瓣能量化，从对称中心开始w(n)应是非1四递增的；为了保证滤波器的通带增益为1，应有：(0) = 土 j W(eje虬=1-兀为了保证滤波器的相位特性不因加窗而改变，一般要求W(ejw)是恒正的； 这里给出如下三个频域指标作为窗函数性能的性能参数：3dB带宽B，它是主瓣归一化幅度(20log览J)下降到-3dB时的、 、，, 一, ，一、，、2兀,、，带宽。当数据长度为N时，最大可能的频率分辨率是=竺，则B的单位N可以是宜；最大旁瓣峰值A(dB)。A越小，由旁瓣引起的振荡幅度越小；旁瓣峰值渐进衰减速度D(dB/oct)；矩形窗(Rectangle Window)由于矩形窗时讲理想

14、的通带所对应的时域部分，无限长截断到有限长。 (n) = 1,n=0,1,N-1它的波形及其幅频响应曲线如3.1示:Time domain5Eples 15208 6 4 2 0 o o o OeauT-DmANormalized Frequency (x兀 rad/sample)图3.1矩形窗波形及其幅频响应曲线（N=21）它的窗口形式Wn = 1*础，这样在转换成频域是，理想的矩形通带会产生畸变，畸 变情况如下图3.3所示40200-20-40-60-80Frequency domain0.20.40.60.8Normalized Frequency (x兀 rad/sample)10.8

15、Time domain0.25Eples 15206 4 o o pduMDmA图3.3经矩形窗截短后的畸变由上图可见，原本十分理想的矩形窗的传递函数现在变成了目前的形状。三角窗（Triangular）它的波形及其幅频响应曲线如3.4所示：-1000图3.4三角窗的波形及其幅频响应曲线（N=21）它的窗的形式为：w (n) = 1 n = -?,.-1,0,1,., J2.,2 .、幅频特性为(ej) = Ne _ j ( n sin( N /4)sin( / 2)3.2.4 汉宁窗(Hanning Window)10.8Time domain0.205Eples 15206 4 o o pd

16、uMDmANormalized Frequency Q兀 rad/sample)Da(pduTnQa它的波形及其幅频响应曲线如图3.5所示:图3.5汉宁窗的波形及其幅频响应曲线(N=21)由得到的频谱图，可以发现旁瓣大大的降低了，从而使能量集中在主瓣中， 但是主瓣的宽度增大了一倍。 1Rn (n)它的窗的形式是：(n) = -1-cos2利用傅氏变换的调制特性(ejdnx(n) o X(ej(-/)，可以得到矩形窗的幅频特性WR()表示的升余弦窗的幅频特性:W() = 0.5W ()+ 0.25 Wr -3.2.5 海明窗(Hamming Window)若对海明窗再做一些调整就可海明窗(ham

17、ming)也称为改进升余弦窗 以得到使旁瓣最小的窗函数：(n)=0.54 - 0.46cosRn (n)同时它的幅频特性为:W() = 0.54W () + 0.23 W (+ Wr +它的波形及其幅频响应曲线如3.6所示:Normalized Frequency S rad/sample)图3.6海明窗的波形及其幅频响应曲线（N=21）这将使99.96%的能量集中在主瓣内。这样，在与升余弦窗相等的主瓣宽度下 获得了更好的旁瓣抑制效果。3.2.6 布莱克曼窗（Blackman Window）它的窗的形式为：(n)=(2m )(4冗n )0.42 - 0.5cos+ 0.08 cos1N -1N

18、 -1)_Rn (n)幅频响应为:U(2兀)+ U3 +2兀、1+ 0.04U (3)+ U（4兀）3 + 一1 N J1N Jk N ；W ()=0.42U ()+ 0.23它的波形及其幅频响应曲线如下图3.7所示:/i/Time domain5Eples 15208 6 4 2 o o o O eauT-DmANormalized Frequency (x兀 rad/sample)图3.7布莱克曼窗的波形及其幅频响应取下（N=21）比较各中窗函数可知，矩形窗函数具有最窄的主瓣B,但也有最大的旁瓣 峰值A和最慢的衰减速度D。汉宁窗的主瓣宽度稍宽，但有着较小的旁瓣和较 大的衰减速度，因而是窗函

19、数设计法中比较好的选择。将上述各种窗的指标列 于下表3.1中名称近似过渡带宽精确过渡带宽最小阻带衰减矩形4n /N1.8n /N21dB三角8n /N6.4n /N25 dB汉宁8n /N6.2n /N44 dB海明8n /N6.6n /N51 dB布莱克曼12/N11n /N74 dBDSP中的McBsp接口电路基本介绍McBSP 通过 7 个引脚（DX、DR、CLKX、CLKR、FSX、FSR 和 CLKS）与外设接 口。DX和DR引脚完成与外部设备进行通信时数据的发送和接收，由CLKX、CLKR、 FSX、FSR实现时钟和帧同步的控制。由CLKS来提供系统时钟。发送数据时， CPU和DM

20、A控制器将要发送的数据写到数据发送寄存器DXR，在FSX和CLKX 作用下，由DX引脚输出。接收数据时，来自DR引脚的数据在FSR和CLKR作 用下，从数据接收寄存器DRR中读出数据。接收和发送帧同步脉冲即可以由内 部采样速率产生器产生，也可以由外部脉冲源驱动，McBSP分别在相应时钟的 上升沿和下降沿进行数据检测。串行口的操作由串行口控制寄存器2SPCR和引脚控制寄存器PCR来决定； 接收控制寄存器RCR和发送控制寄存器XCR分别设置接收和发送的各种参数， 如帧长度等。AD50内部结构介绍AD50是TI公司生产的一个16位、音频范围（采样频率为2K22.05KHZ）、内含抗混叠滤波器和重构滤

21、波器的模拟接口芯片，它有一个能与 许多DSP芯片相连的同步串行通信接口。AD50C片内还包括一个定时器（调 整采样率和帧同步延时）和控制器（调整编程放大增益，锁相环PLL，主从模 式）。AD50有28脚的塑料SOP封装（带DW后缀）和48脚的塑料扁平封装 （带PT后缀），体积较小，适应于便携设备。AD50的工作温度范围是070C， 单一 5V电源供电或5V和3.3V联合供电，工作时的最大功耗为120 mW。AD50的内部结构简图如下：最上面第一通道为模拟信号输入监控通道， 第二通道为模拟信号转化为数字信号（A/D）通道，第三通道为数字信号转化 为模拟信号（D/A）通道，最下面一路是AD50的工

22、作频率和采样频率控制通 道。图1.4 AD50的内部结构简图3.3.2.1 AD50内部寄存器及其作用AD50内部有7个数据和控制寄存器，用于编程控制它们的工作状态。寄存器0：空操作寄存器。寄存器1：软件复位软件掉电选择16位或15位工作方式硬件或软件二次通信请求方式的选择寄存器2：使能ALTDATA输入端为ADC选择16/15位方式寄存器3：选择FS与FSD之间延迟SCLK的个数告诉主机有几个从机被联上寄存器4：为输入和输出放大器选择放大器增益选择N来设置采样频率，fs=MCLK/ (128*N)或MCLK/(512*N) 在MCLK输入端使能外部时钟输入，并旁通内部的PLL 寄存器5， 6

23、：保留AD50与DSP的接口电路图1系统硬件接口音拟号入 SSAD50与TMS320VC5402是以SPI方式连接的。AD50工作在主机模 (M/S=1),提供SCLK (数据移位时钟)和FS (帧同步脉冲)。TMS320VC5402 工作于SPI方式的从机模式，BCLKX1和BFSX1为输入引脚，在接数据和发 数据时都是利用外界时钟和移位脉冲。3.3.2.3 AD50的通讯方式AD50有两个通信方式一种是15+1方式软件申请第二串行通信的；另一种 是用FC来切换通信方式。软件不太可靠而且15位精度小，对于音频数据处理 麻烦。故我们采用拉高FC以达到切换通信方式。普通AD/DA语音模块控制及原

24、理图3.4 MCBSP初始化步骤：设 SPCR (12)中 XRST=RRST=FRST=0,使 MCBSP 口处于复位状态;设置各个相关寄存器对应位；等待两个数据周期，确保MCBSP内部同步；向发送数据寄存器DXR写入数据；设 SPCR (12)中 XRST=RRST=FRST=1，使 MCBSP 口离开复位状态;信号处度饬率放为(6)等待两个数据周期，使MCBSP收发器启动。3.5软件设计3.5.1程序如下#include math.h#include stdlib.h#define pi 3.14159#include 5410reg.h#include stdio.h#define f

25、c_clock port2000#define LowPass 0#define HighPass 1#define BandPass 2 ioport unsigned fc_clock;ioport unsigned temp;/McBSP2初始化表const unsigned short init_tblx2=0,0 x0000,1,0 x0200,上两目使McBSP1处于复位状态2,0 x0040,3,0 x0000,4,0 x0040,5,0 x0000,6,0 x0101,7,0 x2000,8,0 x0000,9,0 x0000,10,0 x0000,11,0 x0000,12,

26、0 x0000,13,0 x0000,14,0 x000c,卜两目使McBSP开始工作0,0 x0001,1,0 x0043；short *init_tbl=(short*)init_tblx;int mtmp;int Type;void Delay()int temp,ft;for(ft=0;ft100;ft+)for (temp=0;temp1000;temp+);void ShortDelay()int tmp;for(tmp=0;tmp100;tmp+);void MainDelay(unsigned int count)int tmp;for(tmp=0;tmpcount;tmp+)

27、;void initMCBSP() int i;for(i=0;in */float WC;int n;int m;int l;float im,re;float a,b,p,wf,d;int k,i;/*/ main() for(m=1; m6; m+)/*参数初始化*/ for(i=0;i57;i+)HDi=0;Hi=0;Wi=0;/ Xi=0;for(i=0;i300;i+)DBi=0;n=21;/* n=57,为奇数 */WC=0.2; /* wc=0.10-0.90 */ l=300;a=(n-1)/2;for(i=0;in;i+) /*理想矩形函数(频域)的时域形式*/ if(i=

28、a)HDi=WC;elseb=i-a;HDi=sin(pi*b*WC)/(pi*b);switch(m) /*根据m的值选择不同的窗函数*/ /* m=1时用矩形窗*/case 1:for(i=0;in;i+) Wi=1.0;break;/* m=2时用三角窗(巴特莱特窗)*/ case 2:for(i=0;i=a)Wi=2.0-2.0*i/(n-1);elseWi=2.0*i/(n-1); break;/* m=3时用汉宁窗*/case 3:for(i=0;in;i+)Wi=0.5*(1.0-cos(2.0*pi*i/(n-1);break;/* m=4时用汉明窗*/case 4:for(i

29、=0;in;i+)Wi=0.54-0.46*cos(2.0*pi*(float)i/(n-1); break;/* m=5时用布莱克曼窗*/case 5:for(i=0;in;i+)Wi=0.42-0.5*cos(2.0*pi*i/(n-1)+0.08*cos(4.0*pi*i/(n-1); break;for(i=0;in;i+)Hi=HDi*Wi; /*对理想冲击响应加窗(时域相乘)*/p=pi/l;/*对加窗后的响应序列进行傅里叶变换，观察具幅频响应*/for(k=0;k=l-1;k+)wf=(pi*k)/l;re=0.0;im=0.0;for(i=0;in;i+)re=re+Hi*co

30、s(float)i*wf);im=im+Hi*sin(float)i*wf);d=sqrt(pow(re,2)+pow(im,2); /*求模(幅度)*/DBk=20.0*log10(d);/*转换为对数表示形式*/*对心电图进行滤波*/for(i=0;i57;i+)Yi=0;for(k=0;(k21&k!*-rJinraRe540-01Bnnlcpa 1 nt KcMgHr:praba peLntv ara ant at valid jinan-xlW邓1-SiGI 心相 Tb5* pIDOiM gL口E TO g I.41T *4114 I .皿 i Dt-5m peIm? points

31、/stv let flt 叫Lid Hdi?:-Ccd Canpoanr 勤占thai* pret-a psmu: to eKh1 i.d. LinrB-18 b I ft BJFX i1Faf shFlUJJ 夺t_. JU j二、二，.: ”ig. .M=1,矩形窗时滤波后心电图如下:M=2,三角窗时滤波后心电图如下:WxHtrji&J IJd CowpLcw ErTDFXP Wmm.L O RznarlcE-M=3,汉宁窗时滤波后心电图如下:UJu !_ L-fi . t r*T F. JJ |四 A 的 E 普B： t*v; g E 呼*： gk * 占3 A- o 口 EH(3 日画

32、担M=4,汉明窗时滤波后心电图如下:色 X 的 E 村：g fipw E*皿fjL1?la. I b b h E加A.Mmm fil4 CoipLrWn Q ErTDE-a:P .jtriiirfTEFrnJjrOfdE8 LUfcWe- J Sam =*xiBffVEDUy4QQgiaDLbLnySa口、-f Iwfai. LkfLftiAihhM.M=5,布莱克曼窗时滤波后心电图如下:imposer囚 File Edit View Project Debug Profiler GEL Option Tools DSP/BIOS Window HelpFIR. pjt愿1 |阳T |Pebu回座国国叵f FilesE-O GEL filesProjects白宙 FIR.pjt (Debug): |_| Dependent Projed? LJ DSP/BIOS Config: 口 Generated Files: |_J IndudeLibraries白-包 Sourcej j