ICA负熵万能算法详细说明

ICA 负熵 梯度下降 信号分解重建 MATLAB 稍加修改可以应用自己问题

function fushangST1（）
%% 清空环境
clear;clc;
%% 参数设置
fs=1000;
N=100; % 时间横坐标轴
t=（0:N-1）/fs;
f1=60;
s1=sin（2*pi*f1*t）;size（s1） % 正弦信号
s2=2*square（100*t50）; size（s2） %方波信号
s3=randn（size（t））;size（s3） % 随机信号
%%
%% 主程序
% 输出原始源信号
subplot（331）plot（s1）title（‘s1‘）axis（[0100-44]）; % （331）表示整幅图分为3行3列共9幅小图，1表示序号
subplot（332）plot（s2）title（‘s2‘）axis（[0100-44]）;
subplot（333）plot（s3）title（‘s3‘）axis（[0100-44]）;

% 将其组成矩阵
S=[s1;s2;s3]; size（S）;    % 信号个数即为变量数，因此S是一个变量个数＊采样个数的矩阵，这里3个信号，每个信号100列，则S为3*100
Sweight=rand（size（S1））;  % 取一随机矩阵，作为信号混合的权矩阵
MixedS=Sweight*S;         % 得到三个信号的混合信号矩阵

% 将混合矩阵重新排列并输出
subplot（334）plot（MixedS（1:））title（‘mixs1‘）axis（[0100-44]）;
subplot（335）plot（MixedS（2:））title（‘mixs2‘）axis（[0100-44]）;
subplot（336）plot（MixedS（3:））title（‘mixs3‘）axis（[0100-44]）;
MixedS_bak=MixedS;                         % 将混合后的数据备份，以便在恢复时直接调用

% 标准化，预处理步骤之一
MixedS_mean=zeros（31）;
for i=1:3
    MixedS_mean（i）=mean（MixedS（i:））;
end                                        % 计算MixedS的均值
for i=1:3
    for j=1:size（MixedS2）
        MixedS（ij）=MixedS（ij）-MixedS_mean（i）;
    end
end

% 白化，预处理步骤之一
MixedS_cov=cov（MixedS‘）;                    % cov为求协方差的函数
[ED]=eig（MixedS_cov）;                      % 对信号矩阵的协方差函数进行特征值分解
white=inv（sqrt（D））*（E）‘;                        % Q为白化矩阵
MixedS_white=white*MixedS;                      % MixedS_white为白化后的信号矩阵MixedS为混合信号矩阵
    
% 优化迭代，这里可以用多种方法
X=MixedS_white;                            % 以下算法将对X进行操作X=MixedS_white为白化后的信号矩阵就是Z
[VariableNumSampleNum]=size（X）;
numofIC=VariableNum;                       % 在此应用中，独立元个数等于变量个数信号个数，X=MixedS_white白化后信号矩阵的行数
B=zeros（numofICVariableNum）;             % 初始化列向量w的寄存矩阵B=[b1  b2  ...   bd]，初始化权值W矩阵，n*n，n为信号个数3个信号，3行3列
for r=1:numofIC
    i=1;
    maxIterationsNum=100;               % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数）
    IterationsNum=0;
    b=rand（numofIC1）-.5;                  % 随机设置b初值3行1列
    b=b/norm（b）;                           % 对b标准化 norm（b）:向量元素平方和开根号相当于单位化向量b
    while i<=maxIterationsNum+1
        if i == maxIterationsNum           % 循环结束处理
            fprintf（‘\n第%d分量在%d次迭代内并不收敛。‘ rmaxIterationsNum）;
            break;
        end
        bb=b;  % 3*1                       
        a2=1;
        u=1;
        t=X‘*b; % 3*1
        g=t.*exp（-a2*t.^2/2）;%相当于g=y.e（-y2/2）
        dg=（1-a2*t.^2）.*exp（-a2*t.^2/2）;%对g的求导
        b=（（1-u）*t‘*g*b+u*X*g）/SampleNum-mean（dg）*b;%对应公式3-11
                                           % 核心公式，参见理论部分公式2.52
        b=b-B*B‘*b;                        % 对b正交化
        b=b/norm（b）; 
        if abs（abs（b‘*bb）-1）<1e-19        % 如果收敛，则
             B（:r）=b;                     % 保存所得向量b
             break;
         end
        i=i+1;        
    end
%    B（:r）=b;                                % 保存所得向量b
end

%%  信号重建，意思是将ICAedS与原始源信号S比较，会有差异，但是差异越小越好
ICAedS=B‘*white*MixedS_bak;      % 计算ICA后的矩阵

% 将混合矩阵重新排列并输出
subplot（33

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3707  2013-09-24 13:03  fushangST1.m