自适应滤波器的实现

自适应滤波器的实现0 简介所谓自适应滤波器，就是利用前一时刻己经获得的滤波器参数的结果，自动地调节现在时刻的滤波器参数，以适应信号和噪声未知的或随时间变化的统计特性，从而实现最优滤波。该系统就是利用前一时刻已经获得的滤波器参数等结果，自动的调节现在有的时刻的滤波器参数，以适应信号或噪声或未知的或随时间变化的统计特性，从而实现最优滤波，如图1所示，自适应滤波器的相关输入信号、输出信号、期望响应和误差信号构成，一旦输入信号发生变化，自适应滤波器能够跟踪其变化并自动调节参数，使滤波器性能重新达到最佳。因此自适应滤波器具有“自我调节”和“跟踪”能力12。未知系统自适应IIR滤波器PSO算法图1 自适应滤波器结构图根据环境的改变，使用自适应算法来改变滤波器的参数和结构，这是自适应滤波器自身最明显也是其最大的优势。一般情况下，不改变自适应滤波器的结构，而自适应滤波器的系数是由自适应算法更新的时变系数，即其系数自动连续地适应于给定信号，以获得期望响应。自适应滤波器的最重要的特征就在于它能够在未知环境中有效工作，并能够跟踪输入信号的时变特征。自适应滤波器的输入输出关可用如下的差分方程来表示： （8）公式中，分别表示滤波器的输出和输入，并且M（）是滤波器的阶数。在系统识别领域中，上述差分方程对应的传递函数为： （9）可以将自适应IIR滤波器的设计当作优化问题，将均方误差（MSE）作为目标函数，即 （10）公式中：N是该滤波器的评价长度，是对象和滤波器的输出13。1 系数编码在自适应滤波器的应用中，系统识别是一个研究的热点，即：利用某种算法来确定滤波器的系数，（本文中使用PSO算法作为应用算法），使得（8）式中的均方误差函数 （15）为最小。显然这是一个多维变量的寻优问题，因此使用PSO算法来求解上述优化问题。在应用PSO算法设计滤波器的时候，适应度函数的求解是一个很重要的问题，该算法通过适应度函数来确定粒子当前位置的优劣，因此粒子的适应度函数值越小表明该粒子的确定的参数越优，算法中的粒子就对应较优的滤波器系数。算法结束之后，整个运行期间适应度之最小的粒子所代表的参数值就是算法得到的最优解，也就是最终得到的滤波器系数。这里选择的函数式（8）作为IIR滤波器设计的适应度函数。在应用PSO算法设计自适应滤波器的时候使用粒子所代表的参数值作为目标滤波器的系数。这样就存在了这样的一个问题，就是如何分配这些被确定的参数作为系数，也就是一个将滤波器系数编码为粒子参数的问题。为了较好的运用PSO算法通过适应度函数来确定目标滤波器系数，将这些系数进行了编码，以形成PSO算法中粒子的系数。根据PSO算法的本身的特点，可以直接使用实数来表示各个参数因此可以直接将滤波器系数作为参数应用到自适应IIR滤波器的设计中来。从而粒子可以表示成如下所示： （16）采用编码结构的粒子是一个维的向量。最后会得到输出的gbest，也就是全局最优的值，这样就会得到期望的自适应的滤波器的系数6,7。2 PSO算法实现模块模块（）：设定粒子群最优化算法的参数，包括群体大小，参数维数，加权因子等；虽然粒子群优化算法并没有许多需要调节的参数，但是这些参数直接影响着算法的性能以及收敛问题。下面分别介绍这些参数及其设置。加权参数：,一般设定在0.9到0.4之间。通过大量实验证明：如果算法的进行而线性减小那么将显著改善算法的收敛性能。为最小权值，为最大权值。加速常数：c1和c2代表了粒子向自身极值pbest和全局极值gbest推进的随机加速权值。小的加速常数值，可使粒子在远离目标区域内振荡，而大的加速常数值可使粒子迅速向目标区域移动，甚至又离开目标区域。如果c1=c2=0，则粒子将以当前速度飞行，直到边界。此时，由于粒子只能搜索有限的区域，故很难找到好的解。当c1=0，则粒子没有认知能力，亦即“只有社会（social-only）”的模型。在粒子相互作用下，算法有能力达到新的搜索空间，其收敛速度比标准算法更快，但碰到复杂问题，比标准算法更容易陷入局部极点。当c2=0，则粒子之间没有社会信息共享，亦即“只有认知（cognition）”的模型。由于个体之间没有交互，一个规模为m的群体等价于m个单个粒子的运行，因而得到最优解的概率非常小。在实际应用中根据一些学者的研究设定c1=c2为2或者1.4。后来Clerc和其他一些研究者提出了一些其它值，但实际上这些也只是局限于某些问题，无法推广应用。粒子数目：即种群的大小一般取20至200。实验表明，对于大多数问题来说，20个粒子就可以取得很好的结果，不过对于比较难的问题或者特殊类别的问题，粒子数目可以取到100或200。另外，粒子数目越多，算法搜索的空间范围就越大，也就更容易发现全局最优解。当然，算法运行的时间也较长；粒子长度：粒子长度就是问题的维数，它由具体优化问题确定。比如，对于一个具有10个参数的寻优问题，粒子长度就可以设置为10；粒子范围：粒子范围同样由具体优化问题确定，通常把问题的参数取值范围设置为粒子的范围。另外，粒子每一维可以设置不同的范围；粒子最大速率：粒子最大速率决定粒子在一次飞行中可以移动的最大距离，因此必须限制粒子最大速率，否则粒子就可能跑出搜索空间。粒子最大速率通常设定为粒子范围的宽度。比如，粒子范围为（-10，10），则粒子最大速率就可以设置为2；算法终止条件：与遗传算法相似，粒子群优化算法的终止条件一般可以设置为达到最大进化世代数或者满足一定的误差准则；适应度函数：粒子群优化算法的适应度函数选择比较简单，通常可以直接把目标函数作为适应度函数。当然，也可以对目标函数进行变换，变换方法可以借鉴遗传算法中的适应度函数变换方法，比如线性变换、幂函数变换等16。为了充分的验证PSO算法在IIR滤波器设计中的性能，对以下的例子进行了随机试验仿真。输入信号采用随机生成的数组，暂时不考虑外界的干扰影响，参数的搜索范围是-2，+2。PSO算法的初始参数设置为：种群的个数为80，最大进化代数也是400，其他的设置按照PSO算法经典设置进行配置17,18。算例（1）：未知系统和自适应IIR滤波器的传递函数分别为6：因此通过PSO算法进行最优化求解然后最后获得的具体最优解。下面就是进行试验的具体流程和步骤：该步骤是应用PSO算法来设计自适应滤波器的一个基本步骤，只有将PSO算法中各个参数确定下来才能保证PSO算法在应用过程中顺利的进行而不会因为自身参数问题使粒子寻优过程中出现飞离限定范围等问题（具体应用到本次试验的参数设计和函数确定在程序中有清楚的表达出）。相应的MATLAB程序如下：初始化PSO算法的各个参数swarmSize=80;%群体大小dimSize=6;%粒子维数upbnd=2;%参数的取值范围上限lowbnd=-2;%参数的取值范围下限wmax=0.9;%最大加权因子wmin=0.4;%最小加权因子vMax=2;%最大速率%vMin=-2;%这个参数是没有意义的itmax=400; %Maximum iteration numberDA=1e-90; %程序跳出的控制条件N=200; %评价长度在本文中，结合实际情况我将参数设定的比较简单，然而，在不同的应用环境下可以随机改变这些参数的具体状况。设定参数为最大加权因子wmax=0.9，最小加权因子wmin=0.4，群体的大小设定为dimsize=6，swarmsize=80；这样在运算仿真的时候等待的时间会比较短一点。c1=c2=2，最大的迭代次数为itmax=400，为了方便将其他的参数取值范围设定为-2，2。加权系数变化按照下面的公式进行：其中iter表示当前的迭代次数。这样就将PSO算法部分的参数进行了初步的定义，之后就可以进行算法的运算来确定滤波器的系数。模块（2）：对粒子的速度和位置进行初始化并计算粒子的适应度值。由图1、图2可知，在设计自适应滤波器的时候必须给粒子一个初始的参量，即：粒子的初始化速度和初始的位置，然后以这个初始的速度和位置进行依次迭代寻优。相应MATLAB程序如下：xin=rand(1,N);%输入粒子的初始值pos = rand(swarmSize, dimSize) ;for i=1:swarmSizemeanpos=mean(pos(i,:);pos(i,:)=pos(i,:)-meanpos;%对粒子进行初始化endvel=(upbnd-lowbnd)*rand(swarmSize,dimSize)+ones(swarmSize,dimSize)*lowbnd;%对速度进行初始化通过这段程序可以给粒子的速度和位置分别定义一组初值，因为通过指令：rand，可以将粒子的初始速度和位置进行预定义，就是给他们的速度和位置一组随机数，然后就能够计算他的最初适应度，检测它是否符合PSO算法的约束条件，这样才能够将算法引入下一个循环，即：for swarmIndex=1:swarmSizebx=pos(swarmIndex,1:3);ax=1 pos(swarmIndex,4:6);dk=filter(bx,ax,xin);fpBest(swarmIndex,1)=sum(dn-yn).2)/N;%计算最初的粒子适应度值end当得到的粒子满足上面的f函数，即满足适应度函数时，就可以进行下一步的运算，计算粒子当前最优值pbest，和全局当前最优值gbest，就是进入下一个运算模块。模块（3）：根据式子（3）和（4）来更新当前粒子的速度和位置，并重新计算粒子的适应度值。相应的程序如下所示：更新粒子的速度：for dimIndex=1:dimSizeif V(dimIndex) vMaxvel( swarmIndex, dimIndex) = vMax;elseif V vMaxpos(swarmIndex,dimIndex)=vMax;elseif P(dimIndex)-vMaxpos(swarmIndex,dimIndex)=-vMax;elsepos(swarmIndex,dimIndex)=P(dimIndex);end模块（4）判断是否更新粒子的个体极值pbest，和粒子群全局的极值gbest。相应的程序如下：for swarmIndex = 1: swarmSizebx=pos(swarmIndex,1:3);ax=1 pos(swarmIndex,4:6);dk=filter(bx,ax,xin);B(1, swarmIndex)=sum(dk-yk).2)/N;if B(1,swarmIndex)fpBest(swarmIndex)fpBest(sw