past2.m

基于投影逼近的子空间跟踪（PAST）方法

clear all;clc;
n=0:0.0001:0.3999;                   %采样频率为10kHz，4000个点
s1=sign(cos(2*pi*155*n));
s2=sin(2*pi*800*n);
s3=sin(2*pi*90*n);
s4=sin(2*pi*300*n+6*cos(2*pi*60*n));
s5=-2*rand(1,4000)+1;               %仿真5个源信号
S=[s1;s2;s3;s4;s5];
beta=0.99;                          %参数beta接近于1
lamda=0.001;                        %参数lamda比较小的步长
N=200;                                %独立进行200次实验

for t=1:N
    A=randn(5,5);                   %随机的混合矩阵
    Wt=eye(5,5);                    
    W=eye(5);
    P=eye(5);                       %几个矩阵自适应迭代初值设为单位阵
    for i=1:4000                    %仿真的4000个点
        x=A*S(:,i);                 %5个信号混合
        y=Wt*x;                     
        Wt=Wt-lamda*(y*y'-eye(5,5))*Wt;%预白化即白化矩阵更新
        z=tanh(W'*y);               %PAST自适应的几个更新方程
        h=P*z;
        m=h/(beta+z'*h);
        Ptri=P-m*h';
        for r=1:5
            for s=r:5
                Ptri(s,r)=Ptri(r,s);%将上三角的值复制到下三角成为对称矩阵
            end
        end
        P=Ptri/beta;
        e=y-W*z;
        W=W+e*m';                   %分离矩阵更新
        Se(:,i)=W'*y;               %利用分离矩阵对源信号进行估计
        
        C=W'*Wt*A;                  %误差计算
        Esum=0;
        for p=1:5
            maxr=max(abs(C(p,:)));
            maxc=max(abs(C(:,p)));
            for q=1:5
                Esum=Esum+abs(C(p,q))/maxr+abs(C(p,q))/maxc;
            end
        end
        E(t,i)=sqrt(Esum-10);
    end
end
Em=sum(E)/N;                    %计算200次实验的均值
figure;
plot(Em,'b');hold on;                      %输出误差曲线
% xlabel('采样点数');
% ylabel('串音误差均值');
% legend('g(y)=tanh(y)','g(y)=y','g(y)=y^3');

st=3001;ed=3200
figure;
subplot(5,1,1)
plot(Se(1,st:ed));hold on;plot(S(2,st:ed),'r--');
ylabel('Se1')
title('原始信号(红)与分离后信号对比图(蓝)')

subplot(5,1,2)
plot(Se(2,st:ed));hold on;plot(S(3,st:ed),'r--');
ylabel('Se2')

subplot(5,1,3)
plot(Se(3,st:ed));hold on;plot(S(1,st:ed),'r--');
ylabel('Se3')

subplot(5,1,4)
plot(Se(4,st:ed));hold on;plot(S(5,st:ed),'r--');
ylabel('Se4')

subplot(5,1,5)
plot(Se(5,st:ed));hold on;plot(S(4,st:ed),'r--');
ylabel('Se5')