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