www.pudn.com > ln.rar > b_20_noisy_free.m

clc;clear
load ABio7.mat;
%%LP的阶数P=20；
P=2;
N=length(benchcicbar);

b=[0.35 0.89];  %%neng ti qu chu ,buguo you zaosheng.
%b=randn(1,P);
for i=2:5
    figure(1);
    s(i-1,:)=benchcicbar(i,:);
    subplot(4,1,i-1);plot(s(i-1,:)); 
    %%calculate the NMSPE:
    for j=1:P
        sum=0;
        %temp=s(i-1,[zeros(1,j) length(benchcicbar)-j]);
        dz(i-1,:)=[zeros(1,j) s(i-1,1:(N-j))];
        sum=sum+b(j)*dz(i-1,:);
    end    
    e(i-1,:)=s(i-1,:)-sum;
    
    eMSPE(i-1)=mean(e(i-1,:)*e(i-1,:)');
    eNMSPE(i-1)=mean(e(i-1,:)*e(i-1,:)')/mean(s(i-1,:)*s(i-1,:)');
end;    

%%mixing matrix:A
A=[-0.9370 0.9448 0.3651 -0.2617;-0.8263 -0.3478 -0.2872 0.6279;
    0.3730 0.1975 0.9953 0.8325;0.8955 0.3578 0.5649 0.6257];

X=A*s;    %%% 观测信号X;四个
for i=1:4
    figure(2);
    subplot(4,1,i);plot(X(i,:));title('观测信号X(i)');  %%.....
end;    

%%%进行BSE。
for i=1:4
    for j=1:2
        sum=0;
        dX(i,:)=[zeros(1,j) X(i,1:(N-j))];
        sum_X=sum+b(j)*dX(i,:);
    end    
    e_X(i,:)=X(i,:)-sum_X;
end; 
%%%首先是没有进行白化的BSE
    w=ones(1,4)';
    w=w/norm(w);
    y(1)=w'*X(:,1); %%........
    y(2)=w'*X(:,2); %%........
    y(3)=w'*X(:,3); %%........
    
for k=3:N
    y(k)=w'*X(:,k); 
    miu=0.00025;    
    e(k)=y(k)-b(1)*y(k-1) - b(2)*y(k-2);
    
    belta_e=0.98;
    belta_y=0.98;
    sigma_e2=0.02;
    sigma_y2=0.02;
    sigma_e2=belta_e*sigma_e2+(1-belta_e)*(e(k)^2);
    sigma_y2=belta_y*sigma_y2+(1-belta_y)*(y(k)^2);
    
    w=w-miu*(e(k)*e_X(:,k)-sigma_e2*y(k)*X(:,k)/sigma_y2)/sigma_y2;
    w=w/norm(w);
end; 

y=w'*X;
w'*A
figure(4);
subplot(411),plot(y);title('the extracted signal wihtout whiten');
subplot(412),plot(s(4,:));title('the source signal');


%%%%%进行白化后的BSE
%%白化
for i=1:4
    X(i,:)=X(i,:)-mean(X(i,:));
end;  
%[U0,S0,V0]=svd(X);    
%[U1,S1,V1]=svd(X*X');
%X=S1^(-1/2)*U0'*X;

[E,D]=eig(X*X');
X=D^(-1/2)*E'*X;   %%与上一样。

for i=1:4
    for j=1:2
        sum=0;
        dX(i,:)=[zeros(1,j) X(i,1:(N-j))];
        sum_X(i,:)=sum+b(j)*dX(i,:);
    end    
    e_X(i,:)=X(i,:)-sum_X(i,:);
end; 

    w=ones(1,4)';
    w=w/norm(w);
    y(1)=w'*X(:,1); %%........
    y(2)=w'*X(:,2); %%........
    y(3)=w'*X(:,3); %%........
    
for k=3:N
    y(k)=w'*X(:,k);     
    e(k)=y(k)-b(1)*y(k-1)-b(2)*y(k-2);
    miu2=0.4;
    w=w-miu2*e(k)*e_X(:,k);
    w=w/norm(w);
    %w=w/sqrt(w'*(eye(4)-sigma_v2*D^(-1/2)*E'*E*(D^(-1/2))')*w);
end;    
%y=w'*X;

figure(5);
subplot(411),plot(y);title('the extracted signal with whiten');
subplot(412),plot(s(4,:));title('the source signal_4');
subplot(413),plot(s(3,:));title('the source signal_3');