DSP的FIR设计(低通滤波)C语言编写

一、设计目的

低通滤波器设计。

本设计中使用的信号为

信息信号： signal=sin(2*pi*sl*n*T)

高频噪声1：noise1=0.7*sin(2*pi*ns1*n*T) 高频噪声2：noise2=0.4*sin(2*pi*ns2*n*T) 混合信号： x=(signal+noise1+noise2)

其中sl=500Hz ，ns1=3000Hz ，ns2=8000Hz ，T=1/20000。混合信号波形为滤波器输入信号波形，信息信号波形为输出信号波形，滤波器的效果为滤除两个高频噪声。

二、FIR 滤波器基本理论

（1）FIR 滤波器的特点

数字滤波器的功能，就是把输入序列通过一定的运算变换成输出序列。它的实现方法有很多，其中比较常用到的是无限长脉冲响应滤波器 IIR 和有限长脉冲响应滤波器FIR 两种。

在计算量相等的情况下，IIR 数字滤波器比FIR 滤波器的幅频特性优越，频率选择性也好。但是，它有着致命的缺点，其相位特性不好控制。它的相位特性

)argH( )f(ω

ωj e

=是使频率产生严重的非线性的原因。但是在图像处理、数据传

输等波形传递系统中都越来越多的要求信道具有线性的相位特性。在这方面 FIR 滤波器具有它独特的优点，设FIR 滤波器单位脉冲响应h(n)长度为N ，其系统函数H(z)为

∑-=-=

1

)()(N n n

z

n h z H

H(z)是1

-z 的(N-1)次多项式，它在z 平面上有(N-1)个零点，原点z=0是(N-1)阶重极点。因此，H(z)永远稳定，它可以在幅度特性随意设计的同时，保证精确、严格的线性相位。

（2）FIR 滤波器的基本结构

数字滤波是将输入的信号序列，按规定的算法进行处理，从而得到所期望的输出序列，FIR 滤波器的差分方程为：

∑-=-=

1

)()(N k k

k n x a

n y

对上式进行Z 变换得到FIR 滤波器的传递函数为：

()()

()

∑-=-=

=

1

N i k

k

z

b

z X z Y z H

由上式可以看出，H(z)是1-z 的N-1次多项式，它在z 平面内有N-1个零点，同时在原点处有N-1个重极点。N 阶滤波器通常采用N 个延迟单元、N 个加法器与N+1个乘法器，取图中(a)、(b)两种结构。

图 FIR 滤波器的一般结构

因为FIR 滤波器的单位抽样响应是有限长的，所以它永远是稳定的。另外，若对 h(n)提出一些约束条件，那么可以很容易地使 H(z)具有线性相位，这在信号处理的很多领域是非常重要的。FIR 滤波器的设计任务，是要决定一个转移函数H(z)，使它的频率响应满足给定的要求。这里所说的要求，除了通带频率p ω、阻带频率及两个带上的最大和最小衰减p ?和s ?外，很重要的一条是保证H(z)具有线性相位。 （3）Chebyshev 逼近法

窗函数法和频率采样法设计出的滤波器的频率特性都是在不同意义上对所给理想频率特性()jw d e H 的逼近。由数值逼近理论可知，对某个函数f(x)的逼近一般有以下三种方法：

插值法(Interpolating Way)

最小平方逼近法(Least Square Approaching Way) 一致逼近法(Consistent Approaching Way)

切比雪夫最佳一致逼近的基本思想是，对于给定区间[a ，b]上的连续函数

()x f ，在所有n 次多项式的集合n ?中，寻找一个多项式 p(x)，使它在[a ，b]上

对()x f 的偏差和其它一切属于n ?的多项式 p(x)对f(x)的偏差相比是最小的，即

()()()()(){}x f x p x f x p

-=-max min ?max 切比雪夫逼近理论，这样的多项式是存在的，且是唯一的，并指出了构造这种最佳一致逼近多项式的方法，就是有名的“交错点组定理”。

切比雪夫逼近理论解决了p(x)的存在性、唯一性和如何构造等问题。J.H.McClellan 、T.W.Parks 、L.R.Rabiner 等人应用切比雪夫逼近理论提出了一

种设计FIR 滤波器的计算机辅助算法。这种算法由于是在一致意义上对()jw d e H 作最佳逼近，因而获得了较好的通带和阻带性能，并能准确地指定通带和阻带的边缘。但它的效率依赖于初始极值频率点的估计，且通带和阻带内波纹数较多，这是Chebyshev 方法的两个主要缺点。

三、FIR 滤波器的MA TLAB 实现

MATLAB 辅助DSP 实现FIR ,其总体过程为在DSP 中编写处理程序，在MATLAB 中利用滤波器设计、分析工具( FDATOOL) ,根据指定的滤波器性能快速设计一个FIR ,再把滤波器系数以头文件形式导入CCS 中,头文件中MATLAB 辅助DSP 实现FIR 数字滤波器含滤波器阶数和系数数组,在MATLAB 中调试、运行DSP 程序并显示、分析处理后的数据。使用该方法,便于采用汇编语言来实现程序。头文件名不变,当MATLAB 中设计的滤波器系数改变时,相应头文件中系数也改变,方便程序调试、仿真。 （1）输入信号的产生

首先利用Matlab 产生导入CCS 的dat 文件，具体实现如下代码所示 sl=500; %有效信号 ns1=3000; %高频噪声 ns2=8000; %高频噪声 fs=20000; %采样频率 N=1000; T=1/fs; n=0:N;

signal=sin(2*pi*sl*n*T);

noise1=0.7*sin(2*pi*ns1*n*T); noise2=0.4*sin(2*pi*ns2*n*T);

x=(signal+noise1+noise2);%待滤波信号 figure(1) plot(x) figure(2)

y=abs(fft(x));%待滤波频谱 df=n*(fs/N); plot(df,y) figure(3) plot(signal) figure(4)

ysignal=abs(fft(signal));%滤波后频谱 df=n*(fs/N); plot(df,ysignal) %滤波数据导出

xout=x/max(x); %归一化

xto_css=round(32767*xout);%数据取整 xoutcss=xto_css;

fid=fopen('input.dat','w'); %打开文件 fprintf(fid,'1651 1 0 0 0\n');%输出文件头 fprintf(fid,'%d\n',xoutcss); %输出 fclose(fid);

产生的时域波形如图所示：

020040060080010001200

-2

-1.5

-1

-0.5

0.5

1

1.5

2

图 输入信号波形

频谱如图所示：

00.20.40.60.81 1.2 1.4 1.6 1.8

2

x 10

4

50

100

150

200

250

300

350

400

450

500

图 输入频谱

经过滤波器后的预期时域波形如图：

020040060080010001200

-1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

1

图 输出时域波形

频谱如图所示

00.20.40.60.81 1.2 1.4 1.6 1.8

2

x 10

4

50

100

150

200

250

300

350

400

450

500

图 输出频谱

（2）滤波器的设计

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

a.打开Filter Design & Analysis Tool

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

图fadatool的启动

b.产生滤波器

阶数为81阶，这里应填80，比阶数少1。窗函数选择切比雪夫型（chebyshev），采样频率为20000Hz，通带截止频率为750Hz。

图滤波器的幅频特性

图滤波器的相位特性c.产生滤波器系数和头文件

图滤波器系数

图头文件的产生

最后将产生的头文件中的滤波器系数数组放入编写好的滤波器函数中，并在进行CCS中进行调试与测试。

四、FIR滤波器的Simulink仿真

通过Matlab的FDA TOOL设计滤波器后，为了确认该滤波器是否有效，现通过simulink进行建模仿真观察该滤波器是否有效。

把三个信号通过加法器相加后，再通过FDA TOOL设计的滤波器即可。注意此处设计的是数字滤波器，需要把三个输入信号数字化。

具体仿真模型建立如下：

分别观察滤波前后的时域波形图

图滤波前时域波形

图 滤波后时域波形

五、FIR 滤波器的DSP 实现

（1）DSP 中滤波器的1

-z 算法实现

FIR 滤波器的输出表达式为

1

/21

()()()()()

N N n n y n h i x n i h i s n i --===

-=-∑∑

式中，为滤波器系数；x(n)表示滤波器在n 时刻的输入；y(n)为n 时刻的输出。

它的基本算法是一种乘法-累加运算，即不断地输入样本x(n)，经过延时后，再进行乘法-累加，最后输出滤波结果y(n)。

1）线性缓冲区法：线性缓冲区法又称延迟线法，其特点： （a ）对于N 级的FIR 滤波器，在数据存储器中开辟一个N 单元的缓冲区（滑窗），用来存放最新的N 个输入样本；

（b ）从最老样本开始取数，每取一个样本后，将此样本向下移位； （c ）读完最后一个样本后，输入最新样本存入缓冲区的顶部。 2）循环缓冲区法：循环缓冲区法的特点如下：

（a ）对于N 级FIR 滤波器，在数据存储器中开辟一个N 单元的缓冲区（滑窗），用来存放最新的N 个输入样本；

（b ）从最新样本开始取数；

（c ）读完最后一个样本（最老样本）后，输入最新样本来代替最老样本，而其他数据位置不变；

（d ）用片内BK （循环缓冲区长度）寄存器对缓冲区进行间接寻址，使循环缓冲区地址首尾相邻。本次设计的FIR 滤波器所采用的就是循环缓冲区法。 （2）C 语言实现FIR

采用C语言算法在DSP平台上实现了FIR低通数字滤波器，C语言算法相比于汇编算法可移植性很强。这里是在TMS320VC5510DSP为平台编写的C语言算法, 此算法可以稍加改动用在其他DSP芯片上, 而汇编算法则不然。这种方法具有以下优点：

(a) 程序的入口和出口由C语言自动管理，不必手工编写汇编程序实现。

(b) 程序结构清晰，可读性强。

(c) 程序调试方便。由于C程序中的变量全部由C语言来定义，因此采用C 源码调试器可以方便地观察C语言变量。

(d) 可移植性较强，通用性较好。

具体代码如下：

#include"stdio.h"

#include"fdacoefs.h" //头文件包含滤波器的系数

#define N 81 //定义滤波器的阶数为81阶

#define Length 200 //定义缓冲区数组大小为200

long yn;

int input[Length]; //存放输入数据

int output[Length]; //存放输出数据

void main()

{

int m,n;

int *x;

for(n=0;n

{

x=&input[n]; //指针指向每次导入的数据

yn=0; //每做完一次乘累加后，把值赋给output数组后，从新归0

for(m=0;m

yn+=B[m]*(*(x++)); //做N次的乘累加

output[n]=yn; 把值赋给output数组

}

while(1); //做完滤波后使程序保持在本循环中

}

（3）CSS仿真调试

CCS 是TI 推出的用于开发其DSP 芯片的继承开发调试工具, 集编辑、编译、链接、软件仿真、硬件调试及实时跟踪等功能于一体, 极大地方便了DSP 程序的设计与开发, 此外还提供图形显示功能, 方便用户观察特定地址的波形。此外, 还需向工程中添加Link. cmd文件（源码见附录）。

在CCSV3.3中建立工程，把c源代码和.cmd文件导入后，外加rts. lib文件，它是TI提供的运行时支持库, 如果是C代码写的源程序, 必须要包含该库,该库

由TI公司做好放在CCS\ cgtools\ lib中, 源代码TI网站可以下载。添加完成后，编译通过后，下载.out文件，导入.dat文件后运行程序。

采用CCS 的图形显示功能, 分别观察输入信号x ( n)、输出信号y ( n)的时域波形和频域波形, 输入信号波形如图所示。

图滤波前时域波形

图滤波前频谱图

经过滤波后，观察输出波形如下

图滤波后时域波形

图滤波后频谱图

为了更加直观的观察滤波器的性能和滤波效果，把滤波前后的时域波形和频谱图进行对比和分析，具体如下：

图滤波前后对比时域和频谱

从上图中可以看出，输入信号明显有高频噪声，波形有很明显的失真。从输入频谱分析，可以看出除了有500Hz的信息信号还有3000Hz和8000Hz的高频噪声。

经过滤波器后，输出的时域波形有了明显改善，几乎接近信息信号。而频谱图中也可以看出3000Hz和8000Hz的高频噪声已经被滤除，而且对原来的信息

信号几乎没有影响。

综上分析，滤波器实现滤除高频噪声，而且保证原信息信号基本不失真的功能，滤波效果较好，性能优越。

（4）硬件仿真

通过CCS的软件仿真调试，发现滤波器性能符合要求，有较好的滤波效果。但是软件仿真过于理想，不过精确，因此在硬件平台上进行仿真调试，观察滤波器性能指标。

本仿真使用的C5509A硬件平台进行仿真，由于C语言的移植性较好，所以本滤波程序可以直接移植到C5509A平台上使用。首先创建工程，添加源文件、.cmd文件、csl5509a.lib和rts55.lib库文件，如下图所示;

图硬件工程

然后编译，导入数据和运行程序，观察输入输出波形如图：

件仿真效果差异不大，从而也可以发现采用C编写的程序通用性较好。