盲源分离（BSS）的学习总结（PCA、ICA）

盲源分离BSS(Blind Signal Separation)

从多维观测信号中分离出源信号，除去混叠与噪声的过程。可以用于麦克风阵列的信号分析、生理电信号（EEG）等多输入多输出的采集场景，多数据指标融合分析。一般观测的通道数M会大于信号源R的数目，信号处理效果较好，称其超定模型。反之，为欠定。

主要的分析方法：

1.主成分分析（PCA）

2.独立成分分析（ICA）

1.主成分分析

PCA本质是一个降维的过程，目的寻找变量的冗余度更小的子集。PCA是ICA的一个非常必要的预处理：可以降低信号的维度，ICA要求输入维度与输出维度一致；当通道数等于源信号数，PCA做预处理，可以将ICA的估计参数减少一半，降低计算复杂度。

但PCA（或SVD）方法保证分解出来的信号分量不相关，不能保证这些分量独立。而ICA就是针对独立分量分解，盲源分离与独立分量分析有很强的联系，因此常采用ICA作为BSS的分析方法。

参考文献链接：CodingLabs - PCA的数学原理、http://blog.codinglabs.org/articles/pca-tutorial.html(54条消息) PCA的数学原理（经典，正常公式版）_Tele_Anti_Nomy的博客-CSDN博客https://blog.csdn.net/Tele_Anti_nomy/article/details/79676017
PCA算法

总结一下PCA的算法步骤：

设有x=m*n数据，降至m*k维数据。

1）将x的每一列进行零均值化

2）求出协方差矩阵

3）求出协方差矩阵的特征值及对应的特征向量

4）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

5）即为降维到k维后的数据

用一个实例很好的说明这个算法的含义

为例，我们用PCA方法将这组二维数据其降到一维。

因为这个矩阵的每行已经是零均值，这里我们直接求协方差矩阵：

然后求其特征值和特征向量，具体求解方法不再详述，可以参考相关资料。求解后特征值为：

其对应的特征向量分别是：

其中对应的特征向量分别是一个通解，c1和c2可取任意实数。那么标准化后的特征向量为：

因此我们的矩阵P是：

可以验证协方差矩阵C的对角化：

最后我们用P的第一行乘以数据矩阵，就得到了降维后的表示：

这里原文的计算Y有错误,实际

Y= $\left ( -3\sqrt{2}\; -1/\sqrt{2}\; \; 0\; \; 3/\sqrt{2}\; \; 1/\sqrt{2}\right )$

降维投影结果如下图：

那么原先的点是一个二维点再蓝色线上的投影之后，变成以蓝色坐标轴的一维向量，其意义就是由二维空间点变成一维长度点。这跟SVD奇异值分解类似。

Matlab 算法:

1.根据上面的流程算法编写Matlab算法，参考：https://www.pianshen.com/article/6398423334/

算法1：


%% 导入上述数据
clc
clear all
FileName1 = ['C:\Users\Administrator\Desktop\PCA_BSS\pca.xlsx'];
A=xlsread(FileName1);
%% 数据标准化处理
a = size(A,1);
b = size(A,2);
for i = 1:bSA(:,i) = (A(:,i) - mean(A(:,i)))/std(A(:,i));
end%% 计算相关系数矩阵的特征值和特征向量
CM = corrcoef(SA);                          %计算相关系数矩阵
[V,D] = eig(CM);                            %计算特征值和特征向量
for j = 1:bDS(j,1)=D(b+1-j,b+1-j);                 %对特征值按降序排列
endfor i = 1:bDS(i,2) = DS(i,1)/sum(DS(:,1));         %贡献率DS(i,3) = sum(DS(1:i,1))/sum(DS(:,1));  %累计贡献率
end%% 选择主成分及对应的特征向量
T = 0.9;    %主成分保留率
for K = 1:bif DS(K,3) >= TCom_num = K;breakend
end%% 提取主成分对应的特征向量
for j = 1:Com_numPV(:,j)=V(:,b+1-j);
end%% 计算个评价对象的主成分的分
new_score = SA*PV;
for i = 1:atotal_score(i,1)= sum(new_score(i,:));total_score(i,2)= i;
endresult_report = [new_score,total_score];    %将各主成分的分与总分放在同一个举证中
result_report = sortrows(result_report,-4); %将总分降序排列
%% 输出模型及结果报告
disp('主成分的分及排序（按第四列的总分进行降序排列，前3列为各主成分得分，第五列为企业编号）')
result_report

