BSE.zip

function [We]=DBSE(X,R,taoind) % % Using Joint Generalized Eigenvectors of a Set of Covariance Matrix Pencils for Deflationary Blind Source Extraction % with fixed step size % X .......... observation signals % R .......... number of source signals % taoind .......... a vector of time lags�� e.g., [1 2 3 4 5 6 7 8 9 10] % We............... estimated extraction matrix [J,Lsource]=size(X); for i = 1:J X(i,:) = X(i,:) - mean(X(i,:)); end Rx0 = (1/Lsource)*X(:,1:Lsource)*X(:,1:Lsource)'; Rx0 = (1/2)*(Rx0+Rx0'); [U,D] = eig(Rx0); [puiss,sd] = sort(abs(diag(D))) ; rangeW = J-R+1:J ; % indices to the R most significant directions scales = sqrt(puiss(rangeW)) ; % scales if R < J scales = sqrt( puiss(rangeW) - mean(puiss(1:J-R)) ); end WMatrix = diag(1./scales) * U(1:J,sd(rangeW))' ; % whitener XW = WMatrix*X; %% N=length(taoind); MatrixNum = N; Ht = []; for i = 1:N sti = taoind(i)+1; slen = Lsource - sti + 1; Rx = (1/slen)*XW(:,sti:Lsource)*XW(:,1:Lsource-sti+1)'; Rx = (1/2)*(Rx + Rx'); Ht = [Ht Rx]; end %% Hvec = []; cornum = 1; for i = 1:N-1 for j = i+1:N mrange = 1:R; mrange= mrange + (i-1)*R; H1 = Ht(:,mrange); mrange = 1:R; mrange= mrange + (j-1)*R; H2 = Ht(:,mrange); for s = 1:R-1 for t = s+1:R H12 = H1(s,:).'*H2(t,:) - H1(t,:).'*H2(s,:); H21 = H2(t,:).'*H1(s,:) - H2(s,:).'*H1(t,:); Hvec = [Hvec vecsymm(H12+H21)]; end end end end %% [a,NumH] = size(Hvec); C = Hvec*Hvec'; %--------------------------------------- iteNum = 1000; tao = 0.005; %the step size should be chosen according to the variance of observation signals threshold = 0.0005; epsilon = 1.0e-25; %the threshold should be chosen according to the variance of observation signals Ae=[]; %First extraction vector x = ones(R,1); x = x/sqrt(x'*x); X = x*x'; error=1; for ite = 1:iteNum vecX = vecsymmx(X); vecU = vecX; g = 2*C*vecX; vecX = vecX - tao*g; X = unvecsymmx(vecX,R); Y = X; Y = Y/sqrt(trace(Y*Y')); [X,u,S] = shrinkage(Y,threshold); if S(2,2) == 0 [X,u] = shrinkage_dominant(Y); %error=1-abs(x'*u)/(sqrt(x'*x)*sqrt(u'*u)); error=1-abs(vecX'*vecU)/(sqrt(vecX'*vecX)*sqrt(vecU'*vecU)); end x = u; if(error<=epsilon) break; end end Ae=[Ae x]; %%Extraction vectors 2...R-1 for num=1:R-2 [Vx,Dx] = eig(Ae*Ae'); [Yx,Ix] = sort(abs(diag(Dx))); Vx = Vx(:,Ix); a = ones(R-num,1); x = Vx(:,1:R-num)*a; x = x/sqrt(x'*x); X = x*x'; vecx1 = vecsymm(Ae*Ae'); C = C + vecx1*vecx1'; error=1; for ite = 1:iteNum vecX = vecsymmx(X); vecU = vecX; g = 2*C*vecX; vecX = vecX - tao*g; X = unvecsymmx(vecX,R); Y = X; Y = Y/sqrt(trace(Y*Y')); [X,u,S] = shrinkage(Y,threshold); if S(2,2) == 0 X = shrinkage_dominant(Y); %error=1-abs(x'*u)/(sqrt(x'*x)*sqrt(u'*u)); error=1-abs(vecX'*vecU)/(sqrt(vecX'*vecX)*sqrt(vecU'*vecU)); end x = u; if(error<=epsilon) break; end end Ae = [Ae x]; end %%The last extraction vector [Vx,Dx] = eig(Ae*Ae'); [Yx,Ix] = sort(abs(diag(Dx))); Vx = Vx(:,Ix); x = Vx(:,1:1); x = x/sqrt(x'*x); X = x*x'; vecx1 = vecsymm(Ae*Ae'); C = C + vecx1*vecx1'; error=1; for ite = 1:iteNum vecX = vecsymmx(X); vecU = vecX; g = 2*C*vecX; vecX = vecX - tao*g; X = unvecsymmx(vecX,R); Y = X; %Y = Y/sqrt(trace(Y*Y')); %[X,u,S] = shrinkage(Y,threshold); [X,u] = shrinkage_dominant(Y); %error=1-abs(x'*u)/(sqrt(x'*x)*sqrt(u'*u)); error=1-abs(vecX'*vecU)/(sqrt(vecX'*vecX)*sqrt(vecU'*vecU)); x = u; if(error<=epsilon) break; end end Ae = [ Ae x ]; We = Ae'*WMatrix; return