接着上一节PCA的内容，我们来谈一下ICA为何同样重要。
目录
- 引入：鸡尾酒会问题
- 独立成分分析法
- 中心极限定理
- 中心极限定理的MATLAB代码实现
- 测量非高斯程度(Kurtosis峰度)
- 独立成分分析法处理流程
- 使用FASTICA找出独立成分
- MATLAB代码实现
引入：鸡尾酒会问题
我们可以看到如下图所示的鸡尾酒会模型：在一场酒会上，有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:NumofVariableX(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为高斯分布/正态分布，改变式子中的值，图象也会发生变化。
独立成分分析法处理流程
- 均值化(Centering data remove mean)
- 白化(球形化) Whitening process(sphere data)
Decorrelate variables
Scale variables so that their variance = 1 - 优化算法使得每个源的非高斯性最大 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 pathsamplerate=500; % in Hz
N=1000; % data lengthfreq1=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 resultsfigure,
for i=1:size(icasig,1)subplot(size(icasig,1),1,i),plot(icasig(i,:))
endreturn
%% 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的原理，到时候重新再写一篇文章讲一下。