本文作者对上述算法进行简单解析，计算特征值降序，求出累计贡献率大于0.9时需要的维度，案例中k为3，对应3个特征向量，组成PV降维矩阵，因为我们是列的降维，所以SA右乘以PV矩阵得到降维之后的矩阵。按行求和，得出15组数据的评价指标，最终按照得分高低排序(result_report第4列),第5列为1-15组数据的排名

PV =

0.3334 0.3788 0.3115
0.3063 0.5562 0.1871
0.3900 -0.1148 -0.3182
0.3780 -0.3508 0.0888
0.3853 -0.2254 -0.2715

result_report =

5.1936 -0.9793 0.0207 4.2350 9.0000
0.7662 2.6618 0.5437 3.9717 1.0000
1.0203 0.9392 0.4081 2.3677 8.0000
3.3891 -0.6612 -0.7569 1.9710 6.0000
0.0553 0.9176 0.8255 1.7984 5.0000
0.3735 0.8378 -0.1081 1.1033 13.0000
0.4709 -1.5064 1.7882 0.7527 15.0000
0.3471 -0.0592 -0.1197 0.1682 14.0000
0.9709 0.4364 -1.6996 -0.2923 2.0000
-0.3372 -0.6891 0.0188 -1.0075 10.0000
-0.3262 -0.9407 -0.2569 -1.5238 7.0000
-2.2020 -0.1181 0.2656 -2.0545 4.0000
-2.4132 0.2140 -0.3145 -2.5137 11.0000
-2.8818 -0.4350 -0.3267 -3.6435 3.0000
-4.4264 -0.6180 -0.2884 -5.3327 12.0000

算法2：利用Matlab的 pcacov(X)函数

%% 数据导入处理
clc
clear all
FileName1 = ['C:\Users\CNWEWAN35\Desktop\PCA主成分分析\PCA_BSS.xlsx'];
A = xlsread(FileName1);
%% 数据标准化处理
a = size(A,1);
b = size(A,2);
for i = 1:bSA(:,i) = (A(:,i) - mean(A(:,i)))/std(A(:,i));
end
%% 计算相关系数矩阵主成分、贡献率
CM = corrcoef(SA);                          %计算相关系数矩阵
[PC,L,ex]=pcacov(CM)
for i = 1:bEx(i) = sum(ex(1:i,1))/sum(ex);  %累计贡献率
end
%% 选择主成分及对应的特征向量
T = 0.9;    %主成分保留率
for M = 1:bif Ex(M) >= TCom_Num2 = M;breakend
endPC=PC(:,1:M);
PC(:,3)=-PC(:,3); %这里为什么在第三列取相反数呢，因为我们的计算机计算特征值会产生一个符号问题，也是为了方便对比算法一，所以这里对第三列取反.
%% 计算个评价对象的主成分的分new_score2 = SA*PC;
for i = 1:atotal_score2(i,1)= sum(new_score2(i,:));total_score2(i,2)= i;
end
result_report2 = [new_score2,total_score2];    %将各主成分的分与总分放在同一个举证中
result_report2 = sortrows(result_report2,-4)%将总分降序排列

PC =

0.3334 0.3788 0.3115
0.3063 0.5562 0.1871
0.3900 -0.1148 -0.3182
0.3780 -0.3508 0.0888
0.3853 -0.2254 -0.2715
0.3616 -0.4337 0.0696
0.3026 0.4147 -0.6189
0.3596 -0.0031 0.5452

result_report2 =

两种方法结果一致。

2.独立成分分析（ICA）

