盲源分离【文档+程序】

【实例简介】参考哈尔滨工程大学2009-柴子娟博士论文《复杂电磁环境下雷达信号的分选识别技术研究》进行的相关仿真，包含两个仿真，一是已知多通道数先验下的盲源分离，二是EMD重构多通道条件下的盲源分离，信号样式为四类脉冲体制雷达信号和高斯白噪声信号

【实例截图】

【源码目录】

盲源分离

├── 盲源分离-仿真测试.pdf
└── 混合调制识别及信号分离程序
    ├── ICA.m
    ├── ICAfs.m
    ├── ICAtest.m
    ├── informaxICA.m
    └── lw-cjf
        ├── lw_cjf_modtest1.m
        ├── lw_cjf_modtest12.m
        ├── lw_cjf_modtest2.m
        ├── lw_cjf_modtest22.m
        └── lw_cjf_modtest3.m

2 directories, 10 files

【核心代码】

function [Z,W] = ICAfs(X,num)

%function:基于负熵固定点算法的逐步剥离步骤
%parameter:X为输入混合信号矩阵，每行包含N个采样点，行数应该大于等于num且秩为num-正定
%正定或超定时：num为源个数，M〉=num--若考虑分离出噪声，则秩为num 1．需提前进行源估计
%实际情况下欠定如何处理？
%reference：2009D-复杂电磁环境下雷达信号的分选识别技术研究-柴娟芳-哈工程
%y^3优化函数更实用于test 1-------牛顿迭代法终止条件

%1-去均值
[M,T] = size(X);      %获取输入矩阵的行列数，行数为观测数据的数目，列数为采样点数
average = mean(X')';  %均值
for i=1:M
    X(i,:) = X(i,:)-average(i)*ones(1,T);
end
 
%1-白化/球化
Cx = cov(X',1);  %计算协方差矩阵Cx
[eigvector,eigvalue] = eig(Cx);   %计算Cx的特征值和特征向量
W = eigvalue^(-1/2)*eigvector';   %白化矩阵
% W=(eigvector*(eigvalue^(-1/2)))';
Z = W*X;    %正交矩阵
 
%2-独立分量数num，p=1,逐个源分离
p = 1;

%迭代-俩个终止条件
Maxcount = 10000;  %最大迭代次数
Critical = 0.00001;  %判断是否收敛
m = num;
W = rand(m);%初始化分离矩阵
for n=1:m
    WP = W(:,n);  %初始权矢量（任意）
    %Y=WP'*Z;
    %G=Y.^3;%G为非线性函数，可取y^3/tanh(ay)&1<a<2、ye(-y^2/2)
    %GG=3*Y.^2；   %G的导数
    count = 0;
    LastWP = zeros(m,1);%初始残差为0-牛顿迭代法终止条件
    W(:,n) = W(:,n)/norm(W(:,n));%欧式距离-按列归一化
    while abs(WP-LastWP)&abs(WP LastWP)>Critical
        count = count 1;  %迭代次数
        LastWP = WP;      %上次迭代的值
        %WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP;
        for i=1:m
%              WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);%tanh(ay)-why LastWP(i)未入均值计算
               WP(i)=mean(Z(i,:).*(((LastWP)'*Z)).^3)-mean(3*((((LastWP))'*Z).^2)).*LastWP(i);%y^3-函数选择决定算法-峭度或者定点迭代
%              WP(i)=mean(Z(i,:).*(((LastWP)'*Z)).^3)-mean(3*((((LastWP))'*Z).^2).*LastWP(i));%y^3
        end
        WPP=zeros(m,1);
        for j=1:n-1
            WPP=WPP (WP'*W(:,j))*W(:,j);
        end
        WP = WP-WPP;
        WP = WP/(norm(WP));
        if count==Maxcount
            fprintf('未找到相应的信号');
            return;
        end
    end
    W(:,n) = WP;
end
     Z = W'*Z;
%     Z = W'*Z W'*average*ones(1,T);%加上均值-误差较大
end