JADE盲分离算法仿真

JADE算法原理

JADE 算法首先通过去均值预白化等预处理过程得到解相关的混合信号，预处理后的信号构建的协方差矩阵变为单位阵，为后续的联合对角化奠定基础；其次，通过建立四阶累积量矩阵，利用高阶累积量的统计独立性等性质从白化后的传感器混合(观测)信号中得到待分解的特征矩阵；最后，通过特征矩阵联合对角化和Givens 旋转得到酉矩阵$U$，从而获得盲源分离算法中混合矩阵$A$ 的有效估计，进而分离出需要的目标信号。
JADE算法的流程图如下：

下面是JADE算法的公式推导，从论文中截的图

JADE仿真程序

JADE算法的函数：

function [A,S]=jade(X,m) 
% Source separation of complex signals with JADE. 
% Jade performs `Source Separation' in the following sense: 
%   X is an n x T data matrix assumed modelled as X = A S + N where 
%  
% o A is an unknown n x m matrix with full rank. 
% o S is a m x T data matrix (source signals) with the properties 
%    	a) for each t, the components of S(:,t) are statistically 
%    	   independent 
% 	b) for each p, the S(p,:) is the realization of a zero-mean 
% 	   `source signal'. 
% 	c) At most one of these processes has a vanishing 4th-order 
% 	   cumulant. 
% o  N is a n x T matrix. It is a realization of a spatially white 
%    Gaussian noise, i.e. Cov(X) = sigma*eye(n) with unknown variance 
%    sigma.  This is probably better than no modeling at all... 
% 
% Jade performs source separation via a  
% Joint Approximate Diagonalization of Eigen-matrices.   
% 
% THIS VERSION ASSUMES ZERO-MEAN SIGNALS 
% 
% Input : 
%   * X: Each column of X is a sample from the n sensors 
%   * m: m is an optional argument for the number of sources. 
%     If ommited, JADE assumes as many sources as sensors. 
% 
% Output : 
%    * A is an n x m estimate of the mixing matrix 
%    * S is an m x T naive (ie pinv(A)*X)  estimate of the source signals 
[n,T]	= size(X); 
 
