基于Matlab的FIR数字滤波器问题的优化设计

基于Matlab的FIR数字滤波器问题的优化设计摘要相对于硬件方式,采用软件方式可以通过对滤器参数的改变来调整滤波器的性能以达到实现滤波目标,本文根据FIR滤波器的原理，运用MATLAB进行了惩罚函数设计法。关键词FIR数字滤波器FFTMATLAB引言数字信号处理中,滤波具有十分重要的地位与作用,FIR数字滤波器在通信、语音、图像、自动控制等众多领域得到了广泛的应用。它具有精度高、可靠性好、灵活性大等特点。FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域经过ADC采样之后，就变成了数字信号。MATLAB是一种面向科学和工程计算的语言,集数值分析、矩阵运算、信号处理和图形显示于一体,具有编程效率高、调试手段丰富、扩充能力强等特点。1FIR数字滤波器的原理数字滤波器原理一般具有如下差分方程：1100()()()kkNNkkynaxnkbynk若设其为冲激响应为h(z)，0≤n≤N-1，则其系统函数为：10()()NnnHzhnzFIR数字滤波器的设计实质是确定能满足要求的转移序列或脉冲响应的常数问题,设计方法主要有惩罚函数法、频率采样法和波纹最佳逼近法等。2惩罚函数法的设计步骤(1)给定所要求的频率响应Hd(ejw);(2)由Hd(ejw)的傅立叶反变换求出hd(n);(3)根据过渡带宽及阻带最小衰减的要求,选择窗的形状并估计窗口长度N;(4)计算滤波器的单位脉冲响应:()()()dhnhnwn3FFT概念理解及其功用FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N.由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。3.1FFT采样点数与精度关系说明在采样频率一定的情况下，采样点数越多，FFT的分辨率越高，所以需要调试FFT的分辨率和采样点数的关系。现在VC++中，做以下仿真运算：采样频率是45×32＝1440HZ，设计50HZ陷波器，对带有工频干扰信号的45HZ有用信号进行FIR滤波和FFT频谱特性分析。加载带有干扰的信号：S1(T)=2.0+1024.0*cos(M_PI*2*45*t+M_PI*1/3)+2048*cos(M_PI*2*50.1*t+M_PI*1/3);S2(T)=3.0+800.0*cos(M_PI*2*45*t+M_PI*1/4);注：下划线是有用信号、虚线是干扰信号；干扰信号的幅值是有用信号的2倍；有用信号的频率是45HZ，干扰信号频率是50.1HZ。假设S1是测试的电压信号、S2假设是测试的电流信号，求出S1、S2的在45HZ频点的相位和幅值即为所求的介质损耗的损耗角(没有求tan)和电容量。分别求出采样点数为1024、2048、4096,8192时，FFT的分辨率情况，如下表1：采样点数实际相差（度）仿真相位差（度）误差(%)实际电容量仿真电容值误差(%)10241514.991220-0.05850.0096060.009629+0.2394204814.993467-0.04360.009607+0.0104409615.019326+0.12880.009598-0.0833819215.004272+0.02850.009605-0.0104程序运行流程:通过S1,S2生成待测量数据；通过MARLAB设计50HZ陷波器，把49.5HZ～50.5HZ范围内的干扰彻底滤除掉以便数据下一步处理，使用FIR对数据进行滤波；对已经FIR滤波后的数据进行FFT运算。3.2分析结果如下：（1）如果信号不加干扰（如信号S3,S2），且没有经过FIR滤波，直接对生成的数据进行FFT运算后，分析45HZ频点的相频特性.（2）S3(T)=2.0+1024.0*cos(M_PI*2*45*t+M_PI*1/3);S2(T)=3.0+800.0*cos(M_PI*2*45*t+M_PI*1/4);4MATLAB环境下的设计实例假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30°、幅度为3V的交流信号，以及一个频率为75Hz、相位为90°、幅度为1.5V的交流信号。用数学表达式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)式中cos参数为弧度，所以-30°和90°要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第50个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？图1FFT结果从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：1点：512+0i2点：-2.6195E-14-1.4162E-13i3点：-2.8586E-14-1.1898E-13i50点：-6.2076E-13-2.1713E-12i51点：332.55-192i52点：-1.6707E-12-1.5241E-12i75点：-2.2199E-13-1.0076E-12i76点：3.4315E-12+192i77点：-3.0263E-14+7.5609E-13i很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，结果如下：1点：51251点：38476点：192按照公式，可以计算出直流分量为：512/N=512/256=2；50Hz信号的幅度为384/(N/2)=384/(256/2)=3；75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来的幅度是正确的。然后再来计算相位信息。直流信号没有相位可言，不用管它。先计算50Hz信号的相位，atan2(-192,332.55)=-0.5236,结果是弧度，换算为角度就是180*(-0.5236)/pi=-30.0001。再计算75Hz信号的相位，atan2(192,3.4315E-12)=1.5708弧度，换算成角度就180*1.5708/pi=90.0002。可见，相位也是对的。根据FFT结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。本测试数据使用的matlab程序如下：closeall;Adc=2;%直流分量幅度A1=3;%频率F1信号的幅度A2=1.5;%频率F2信号的幅度F1=50;%信号1频率(Hz)F2=75;%信号2频率(Hz)Fs=256;%采样频率(Hz)P1=-30;%信号1相位(度)P2=90;%信号相位(度)N=256;%采样点数t=[0:1/Fs:N/Fs];%采样时刻%信号;S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);%显示原始信号plot(S);title('原始信号');figure;Y=fft(S,N);%做FFT变换Ayy=(abs(Y));%取模plot(Ayy(1:N));%显示原始的FFT模值结果title('FFT模值');figure;Ayy=Ayy/(N/2);%换算成实际的幅度Ayy(1)=Ayy(1)/2;F=([1:N]-1)*Fs/N;%换算成实际的频率值plot(F(1:N/2),Ayy(1:N/2));%显示换算后的FFT模值结果title('幅度-频率曲线图');figure;Pyy=[1:N/2];fori=1:N/2Pyy(i)=phase(Y(i));%计算相位Pyy(i)=Pyy(i)*180/pi;%换算为角度end;plot(F(1:N/2),Pyy(1:N/2));%显示相位图title('相位-频率曲线图')图2滤波器特性5总结使用FIR数据滤波和FFT变换可以满足介质损耗测试仪的理论算法:（1）设计FIR滤波器时，当待处理的数据频率一定，而采样频率过高时，过渡带将严重偏离幅频特性坐标轴的中间位置，如图4所示，这样的后果是将大幅度提高FIR滤波器的阶数，最终后果是运算量将增大。FIR运算量计算是FIR_NUM*DATA_NUM次乘加运算。例如当滤波器为51阶、采样1024点数据，则需要51×1024＝52224次乘加运算；（2）根据需求设计的是50HZ陷波器，对50HZ的工频干扰衰减65DB，约衰减1778倍；对49.5HZ和50.5HZ（既是图9中标注的位置0.06875=49.5/[45*32/2]）的衰减为47.6DB,约240倍；对46.5HZ和53.5HZ,衰减3DB，约衰减1.4倍。对小于46HZ和大于54HZ的频率则完全可以通过，所设计的陷波器的阶数是801阶。(3)由于测试频率是由自己设计的变频电源产生，一般变频电源的精度可以达到0.02HZ的频偏，既是只需设计阻带在45.5HZ~54.5HZ的陷波器即可，这样可以大量降低陷波器的阶数。(4)经过测试使用整形查表法，进行FFT计算存在精度丢失等问题，如果为了保证精度需要在计算过程中不断判断数据的大小，以便适当的对数据进行缩放，则增加相应的运算量不可能大幅度减少FFT运算时间，所以打算不再使用FFT整形查表法进行FFT计算；(5)经过测试发现FFT浮点查表法和FFT浮点计算法在计算结果上差别不大，约在差别0.001%。且FFT浮点查表的时间远低于FFT浮点计算法，所以打算采用浮点查表法进行FFT计算；(6)经过测量发现：使用C8051F060单片机时，如使用其通过浮点查表法计算FFT，如在一个周期20ms