【医学信号处理与MATLAB（5）】独立成分分析法-ICA

接着上一节PCA的内容，我们来谈一下ICA为何同样重要。

引入：鸡尾酒会问题

我们可以看到如下图所示的鸡尾酒会模型：在一场酒会上，有3个说话者(si)和三个麦克风(xi)，麦克风距离说话者的远近不同，因此收音效果也不同，但实际都可看作是线性混合的，如图右边的式子所示。
将线性方程组用矩阵的方式表示，如图右下方，A为3乘3的混合矩阵。
如何根据麦克风里的内容分别推测出三个说话者的声音呢？这便是ICA需要处理的问题。

独立成分分析法

与PCA不同，目标不在降低资料变量的维度，而是尽可能从混合信号中找出更具有生理或物理意义的信号来源。

假设：

• 信号源s都是各自独立的变量
• 各个信号源变量的分布都是非高斯分布
  (会让复杂的证明稍微简单一点233)

限制：

• 无法计算出信号源s的方差或者能量强度 （如：可求解混合比例，但不知麦克风与人的距离）
• 无法解出信号源s的正负号

中心极限定理

中心极限定理(Central Limit Theorem)：数个非高斯分布的独立变量，其混合信号的分布会趋向高斯分布。

中心极限定理的MATLAB代码实现

clear,close all
clc

%% initialize parameters
NumofVariable=20;
N = 10000;  %data points

%generate rand data
X=[];
for i=1:NumofVariable
    X(i,:) = rand(1,1000)-0.5;    % uniform rand noise 
%     X(i,:) = sin(2*pi*i*[1:1000]/500);   % sine wave
end

%% mix variables and plot histogram
mixX = mean(X,1);
mixX_kurt=kurtosis(mixX)-3;  % help kurtosis for details

[count,bin]=hist(mixX,100);
figure,plot(bin, count/max(count))
title(['kurtosis (non-Gaussianity) = ' num2str(mixX_kurt)])

测量非高斯程度(Kurtosis峰度)

目的：我们想让变量尽量做到非高斯分布，因此需要一个值去测量非高斯程度。

• Kurtosis = 0 → Gaussian
• Kurtosis ＞ 0 → Super-Gaussian
• Kurtosis ＜ 0 → Sub-Gaussian

下图中黑线=0为高斯分布/正态分布，改变式子中的值，图象也会发生变化。

独立成分分析法处理流程

1. 均值化(Centering data remove mean)
2. 白化(球形化) Whitening process(sphere data)
   Decorrelate variables
   Scale variables so that their variance = 1
3. 优化算法使得每个源的非高斯性最大 Optimization algorithm to maximize non-Gaussianity of each source

均值化和PCA一样，就是我们会将数据减去它的均值。
白化通俗点解释就是，我们在PCA里面取得了主成分之后，可以发现主成分的方差其实相差还是很大的，这样也能够让我们取得一些我们需要的成分。但是在ICA中，我们想将主成分的方差变得一样(=1)，比如将方差放大或者缩小等等使其达到一致。这便是球形化的概念

优化算法：在基础阶段，我们会使用FASTICA就足够了。

使用FASTICA找出独立成分

还是非常熟悉的图，我们在上一节PCA中也用到的：

只不过我们这次需要做的，是将左边的正弦信号和三角波信号给分离出来。

MATLAB代码实现

clear, close all
clc

%% initialize parameters
path(path,'.\fastica_25');  % add FastICA path

samplerate=500; % in Hz
N=1000; % data length

freq1=5; % in Hz
freq2=7; % in Hz
taxis=[1:N]/samplerate;

ICNo=3;  % the specified number of independent components
PCNo=3;  % the preserved number of PCs 

%% generate test signals 
C1 = 0.75*sin(2*pi*freq1*taxis);        % 1st component: a sine wave
C2 = sawtooth(2*pi*freq2*taxis,0.5);    % 2nd component: a triangular wave

% Combine data in different proportions
X(1,:) = 0.5*C1 + 0.5*C2 + 0.1*rand(1,N);
X(2,:) = 0.7*C1 + 0.2*C2 + 0.1*rand(1,N);
X(3,:) = 0.2*C1 + 0.7*C2 + 0.1*rand(1,N);
X(4,:) = -0.3*C1 - 0.6*C2 + 0.3*rand(1,N);
X(5,:) = 0.6*rand(1,N);    % Noise only

% Center data by subtracting mean
X = X - mean(X,2)*ones(1,N);

figure, 
for i=1:size(X,1)
    subplot(size(X,1),1,i)
    plot(taxis,X(i,:)),xlim([taxis(1) taxis(end)])
end

%% Independent Components Analysis using FastICA
[icasig, A, W] = fastica(X,'numOfIC',ICNo,'displayMode','off','firstEig',1,'lastEig',PCNo); % fast ICA

% [icasig, A, W] = fastica(X,'numOfIC',ICNo,'displayMode','off','interactivePCA','on'); % fast ICA
% fasticag(X)  % used the correct number of PCA can always produce good results

figure,
for i=1:size(icasig,1)
    subplot(size(icasig,1),1,i),plot(icasig(i,:))
end

return
%% reconstruct signal  信号的重构（丢掉不需要的成分）
rejectICA=1;
% reconstruct the signal
A2=A;
icasig2=icasig;
A2(:,rejectICA)=[];
icasig2(rejectICA,:)=[];

newX=(A2*icasig2);

figure, 
for i=1:size(newX,1)
    subplot(size(newX,1),1,i)
    plot(taxis,newX(i,:)),xlim([taxis(1) taxis(end)])
end

关于ICA的原理，到时候重新再写一篇文章讲一下。