%%  source detection not implemented yet ! 
if nargin==1, m=n ; end; 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% A few parameters that could be adjusted 
nem	= m;		% number of eigen-matrices to be diagonalized 
seuil	= 1/sqrt(T)/100;% a statistical threshold for stopping joint diag 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% whitening 
% 
if m<n, %assumes white noise 
 	[U,D] 	= eig((X*X')/T);  
	[puiss,k]=sort(diag(D)); 
 	ibl 	= sqrt(puiss(n-m+1:n)-mean(puiss(1:n-m))); 
 	bl 	= ones(m,1) ./ ibl ; 
 	W	= diag(bl)*U(1:n,k(n-m+1:n))'; 
 	IW 	= U(1:n,k(n-m+1:n))*diag(ibl); 
else    %assumes no noise 
 	IW 	= sqrtm((X*X')/T); 
 	W	= inv(IW); 
end; 
Y	= W*X; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% Cumulant estimation 
 
 
R	= (Y*Y' )/T ; 
C	= (Y*Y.')/T ; 
 
Yl	= zeros(1,T); 
Ykl	= zeros(1,T); 
Yjkl	= zeros(1,T); 
 
Q	= zeros(m*m*m*m,1) ; 
index	= 1; 
 
for lx = 1:m ; Yl 	= Y(lx,:); 
for kx = 1:m ; Ykl 	= Yl.*conj(Y(kx,:)); 
for jx = 1:m ; Yjkl	= Ykl.*conj(Y(jx,:)); 
for ix = 1:m ;  
	Q(index) = ... 
	(Yjkl * Y(ix,:).')/T -  R(ix,jx)*R(lx,kx) -  R(ix,kx)*R(lx,jx) -  C(ix,lx)*conj(C(jx,kx))  ; 
	index	= index + 1 ; 
end ; 
end ; 
end ; 
end 
 
%% If you prefer to use more memory and less CPU, you may prefer this 
%% code (due to J. Galy of ENSICA) for the estimation the cumulants 
%ones_m = ones(m,1) ;  
%T1 	= kron(ones_m,Y);  
%T2 	= kron(Y,ones_m);   
%TT 	= (T1.* conj(T2)) ; 
%TS 	= (T1 * T2.')/T ; 
%R 	= (Y*Y')/T  ; 
%Q	= (TT*TT')/T - kron(R,ones(m)).*kron(ones(m),conj(R)) - R(:)*R(:)' - TS.*TS' ; 
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%computation and reshaping of the significant eigen matrices 
 
[U,D]	= eig(reshape(Q,m*m,m*m));  
[la,K]	= sort(abs(diag(D))); 
 
%% reshaping the most (there are `nem' of them) significant eigenmatrice 
M	= zeros(m,nem*m);	% array to hold the significant eigen-matrices 
Z	= zeros(m)	; % buffer 
h	= m*m; 
for u=1:m:nem*m,  
	Z(:) 		= U(:,K(h)); 
	M(:,u:u+m-1)	= la(h)*Z; 
	h		= h-1;  
end; 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% joint approximate diagonalization of the eigen-matrices 
 
 
%% Better declare the variables used in the loop : 
B 	= [ 1 0 0 ; 0 1 1 ; 0 -i i ] ; 
Bt	= B' ; 
Ip	= zeros(1,nem) ; 
Iq	= zeros(1,nem) ; 
g	= zeros(3,nem) ; 
G	= zeros(2,2) ; 
vcp	= zeros(3,3); 
D	= zeros(3,3); 
la	= zeros(3,1); 
K	= zeros(3,3); 
angles	= zeros(3,1); 
pair	= zeros(1,2); 
c	= 0 ; 
s	= 0 ; 
 
 
%init; 
encore	= 1; 
V	= eye(m);  
 
% Main loop 
while encore, encore=0; 
 for p=1:m-1, 
  for q=p+1:m, 
 
 	Ip = p:m:nem*m ; 
	Iq = q:m:nem*m ; 
 
	% Computing the Givens angles 
 	g	= [ M(p,Ip)-M(q,Iq)  ; M(p,Iq) ; M(q,Ip) ] ;  
 	[vcp,D] = eig(real(B*(g*g')*Bt)); 
	[la, K]	= sort(diag(D)); 
 	angles	= vcp(:,K(3)); 
	if angles(1)<0 , angles= -angles ; end ; 
 	c	= sqrt(0.5+angles(1)/2); 
 	s	= 0.5*(angles(2)-j*angles(3))/c;  
 
 	if abs(s)>seuil, %%% updates matrices M and V by a Givens rotation 
	 	encore 		= 1 ; 
		pair 		= [p;q] ; 
 		G 		= [ c -conj(s) ; s c ] ; 
		V(:,pair) 	= V(:,pair)*G ; 
	 	M(pair,:)	= G' * M(pair,:) ; 
		M(:,[Ip Iq]) 	= [ c*M(:,Ip)+s*M(:,Iq) -conj(s)*M(:,Ip)+c*M(:,Iq) ] ; 
 	end%% if 
  end%% q loop 
 end%% p loop 
end%% while 
 
%%%estimation of the mixing matrix and signal separation 
A	= IW*V; 
S	= V'*Y ; 
 
return ; 

主程序：

%% JADE算法仿真
% 输入信号为两段语音，混合矩阵为随机数构成，
% 采用基于四阶累计量的特征矩阵联合近似对角化JADE算法对两段语音进行分离，并绘制了源信号、混合信号和分离信号
% Author:huasir 2023.9.19 Beijing
close all,clear all;clc;
%=========================================================================%
%                          读取语音文件，输入源信号                       %
%=========================================================================%
[S1,fs1] = audioread('E:\sound1.wav'); % 读取原始语音信号，需要将两个语音文件放置在相应目录下
[S2,fs2] = audioread('E:\ICA\sound2.wav');
figure;
subplot(3,2,1),plot(S1),title('输入信号1'); %绘制源信号
subplot(3,2,2),plot(S2),title('输入信号2');
s1 = S1'; %一行代表一个信号
s2 = S2';
S=[s1;s2];  % 将其组成矩阵
%=========================================================================%
%                      对源信号进行混合，得到观测信号                     %
%=========================================================================%
Sweight = rand(size(S,1));  %由随机数构成混合矩阵
MixedS=Sweight*S;     % 将混合矩阵重新排列
subplot(3,2,3),plot(MixedS(1,:)),title('混合信号1'); %绘制混合信号
subplot(3,2,4),plot(MixedS(2,:)),title('混合信号2');
%=========================================================================%
%               采用JADE算法进行盲源分离，得到源信号的估计                %
%=========================================================================%
[Ae,Se]=jade(MixedS,2);  %Ae为估计的混合矩阵，Se为估计的源信号
% 将混合矩阵重新排列并输出
subplot(3,2,5),plot(Se(1,:)),title('JADE解混信号1');
subplot(3,2,6),plot(Se(2,:)),title('JADE解混信号2');
%=========================================================================%
%        源信号、混合信号以及解混合之后的信号的播放                       %
%=========================================================================%
% sound(S1,8000); %播放输入信号1
% sound(S2,8000); %播放输入信号2
% sound(MixedS(1,:),8000); %播放混合信号1
% sound(MixedS(2,:),8000); %播放混合信号2
% sound(Se(1,:),8000); %播放分离信号1
% sound(Se(2,:),8000); %播放分离信号2
fprintf('混合矩阵为：\n'); % 输出混合矩阵以及估计的混合矩阵
disp(Sweight);
fprintf('估计的混合矩阵为：\n');
disp(Ae);

然后对其进行混合，混合后调用JADE函数进行解混合，最后对解混合的信号进行绘制并进行读取。
可以听到两段录音的内容不一样，音调也不用，它们满足不相关性，因此能够很好的分离。由下图可以看出，分离后的信号的幅度和真实信号有所不同，并且排序也不同，这是盲分离算法本身的局限性：即幅度模糊性和排序模糊性。但是一般情况下，信号的信息保存在波形的变化中，人们对于其绝对幅度并不敏感。
结果如下图：

图1. JADE算法分离结果 在主程序中，首先是读取语音文件，语音文件由以下链接给出，当然也可以自己生成源信号。

链接：https://pan.baidu.com/s/1DwnZqDBc1sogERcq7RrVqA
提取码：ngk1