DSP课程设计基于DSP的实时FIR滤波器的实现.doc

设计题目基于DSP的实时FIR滤波器的实现设计要求1.熟悉AD50/McBSP的初始化设置及其使用方法，以及实时R滤波器的实现2.了解DSP/BIOS程序设计。3.熟悉DSPLIB中FIR函数的调用及参数设置。4.掌握FIR滤波器的原理与设计方法；5.学会在CCS平台编写C程序实现FIR滤波。设计过程1.设计原理1.1 FIR滤波器设计原理 在数字信号处理应用中往往需要设计线性相位的滤波器，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性。FIR滤波器不断地对输入样本x（n）延时后，再作乘法累加算法，将滤波结果y（n）输出，因此，FIR实际上是一种累加运算。在数字滤波器中，FIR滤波器的最主要的特点是没有反馈回路，故不存在不稳定的为题，同时，可以在幅度特性随意设置的同时，保证精确的线性相位。稳定和线性相位特性是FIR滤波器的突出优点。另外，它还有以下特点：设计方式是线性的；硬件容易实现；滤波器过渡过程具有有限区间；相对IIR滤波器而言，阶次较高，其延迟也要比同样性能的IIR滤波器大得多。FIR数字滤波器系统的传递函数为：通过反z变换，数字滤波器的差分方程为： 由上式可以得出如下图所示的直接型结构，这种结构又可以成为卷积型结构。将转置理论应用于图1.1可以得到转置直接型结构。将式中的系统函数H（z）分解成若干一阶和二阶多项式的连乘积： 则由此式可以构成如图1.2所示的级联型结构。其中为一阶节，为二阶节。每个一阶节、二阶节可用图1.2级联型结构实现。当M1=M2时，即得到图1.3所示的具体结构。这种结构的每一节都便于控制零点，在需要控制传输零点时可以采用。但是它所需要的系数a比直接型的h（n）多。在对滤波器计算时间没有特殊要求的时候可以采用这种形式。若需要严格考虑滤波器的计算时间则需要同时考虑它们的优点及缺点来设计。这在算法设计的时候要使用软件编辑环境来计算运行的时间问题。通常FIR的计算时间都较长。很多时候我们需要牺牲时间来获得想要得到的滤波区功能。 图1.1 FIR滤波器直接型结构图 图1.2 级联型结构图 图1.3 级联型具体结构 FIR滤波器实质上就是一个分节的延迟线，把每一节的输出用滤波器系数进行加权累加，便得到滤波器的输出结果，它总是稳定并且可实现的。在一些工程实际应用（如：图像处理、数据调制解调）中，往往对相位要求较高。FIR滤波器可以实现严格的线性相位，从而得到了广泛应用。它的差分方程数学表达式为：式中，N是FIR 滤波器的抽头系数，x（n）表示在n时刻输入的信号样值，h（n）表示滤波器的第n级抽头系数。横截型FIR滤波器的结构如图1.4所示： X(n)y(n)h(0)h(1)h(2)h(N-2)h(N-1)z -1z -1z-1 图1.4 FIR滤波器的横截型结构2设计方案2.1 方案一: 窗函数法窗函数法也称为傅里叶级数法。理想的数字滤波器频率特性是无法实现的，FIR的设计就是要寻找一个可以得到的频率特性来逼近，这相当于用一个可实现的单位脉冲响应 h（n）去逼近一个理想单位脉冲响应。可由理想频率特性通过傅氏反变换得到，。一般来说，这样得到的理想单位脉冲响应序列是个无限长序列，因而是非因果的。设有一个截止频率为的理想线性相位低通，延时为，其频率特性是：得到： 这是一个以n=为中心偶对称的无限长非因果序列，要想用一个有限长的因果序列去逼近它，最简单的方法是截取从0到N-1的一段来表示它，及=;其他N：。同时，为了保证线性相位，还要满足偶对称。这就像透过一个窗口观看到的一段，因此就表示成和一个窗口函数的乘积，这样对的求解就变为，这里的就被称为窗函数，既然一个频域上的标准的矩形窗口对应于时域是一个无限长的序列，那么在时域上截取一段必然会造成频域的矩形窗口的失真。相应的，截取出的信号也相应失真，从而需要改变原来窗口的形状来修正经过时域截取后的窗口失真。2.2 方案二: 频率抽样法窗函数法是从时域出发，用窗函数截取理想的得到，以此有限长的近似，这样得到的频率响应逼近于理想的频率响应。频率抽样法是从频率出发，将给定的理想频率响应加以等间隔抽样然后以此作为FIR滤波器的频率响应抽样值，再根据DFT（离散傅氏变换）定义由频域这N个抽样值来唯一定义一个有限长序列，同样也可以算出FIR滤波器的系统抽样函数及频率响应，可以推出频率响应是频率抽样值与线性相位因子及如下内插函数的线性组合。所以，在各频率取样点上，实际滤波器的频率响应是严格地和所要求的滤波器的频率响应一致的，逼近误差为零，但在抽样点之间的频率响应是各取样点的内插函数的延伸叠加而成，有一定的逼近误差，误差大小取决于频率响应曲线的圆滑程度和抽样点的密度。为了减少误差就要增加抽样点数目即增大采样频率，抽样点之间的理想频率特性变化越陡，则逼近误差越大，在理想频率特性的不连续点附近会产生肩峰和纹波。频率抽样法的优点是可以在频域设计，适于利用最优化方法，而且这种方法特别适用于窄带选频滤波器，但频率抽样法的抽样频率只能是2/N的整数倍或2/N的整数倍加上/N。不能保证截止频率的准确取值，要实现精确的就必须取N足够大，相应的计算量也很大。此外，它的阻带最大衰减一般，也只有30-50dB左右，很难满足频域特性要求较高的场合。2.3 方案三: Chebyshev逼近法窗函数法和频率采样法设计的滤波器的频率特性都是在不同意义上对所给理想频率特性的逼近。由数值逼近理论可知，对某个函数f（x）的逼近一般有以下三种方法：插值法（Interpolating Way）最小平方逼近法（Least Square Approaching Way）一致逼近法（Consistent Approaching Way）切比雪夫最佳一致逼近的基本思想是，对于给定区间a,b上的连续函数f(x)，在所有n次多项式的集合中，寻找一个多项式p(x)，使它在a,b上对f(x)的偏差和其它一切属于的多项式p(x)对f(x)的偏差相比是最小的，即：切比雪夫逼近理论，解决了p（x）的存在性、唯一性和如何构造等问题。但切比雪夫逼近法的效率初始极值频率点的估计，且通带和阻带内波纹数较多，这是Chebyshev方法的两个主要缺点。3设计步骤3.1 FIR窗函数法的基本步骤：（1）把经傅里叶变换，得；（2）对自然截短到所需的长度；（3）将截短后的右移M个采样间隔，得；（4）将乘以合适的窗函数得到所要滤波器的冲击响应，窗函数以n=M对称。利用所求得的单位抽样响应，即可用硬件构成滤波器的转移函数 H（z），也可利用在计算机上用软件实现滤波。3.2 典型传函数介绍数据窗在FIR滤波器的窗函数设计中起着重要的作用，它的性能的好坏直接影响着滤波器的过渡带宽和衰减的大小。对窗函数总的要求，是希望他的频谱中的主瓣尽量窄，旁瓣幅度尽量小，尤其是第一旁瓣。此外，窗函数还应该满足下列要求以便可以定量地比较各函数的性能。W（n）应是非负的实偶函数，为了使滤波器获得较大的主旁瓣能量化，从对称中心开始w（n）应是非递增的；为了保证滤波器的通带增益为1，应有：为了保证滤波器的相位特性不因加窗而改变，一般要求是恒正的；这里给出如下三个频域指标作为窗函数性能的性能参数：（1）3dB带宽B，它是主瓣归一化幅度 下降到-3dB时的带宽。当数据长度为N时，最大可能的频率分辨率是，则B的单位可以是；（2）最大旁瓣峰值A（dB）。A越小，由旁瓣引起的振荡幅度越小；（3）旁瓣峰值渐进衰减速度D（dB/oct）；3.2.1 矩形窗（Rectangle Window）由于矩形窗时讲理想的通带所对应的时域部分，无限长截断到有限长。,n=0,1,N-1它的波形及其幅频响应曲线如3.1示： 图3.1 矩形窗波形及其幅频响应曲线（N=21）它的窗口形式，这样在转换成频域是，理想的矩形通带会产生畸变，畸变情况如下图3.3所示 图3.3 经矩形窗截短后的畸变由上图可见，原本十分理想的矩形窗的传递函数现在变成了目前的形状。3.2.2 三角窗（Triangular） 它的波形及其幅频响应曲线如3.4所示： 图3.4 三角窗的波形及其幅频响应曲线（N=21）它的窗的形式为： 幅频特性为3.2.4 汉宁窗（Hanning Window） 它的波形及其幅频响应曲线如图3.5所示： 图3.5 汉宁窗的波形及其幅频响应曲线（N=21）由得到的频谱图，可以发现旁瓣大大的降低了，从而使能量集中在主瓣中，但是主瓣的宽度增大了一倍。它的窗的形式是：利用傅氏变换的调制特性，可以得到矩形窗的幅频特性表示的升余弦窗的幅频特性：3.2.5 海明窗（Hamming Window） 海明窗（hamming）也称为改进升余弦窗，若对海明窗再做一些调整就可以得到使旁瓣最小的窗函数：同时它的幅频特性为：它的波形及其幅频响应曲线如3.6所示： 图3.6 海明窗的波形及其幅频响应曲线（N=21）这将使99.96%的能量集中在主瓣内。这样，在与升余弦窗相等的主瓣宽度下获得了更好的旁瓣抑制效果。3.2.6 布莱克曼窗（Blackman Window）它的窗的形式为： 幅频响应为：它的波形及其幅频响应曲线如下图3.7所示： 图3.7 布莱克曼窗的波形及其幅频响应取下（N=21） 比较各中窗函数可知，矩形窗函数具有最窄的主瓣B，但也有最大的旁瓣峰值A和最慢的衰减速度D。汉宁窗的主瓣宽度稍宽，但有着较小的旁瓣和较大的衰减速度，因而是窗函数设计法中比较好的选择。将上述各种窗的指标列于下表3.1中名称近似过渡带宽精确过渡带宽最小阻带衰减矩形4/N1.8/N21dB三角8/N6.4/N25 dB汉宁8/N6.2/N44 dB海8/N6.6/N51 dB布莱克曼12/N11/N74 dB3.3.1 DSP中的McBsp接口电路基本介绍McBSP通过个引脚(DX、DR、CLKX、CLKR、FSX、FSR和CLKS)与外设接口。DX和DR引脚完成与外部设备进行通信时数据的发送和接收，由CLKX、CLKR、FSX、FSR实现时钟和帧同步的控制。由CLKS来提供系统时钟。发送数据时，CPU和DMA控制器将要发送的数据写到数据发送寄存器DXR，在FSX和CLKX作用下，由DX引脚输出。接收数据时，来自DR引脚的数据在FSR和CLKR作用下，从数据接收寄存器DRR中读出数据。接收和发送帧同步脉冲即可以由内部采样速率产生器产生，也可以由外部脉冲源驱动，McBSP分别在相应时钟的上升沿和下降沿进行数据检测。 串行口的操作由串行口控制寄存器2SPCR和引脚控制寄存器PCR来决定；接收控制寄存器RCR和发送控制寄存器XCR分别设置接收和发送的各种参数，如帧长度等。3.3.2 AD50内部结构介绍AD50是TI公司生产的一个16位、音频范围（采样频率为2K22.05KHZ）、内含抗混叠滤波器和重构滤波器的模拟接口芯片，它有一个能与许多DSP芯片相连的同步串行通信接口。AD50C片内还包括一个定时器（调整采样率和帧同步延时）和控制器（调整编程放大增益，锁相环PLL，主从模式）。AD50有28脚的塑料SOP封装（带DW后缀）和48脚的塑料扁平封装（带PT后缀），体积较小，适应于便携设备。AD50的工作温度范围是070，单一5V电源供电或5V和3.3V联合供电，工作时的最大功耗为120 mW。AD50的内部结构简图如下：最上面第一通道为模拟信号输入监控通道，第二通道为模拟信号转化为数字信号（A/D）通道，第三通道为数字信号转化为模拟信号（D/A）通道，最下面一路是AD50的工作频率和采样频率控制通道。 图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：保留3.3.2.2 AD50与DSP的接口电路AD50与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语音模块控制及原理图3.4 MCBSP初始化步骤：（1）设SPCR（12）中XRSTRRSTFRST0，使MCBSP口处于复位状态；（2）设置各个相关寄存器对应位；（3）等待两个数据周期，确保MCBSP内部同步；（4）向发送数据寄存器DXR写入数据；（5）设SPCR（12）中XRSTRRSTFRST1，使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 fc_clock port2000#define LowPass 0#define HighPass 1#define BandPass 2ioport unsigned fc_clock;ioport unsigned temp;/McBSP2初始化表const unsigned short init_tblx2=0,0x0000,1,0x0200,/上两目使McBSP1处于复位状态2,0x0040,3,0x0000,4,0x0040,5,0x0000,6,0x0101,7,0x2000,8,0x0000,9,0x0000,10,0x0000,11,0x0000,12,0x0000,13,0x0000,14,0x000c,/下两目使McBSP 开始工作0,0x0001,1,0x0043;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+);void initMCBSP()int i;for(i=0;in */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=a) HDi=WC; else b=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); else Wi=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=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*cos(float)i*wf); im=im+Hi*sin(float)i*wf); d=sqrt(pow(re,2)+pow(im,2); /*求模(幅度)*/ DBk=20.0*log10(d); /*转换为对数表示形式*/ /*对心电图进行