递归最小二乘法（RLS）自适应滤波算法摘要所谓自适应，从通俗意义上讲，就是这种滤波器能够根据输入信号统计特性的变化自动调整其结构参数，以满足某种最佳准则的要求.自适应滤波器所采用的最佳准则由最小均方误差准则、最小二乘准则、最大信噪比准则和统计监测准则等。
自适应滤波理论和技术是统计信号处理和非平稳随机信号处理的主要内容,它可以在无需先验知识的条件下,通过自学习适应或跟踪外部环境的非平稳随机变化,并最终逼近维纳滤波和卡尔曼滤波的最佳滤波性能。
因而，自适应滤波不但可以用来检测确定性信号，而且可以检测平稳的或非平稳的随机信号。
自适应技术应用包括自适应谱线检测增强与谱估计方法、自适应噪声抵消技术、自适应均衡技术、只适用阵列处理与波束形成以及自适应神经网络信号处理等内容。
关键词：递归最小二乘法；自适应滤波；滤波器设计；自适应算法；1 引言滤波可分为经典滤波和现代滤波。
经典滤波要求已知信号和噪声的统计特性，如维纳滤波和卡尔曼滤波。
现代滤波则不要求己知信号和噪声的统计特性，如自适应滤波。
自适应滤波的原理就是利用前一时刻己获得的滤波参数等结果,自动地调节当前时刻的滤波参数，从而达到最优化滤波。
自适应滤波一般包括3个模块：滤波结构、性能判据和自适应算法。
其中,自适应滤波算法一直是人们的研究热点,包括线性自适应算法和非线性自适应算法，非线性自适应算法具有更强的信号处理能力，但计算比较复杂,实际应用最多的仍然是线性自适应滤波算法，线性自适应滤波算法的种类很多，有LMS自适应滤波算法、R路自适应滤波算法、变换域自适应滤波算法、仿射投影算法、共扼梯度算法等.自适应滤波器主要包括滤波器的结构和自适应算法两部分，这两部分不同的变化与组合,可以导出许多不同形式的自适应滤波器.所谓的自适应滤波，就是利用前一时刻以获得的滤波器参数的结果，自动地调节现时刻的滤波器参数，以适应信号和噪声未知的或随时间变化的统计特性，从而实现最优滤波。
自适应滤波器实质上就是一种能调节其自身传输特性以达到最优的维纳滤波器。
4.信号的函数表达式为：()()()()sin(2100) 1.5sin(2300)sin(2200)x t t t A t t dn t n t πππ=++++，其中，()A t 为一随时间变化的随机过程，()dn t 为经过390—410Hz 带通滤波器后的高斯白噪声，()n t 为高斯白噪声，采样频率为1kHz ，采样时间为2.048s 。
（1）利用现代信号处理的知识进行信号谱估计；（2）利用现代信号处理知识进行信号的频率提取； （3）分别利用Winner 滤波和Kalman 滤波进行去噪； （4）利用Wigner-Ville 分布分析信号的时频特性。
（1）：利用现代信号处理的知识进行信号谱估计：经典谱估计中两种主要的方法为直接法和间接法，其中间接法则先根据N 个样本数据()x n 的样本自相关函数()()()1*01,01N x n R k x n k x n k M N-==+=⋅⋅⋅∑，，，（4.1）其中1M N ≤<，且()()*x x R k R k -=。
计算样本自相关函数的Fourier 变换，得到功率谱()()Mjk x x k MP R k e ωω-=-=∑（4.2）周期图方法估计的功率谱为有偏估计，可通过加窗来减少其偏差。
定义为 ()()()2101N jn x n P x n c n e NWωω--==∑ （4.3）式中()()122112N n W c n C d NNππωωπ--===∑⎰（4.4）式中，()C ω是窗函数()c n 的Fourier 变换。
功率谱估计程序为： clear clcclose all hidden sf=1000;nfft=2048; t=0:1/1000:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410;wc1=2*f1/sf; wc2=2*f2/sf; %归一化频率f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; a(1,:)=y;a(2,:)=y.*sin(y); x=a(1,:);y=a(2,:)-a(1,:);f=0:sf/nfft:sf/2-sf/nfft;w=boxcar(nfft);%加矩形窗 z=psd(y,nfft,sf,w,nfft/2); nn=1:nfft/2;plot(f(nn),abs(z(nn))); xlabel('频率(Hz)'); ylabel('幅值'); grid on;图4.1 功率谱估计结果图(2).信号频率的提取用离散傅立叶算法离散傅立叶算法程序 clear clcclose all hidden sf=1000;nfft=2048; t=0:1/1000:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410; wc1=2*f1/sf; wc2=2*f2/sf;050100150200250300350400450500200400600800频率(Hz)幅值%归一化频率f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通和带阻weigh=[1 1 1 ];%设置带通和带阻权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; t2=(0:nfft-1)/sf;f=(0:nfft-1)*sf/nfft;y1=abs(fft(y));f=f(1:nfft/2);y1=y1(1:nfft/2);plot(t,y);title('原始信号');axis([0 2.047 -6 8]);plot(f,y1);title('fft频率提取');axis([0 500 0 1000]);xlabel('f/Hz');grid on;原信号时间（t）图4.2 原始信号时域图图4.3 信号频谱(3)分别利用Winner 滤波和Kalman 滤波进行去噪；clear all close allM=100;%维纳滤波器阶数 sf=1000;nfft=2048; L=nfft;t=0:1/1000:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410; wc1=2*f1/sf; wc2=2*f2/sf; %归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; phixx=xcorr(y,y); for i=1:M for j=1:MRxx(i,j)=phixx(i-j+L); end ends=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); phixs=xcorr(y,s); for i=1:Mrxs(i)=phixs(i+L); endh1=(inv(Rxx))*rxs';2004006008001000fft 频率提取f/Hz%获得理想FIR滤波器系数h1AA=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); for i=1:Mh(i)=AA(i);end%绘图比较估计滤波器与实际滤波器figurek=1:M;plot(k,h(k),'r',k,h1(k),'b');title('Ideal h(n) & Calculated h(n)');legend('Ideal h(n)',' Calculated h(n)');xlabel('n');ylabel('h(n)');%比较理想输出与实际输出v=D+N;S=conv(h,v);SI(1)=S(1);LL1=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); for i=2:LSI(i)=LL1(i);endfigurek=1:L;plot(k,s(k),'r',k,SI(k),'b');title('s(n) VS. SI(n)');legend('s(n)','SI(n)',0);xlabel('n');ylabel('Ideal Output');hold onSR=conv(h1,y);figurek=1:L;plot(k,s(k),'r',k,SR(k),'b');title('s(n)VS. SR(n)');legend('s(n)去噪前','SR(n)去噪后',0);xlabel('n');ylabel('Actual Output');图4.4 Winner 滤波去噪图Kalman 滤波程序 clear; clc;Fs=1000; nfft=2048;t1=0:1/Fs:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,2); f1=390;f2=410; wc1=2*f1/Fs; wc2=2*f2/Fs; wc2=2*f2/sf; %归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N; x1=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200); a1=-1.352;a2=1.338;a3=-0.662;a4=0.240;A=[-a1 -a2 -a3 -a4;1 0 0 0;0 1 0 0;0 0 1 0];%状态转移矩阵 H=[1 0 0 0];%观测矩阵Q=[1 0 0 0;0 0 0 0;0 0 0 0;0 0 0 0];%状态噪声方差 R=1;%观测噪声方差阵X(:,1)=[x(4);x(3);x(2);x(1)];p(:,:,1)=[10 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];%一步预测误差方针 %开始滤波 for k=2:nfftp1(:,:,k)=A*p(:,:,k-1)*A'+Q;%p1(:,:,k)即是一步预测误差的自相关矩阵，它是4*4的矩阵，取不同的k 值就构成了一个三维矩阵K(:,k)=p1(:,:,k)*H'/(H*p1(:,:,k)*H'+R); %K(:,:,k)是增益矩阵,对于固定的k 值它是4*1矩阵，取不同的k 值就是三维矩阵s(n)VS. SR(n)nA c t u a l O u t p u tX(:,k)=A*X(:,k-1)+K(:,k)*[x(k)-H*A*X(:,k-1)]; %X(:,k)是估计值,4*1矩阵p(:,:,k)=p1(:,:,k)-K(:,k)*H*p1(:,:,k);%p(:,:,k)是估计误差的自相关矩阵，4*4矩阵的三维矩阵end%结束一次滤波%绘图t=1:nfft;figure(2);plot(t,x1,'k-',t,x,'r-',t,X(1,:),'b-.');title('卡曼滤波去噪')legend('真实轨迹','观测样本','估计轨迹');grid on;卡曼滤波去噪n图5 Kalman滤波去噪图(4) 利用Wigner-Ville分布分析信号的时频特性MATLAB程序clear;clc;Fs=1000;nfft=2049;t1=0:1/Fs:2.048;A=normrnd(0,1,1,2049);N=wgn(1,2049,2);f1=390;f2=410;wc1=2*f1/Fs;wc2=2*f2/Fs;%归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通或带阻，1为带通，0为带阻weigh=[1 1 1 ];%设置通带和阻带的权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N; figure(8)tfrwv(x');xlabel('时间t');ylabel('频率f');0.50.450.40.350.30.250.20.150.10.05图6 幅频特性图。
