bss_gbjz_ep_eeg_end.m

语音信号的盲分离

function main()  %主程序

clear all;
%------------------------常量初值--------------------- 
load eeg.txt

 load c07_0.dat;
 epsig=c07_0;

 
 hello=epsig/5;



N=2000;%快拍数

 hello2=hello.';
 hello=zeros(1,N);
 hello(1,1:N)=hello2(1,1:N);

eeg2=eeg.';
eeg=zeros(1,N);
eeg(1,1:N)=eeg2(1,1:N);



M=5;%传感器个数
n=0:N-1;
P=2;%信号源个数

%-----------------------------------------------------

source_angle=-20;%期望信号角度
sita_s=source_angle*pi/180;%期望信号方位

source_angle_2=40;%期望信号角度
sita_s_2=source_angle_2*pi/180;%期望信号方位

 

%-----------------------------------------------------

 power_source=10^(20/20);%期望信号功率

power_noise=10^( 0/20);%干扰信号功率

%-----------------------------------------------------
alpha=1.4;%稳定分布特征指数，0<alpha<2
beta=0;
gama=1;
miu=0;

p=1.2;%分数阶矩的阶数p<alpha<2
%-----------------------------------------------------

 source=power_source*hello(1,1:N);
 source=source-mean(source);

  source_2= power_source*eeg(1,1:N);
 


 noise =power_noise*sasd(N,alpha).'; % 复稳定分布
 %noise =power_noise*casd(N,alpha,beta,gama,miu); % 复稳定分布



 ii=[0:M-1].';
 asita_s=exp(-i*ii*pi*(sin(sita_s)))/sqrt(M);  %  期望信号方向矢量
 asita_s_2=exp(-i*ii*pi*(sin(sita_s_2)))/sqrt(M);  %  期望信号方向矢量

 
 
    X= (asita_s*source  + asita_s_2*source_2)*sqrt(M) + ones(M,1)*noise%+  randn(M,N);  % 阵列接收信号
   % X= (asita_s*source  + asita_s_2*source_2)*sqrt(M)+  randn(M,N);  % 阵列接收信号
  
    
        % for nu=1:N 
        % X(:,nu)=X(:,nu)-mean(X(:,nu));
        %end
 %for nu=1:M 
  %   X(nu,:)=X(nu,:)-mean(X(nu,:));
  %end
     
    X=X./50;  
     
     
     
  %%%%%%---1----%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
     
     
   
      cov_X=covariationmatrix(X,p,alpha);%计算 X 的共变矩阵
    % cov_X=flocmatrix(X,p);%计算 X 的FLOC矩阵
  
 [vec,di ]=eig(cov_X)
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %  V=([di(1,1) 0;0 di(2,2)]^(-1/2))*[vec(:,1) vec(:,2)].';
   % Z=V*X ;
      V=[vec(:,1) vec(:,2)]*([di(1,1) 0;0 di(2,2)]^(-1/2));

      %V=[vec(:,1) vec(:,2) vec(:,3)]*([di(1,1) 0 0;0 di(2,2) 0;0 0 di(3,3)]^(-1/2));
    % V=([di(1,1) 0;0 di(2,2)]^(-1/2))*[vec(:,1) vec(:,2)]';
    %V=[vec(:,M-1) vec(:,M)]*([di(M-1,M-1) 0;0 di(M,M)]^(-1/2));
  
    %  Z=V'*(abs(X).^(p-1).*conj(X)) ;%数据阵----FLOC化
     Z=V'*X ;%
    
    
     B= rand(P,P)+i*rand(P,P) ;
    % B= rand(M,M)+i*rand(M,M) ;
    B0=B;
    
    
    u=0.01;
 
   N_iter=N;%迭代次数
   
   
   
   for k=1:N_iter 
       
        Y(:,k)=B'*Z(:,k);
        
  
        B=B+u*( Z(:,k)-B*g1(Y(:,k)) )*g1(Y(:,k)');
        %B=B+u*( Z(:,k)-B*g(Y(:,k),p)) * g(Y(:,k).',p);
        norm_B_1(k)= abs(log10(norm(B-B0,2)));
    end      
  
  
     
      Y1=B'*Z ;
  
  
      
      
      
      
      
      
      
      
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
   cov_X=covariancematrix(X); %计算 协方差矩阵
  
 [vec,di ]=eig(cov_X) 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %V=[vec(:,1) vec(:,2)]*([di(1,1) 0;0 di(2,2)]^(-1/2));
      V=[vec(:,M-1) vec(:,M)]*([di(M-1,M-1) 0;0 di(M,M)]^(-1/2));
    % V=[vec(:,M-2) vec(:,M-1) vec(:,M)]*([di(M-2,M-2) 0 0;0 di(M-1,M-1) 0;0 0 di(M,M)]^(-1/2));
  
   Z=V'*X ;
    
  
     B= rand(P,P)+i*rand(P,P) ;
    B0=B;
    
    
    u=0.01;
 
   for k=1:N_iter 
       
        Y(:,k)=B'*Z(:,k);
        
  
        B=B+u*( Z(:,k)-B*g2(Y(:,k) ))*g2(Y(:,k)');
        %B=B+u*( Z(:,k)-B*g(Y(:,k),p ))*g(Y(:,k).',p);
        norm_B_2(k)= abs(log10(norm(B-B0,2)));
    end      
  
  
     
      Y2=B'*Z ;
  
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 figure(1);
 MMM=300;
 subplot(3,2,1);
plot(0.001*real(source(1:MMM)-mean(source(1:MMM))),'r');
title('a.source signal 1');
 %subplot(4,2,3);
%plot(real(X(1,1:MMM)),'r');
%title('c.混合信号1');
  subplot(3,2,3);
  plot(real(Y2(1,:)+0.1*rasd(N,alpha,beta,gama,miu)),'r');
   plot(real(Y2(1,1:MMM)),'r');
title('c.separate signal 1(SOS）');
 subplot(3,2,5);
  plot(0.1*real(Y1(1,1:MMM)),'r');
title('e.separate signal 1(FLOS）');

 subplot(3,2,2);
plot(0.01*real(source_2(1:MMM)-mean(source_2)),'b');
title('b.source signal 2');
 %subplot(4,2,4);
%plot(real(X(2,1:MMM)),'b');
%title('d.混合信号2 ');
 subplot(3,2,4);
 %plot(real(Y2(2,:)+0.1*rasd(N,alpha,beta,gama,miu)),'b');
 plot(real(Y2(2,1:MMM)),'b');
title('d.separate signal 2(SOS）'); 
 subplot(3,2,6);
 plot(0.1*real(Y1(2,1:MMM)),'b');
title('f.separate signal 2(FLOS）'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

 
     figure(2);
     MM=500;
     plot(1:MM,abs(norm_B_1(1:MM)-mean(norm_B_1)),'r-',1:MM,abs(norm_B_2(1:MM)-mean(norm_B_2)),'b-.');
     xlabel('iteration nunber');
     ylabel('MSE');
    %legend('SOS','FLOS')
       legend('FLOS','SOS')






 %-------------------------
 
   
     function y=g1(x)
             p=1.2;
             y=abs(x).^(p-2).*conj(x) ;

             
            
      
      
 function y=g2(x )
     y= tanh(x) ;


 
     
  

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %------------% 计算 协方差矩阵--------------------
 function r_X=covariancematrix(X)
 [M,N]=size(X);
 for in=1:M
        
         for jn=1:M
                         
             r_X_a(in,jn)=0;
             
             for nn=1:N
             r_X_a(in,jn)=r_X_a(in,jn)+X(in,nn)*conj(X(jn,nn));
             end
             
             r_X(in,jn)=r_X_a(in,jn)/N;
         end
     end  
     
     
 
     
     
     
 %------------% 计算 FLOC矩阵--------------------
 function r_X=flocmatrix(X,p)
 [M,N]=size(X);
 for in=1:M
        
         for jn=1:M
                         
             r_X_a(in,jn)=0;
             
             for nn=1:N
            % r_X_a(in,jn)=r_X_a(in,jn)+(X(in,nn))*abs(X(jn,nn))^(p-1)*conj(X(jn,nn));
              r_X_a(in,jn)=r_X_a(in,jn)+abs(X(in,nn))^(p-1)*conj(X(in,nn))*abs(X(jn,nn))^(p-1)*conj(X(jn,nn));
             end
             
             r_X(in,jn)=r_X_a(in,jn)/N;
         end
     end  


     
     
 %------计算 X 的共变矩阵------------
  function cov_X=covariationmatrix(X,p,alpha)
  
   [M,N]=size(X);
    for in=1:M
        
         for jn=1:M
                         
             cov_coef_a(in,jn)=0;
             cov_coef_b(in,jn)=0;
             
             for nn=1:N
             cov_coef_a(in,jn)=cov_coef_a(in,jn)+X(in,nn)*(abs(X(jn,nn)).^(p-2))*conj(X(jn,nn));
             cov_coef_b(in,jn)=cov_coef_b(in,jn)+(abs(X(jn,nn)).^(p-1));
             end
             
             c_p_alpha=2^(p+1)*gamma((p+2)/2)*gamma(-p/alpha)/(alpha*gamma(1/2)*gamma(-p/2));
             cov_X(in,jn)=(cov_coef_a(in,jn)/cov_coef_b(in,jn))*(cov_coef_b(in,jn)/(N*c_p_alpha)).^(alpha/p);
                       
         end
     end
     
 
     
     
%----------------产生 复全向 稳定分布信号-------------------
function X=casd(N,alpha,beta,gama,miu)

G1=normrnd(0,1,1,N);
G2=normrnd(0,1,1,N);
AA=rasd(N,alpha/2,1,(cos(pi*alpha/4))^2,miu);
X=AA.^(1/2).*(G1+i*G2);


function Y=rasd(N,alpha,beta,gama,miu)


% ----- 均匀分布 -----

%U=unifrnd(-pi/2,pi/2,1,N);

nuni=rand(1,N);
U=(nuni*pi)-pi/2;


% ----- 指数分布 -----

%W=exprnd(1,1,N);

zhi=rand(1,N);
%W=-log(1-zhi);
W=-log(zhi);




if alpha~=1
   X1=S(alpha,beta);
   X2=sin(alpha*(U+B(alpha,beta)))./(cos(U)).^(1/alpha);
   %X2=sin(alpha*(U-B(alpha,beta)))./(cos(U)).^(1/alpha);
   X3=(cos(U-alpha*(U+B(alpha,beta)))./W).^(1/alpha-1);
   %X3=(cos(U-alpha*(U-B(alpha,beta)))./W).^(1/alpha-1);
   X=X1.*X2.*X3;
else
    X=(2/pi)*( (pi/2+beta*U).*tan(U)-beta*log( 0.5*pi*W.*cos(U)./(pi/2+beta*U) ) );
    %X=(2/pi)*( (pi/2+beta*U).*(U)-beta*log( W.*cos(U)./(pi/2+beta*U) ) );
end



if alpha~=1
    Y=gama*X+miu;

else
    Y=gama*X+miu+(2/pi)*beta*gama*log(gama);
end




function y=B(alpha,beta)
y=atan(beta*tan(pi*alpha/2))/alpha;



function y=S(alpha,beta)
y=(1 + beta^2 * (tan(pi*alpha/2)) ^ 2  ) ^   (1/(2*alpha));






%-----------------------------------------------------------------------
%  This program generate an ARPHA-Stable random 
%  signal.

function [sas]=sasd(N,arpha)



pi2=pi/2;

beta=0;

zl=N;
% ----- 均匀分布 -----

nuni=rand(zl,1);
nunip=(nuni*pi)-pi2;

% ----- 指数分布 -----

zhi=rand(zl,1);
zhi1=-log(1-zhi);

% ----- ARPHA 稳定信号 ------

if arpha==1
    beta_a=beta;
    gamma_a=pi2;
    phi_0=0;
else
    ka=1-abs(1-arpha);
    beta_a=2*atan(beta*tan(pi*arpha/2))/(pi*ka);
    gamma_b=cos(pi*beta_a*ka/2);
    phi_0=-0.5*pi*beta_a*(ka/arpha);
end

ee=1-arpha;
tau=-ee*tan(arpha*phi_0);
a=tan(0.5*nunip);
b=tan(0.5*ee*nunip);

z=(cos(ee*nunip)-tan(arpha*phi_0)*sin(ee*nunip))./(zhi1.*cos(nunip));

if arpha==1
    d=0;
    B=0;
else
    d=(z.^(ee/arpha)-1)/ee;
    B=tan(b)./(0.5*ee*nunip);
end

y1=2*(a-b).*(1+a.*b)-tau*nunip.*B.*(b.*(1-a.*a)-2*a);
y2=(1-a.*a).*(1+b.*b);
sas=(y1./y2).*(1+ee*d)+tau.*d;




alphastring=num2str(arpha);