一、簡介

1 基礎概念
FastICA演算法，又稱不動點(Fixed-Point)演算法，是由芬蘭赫爾辛基大學Hyvärinen等人提出來的。是一種快速尋優迭代演算法，與普通的神經網路演算法不同的是這種演算法採用了批處理的方式，即在每一步迭代中有大量的樣本資料參與運算。但是從分散式並行處理的觀點看該演算法仍可稱之為是一種神經網路演算法。
FastICA演算法有基於峭度、基於似然最大、基於負熵最大等形式，這裡，我們介紹基於負熵最大的FastICA演算法（可以有效地把不動點迭代所帶來的優良演算法特性與負熵所帶來的更好統計特性結合起來）。它以負熵最大作為一個搜尋方向，可以實現順序地提取獨立源，充分體現了投影追蹤（Projection Pursuit）這種傳統線性變換的思想。此外，該演算法採用了定點迭代的優化演算法，使得收斂更加快速、穩健。
因為FastICA演算法以負熵最大作為一個搜尋方向，因此先討論一下負熵判決準則。由資訊理論理論可知：在所有等方差的隨機變數中，高斯變數的熵最大，因而我們可以利用熵來度量非高斯性，常用熵的修正形式，即負熵。根據中心極限定理，若一隨機變數由許多相互獨立的隨機變數之和組成，只要具有有限的均值和方差，則不論其為何種分佈，隨機變數較更接近高斯分佈。換言之，較的非高斯性更強。因此，在分離過程中，可通過對分離結果的非高斯性度量來表示分離結果間的相互獨立性，當非高斯性度量達到最大時，則表明已完成對各獨立分量的分離[1]。

2 負熵的定義：


3 演算法原理

4 演算法步驟

二、原始碼

clc
clear all
%% --------------------------------- Set Parameters
N = 1;                              %The number of observed mixtures
Ns = 2;                             %The number of independent sources
Ls = 1000;                          %Sample size, i.e.: number of observations
finalTime = 40*pi;                  %Final sample time (s)
initialTime = 0;                    %Initial sample time (s)

%% --------------------------------- Generating Data for SSA-ICA
Amix = rand(N,Ns);                      %Amix is a random N x Ns mixing matrix
timeVector = initialTime:(finalTime-initialTime)/(Ls-1):finalTime;  %Vector of time coordinates
source1 = sin(1.1*timeVector);          %Independent source component 1, sin(a * t)
source2 = cos(0.25*timeVector);         %Independent source component 2, cos(b * t)
S = [source1;source2];                  %Source Matrix

figure
plot(timeVector,source1)                    %Plotting the N independent sources vs. time
xlabel('time (s)')
ylabel('Signal Amplitude') 
legend('source 1')

figure
plot(timeVector,source2)                    %Plotting the N independent sources vs. time
xlabel('time (s)')
ylabel('Signal Amplitude') 
legend('source 2')

Yobs = Amix*S;                              %Matrix consisting of M samples of N observed mixtures

figure
plot(timeVector,Yobs)                       %Plotting the observed signal vs. time
xlabel('time (s)')
ylabel('Signal Amplitude') 
legend('observed signal')
function [Sest] = Fast_ICA(Xobs,C)
%% Preprocessing, Centering
SX = size(Xobs);
N = SX(1);
M = SX(2);

X = Xobs';                                  %X is the transpose of the matrix of M samples of N mixtures, used in subsequent calculations

Xmean = mean(X);                            %Xmean is the mean vector of the matrix X

for i = 1:N
    X(:,i) = X(:,i) - Xmean(i);             %The matrix X is centered by subtracting each of the N mixtures by their corresponding sample averages
end

%% Preprocessing, Whitening

ExxT    = cov(X);                           %The covariance matrix of X is computed and stored in ExxT
[E,D]   = eig(ExxT);                        %Eigenvalue decomposition is applied on the covariance matrix of X, ExxT

Z = E*1/sqrt(D)*E'*X';                      %The matrix X is whitened to Z

%% FastICA algorithm

W = 0.5*ones(C,N);                          %Initializing W, a matrix consisting of columns corresponding with the inverse of the (transformed) mixing Amix

iterations = 100;                           %The amount of iterations used in the fastICA algorithm

for p = 1:C
    
wp = ones(N,1)*0.5;
wp = wp / sqrt(wp'*wp);

    for i = 1:iterations
        
        G       = tanh(wp'*Z);
        Gder    = 1-tanh(wp'*Z).^2;
        
        wp      = 1/M*Z*G' - 1/M*Gder*ones(M,1)*wp;
        
        dumsum  = zeros(C,1);
        
        for j = 1:p-1
            dumsum = dumsum + wp'*W(:,j)*W(:,j);
        end
        
        wp      = wp - dumsum;        
        wp      = wp / sqrt(wp'*wp);
    end
    
    W(:,p) = wp; 
end

%% Output Results
W = W/sqrt(2);    %The factor sqrt(2) is an emirical constant added to make the predictions fit the data properly. The source of the factor has yet to be determined.
Sest = W'*Z;
end


 

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423