背景：鸡尾酒会问题：人耳可以在多人说话或者嘈杂的环境中有选择性的接收他感兴趣的说话人的语音，这需要人脑对音频信号的处理，可以说非常复杂，则也是盲源分离思想的初衷。

应用：1.麦克风阵列盲分离 2.生物医学：ECG提取房颤信息，EEG多通道脑电信号 3.图像增强处理、恢复 4.故障诊断 5.雷达传感器融合

介绍三种盲源分离算法：JADE（联合近似对角化法）FastICA(固定点算法)、InforMax(信息量最大化)
(46条消息) ICA算法_heda3的博客-CSDN博客_ica算法

ICA分析与应用 - 百度文库 (baidu.com)

(46条消息) ICA(独立成分分析)在信号盲源分离中的应用_Richard_Cai-CSDN博客_ica盲源分离

几种盲源分离算法的比较 - 道客巴巴 (doc88.com)

做了一个流程图，帮助大家理解：

包括未知的源信号与混合系统，观测信号，分离系统

%主程序
clear all;clc;N=200;n=1:N; % N为采样点数s1=2*sin(0.02*pi*n); % 正弦信号
t=1:N;s2=2*square(100*t,50); % 方波信号
a=linspace(1,-1,25);s3=2*[a,a,a,a,a,a,a,a]; % 锯齿信号
s4=rand(1,N); % 随机噪声
S=[s1;s2;s3;s4]; %四通道源信号
A=rand(4,4);
X=A*S; % 观察信号，四通道观测信号,一般来说观测信号通道数大于等于源信号数，才不会欠分离%源信号波形图
figure(1);subplot(4,1,1);plot(s1);title('源信号');
subplot(4,1,2);plot(s2);
subplot(4,1,3);plot(s3);
subplot(4,1,4);plot(s4);xlabel('Time/ms');%观察信号（混合信号）波形图
figure(2);subplot(4,1,1);plot(X(1,:));title('观察信号（混合信号）');
subplot(4,1,2);plot(X(2,:));
subplot(4,1,3);plot(X(3,:));subplot(4,1,4);plot(X(4,:));Z=ICA(X);figure(3);subplot(4,1,1);plot(Z(1,:));title('解混后的信号');
subplot(4,1,2);plot(Z(2,:));
subplot(4,1,3);plot(Z(3,:));
subplot(4,1,4);plot(Z(4,:));xlabel('Time/ms');

%ICA函数function Z=ICA(X)
%-----------去均值---------
[M,T] = size(X); %获取输入矩阵的行/列数，行数为观测数据的数目，列数为采样点数
average=mean(X')';  %均值
for i=1:MX(i,:)=X(i,:)-average(i)*ones(1,T);
end
%---------白化/球化------
Cx =cov(X',1);    %计算协方差矩阵Cx
[eigvector,eigvalue]= eig(Cx); %计算Cx的特征值和特征向量
W=eigvalue^(-1/2)*eigvector';   %白化矩阵
Z=W*X;   %正交矩阵
%----------迭代-------
Maxcount=10000;        %最大迭代次数
Critical=0.00001;   %判断是否收敛
m=M;                %需要估计的分量的个数
W=rand(m);
for n=1:m WP=W(:,n); %初始权矢量（任意）count=0;LastWP=zeros(m,1);W(:,n)=W(:,n)/norm(W(:,n));while abs(WP-LastWP) & abs(WP+LastWP)>Criticalcount=count+1;   %迭代次数LastWP=WP;      %上次迭代的值for i=1:m    WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);endWPP=zeros(m,1);for j=1:n-1WPP=WPP+(WP'*W(:,j))*W(:,j);endWP=WP-WPP;WP=WP/(norm(WP));       if count==Maxcountfprintf('未找到相应的信号');return;end
endW(:,n)=WP;
end
Z=W'*Z;

上述代码使用四种线性不相关的信号进行混合

混合信号作为观测信号，进行分解后得到解混后的信号（S1,S2,S3,S4，无顺序问题）基本还原了源信号