ICA用于盲源分离.docx
《ICA用于盲源分离.docx》由会员分享，可在线阅读，更多相关《ICA用于盲源分离.docx（18页珍藏版）》请在冰豆网上搜索。
ICA用于盲源分离
盲信号实验报告
ICA算法
******
系别：
电信学院
专业：
电磁场与微波
学号：
**********
*******
2011年07月9日
ICA用于盲源分离
1.ICA模型
设无噪声信号模型为X=As
（1）
A为信号混合矩阵,x是N维观测信号向量,s是M（N>M）维原始信号向量。
由
（1）可见，信号S放大k倍与A的相应列缩小k倍的结果相同,从而决定了ICA得到的信号存在强度的不确定性。
为此,在求解时往往把观测信号先转化为有单位协方差的信号,即在ICA之前先有一个白化过程。
设信号向量y的联合概率密度为p（y），而每一个信号成分的概率密度为p（yi），则信号向量的互信息可以表示为：
（2）
当各个信号成份相互独立时,p（y）=∏Mi=1p（yi）（3）
则I（y）=0
ICA的目的是：
在我们不知道混合矩阵的情况下,寻找线性映射w，从观测信号中提取不能被直接观测的原信号, 这里把它记为:
y=wx=wAs（4）
2.ICA理论、
（1）互信息极小判据
（5）
互信息极小简化成了四阶累积量最大,从而可以通过对四阶累积量的计算,实现独立成份的分离。
（2）信息极大判据
理论分析表明,如果把完成ICA的过程用一个运算网络表示，并在此网络的输出端,引入相应的信源的累积分布函数为变换函数的一个非线性环节,把转化为,则的熵最大就等效于式（5）互信息极小
（3）极大似然估计判据
极大似然估计的目的是通过对观测模型式x=As进行估计，得到潜在的信号S,利用
（6）
当N足够大时，其对数似然概率收敛于它的期望
（7）
式可改写为：
（8）
即，可通过极大似然估计判据提取独立成份
3.常用ICA算法
（1）成对旋转法：
利用Givens旋转,将中的独立成份两两成对旋转直到独立性判据目标函数收敛为止
2）固定点算法:
（i）（四阶累积量）
（9）
（ii）Newton法
（10）
（3）自然梯度学习算法
（11）
4.ICA算法仿真实例：
原始信号混合信号ICA分离信号
5.ICA理论研究方向
（1）固定点算法
（2）非线性ICA研究
（3）噪声ICA研究
（4）overcompleteICA研究
（5）子空间ICA研究
6.ICA应用研究
（1）信号分离
（2）图像去噪
（3）EEG数据分离
（4）fMRI数据处理
下面介绍两种算法进行ICA用于信号分离的实验验证：
第一种：
快速ICA算法
算法步骤：
（1）取任意初始矢量
，要求它满足
（令
=0）。
（2）求
，注意其中k为迭代序号，不是时间序列。
E（.）可以通过对z的各采样时刻值求均值来估计。
（3）将
归一化：
。
（4）如果
不接近1，则令k加1，回到步骤
（2）。
否则迭代结束。
输出最终的
作为
。
（5）提取独立分量：
（注：
具体原理可参考教材第99页。
）
实验程序：
%--------------------产生信号------------------------------
t=linspace（0,100,100）;
s1=sin（t/2）;
plot（s1）
s2=exp（-0.2*t）.*cos（pi*t）;
figure;
plot（s2）;
s3=1/1000*t.^2.*cos（t/4）.^2;
figure;
plot（s3）;
%------------------信号进行混合------------------------------
A=[0.1,0.2,0.7;0.2,0.6,0.2;0.3,0.5,0.2];%混合矩阵
s=[s1;s2;s3];
Y=A*s;%进行信号混合
%-----------------去均值------------------------------------
meanValue=mean（Y'）;
Y=Y-diag（meanValue）*ones（size（Y））;
%-----------------球化-------------------------------------
covarianceMatrix=cov（Y'）;
[U,V,W]=svd（covarianceMatrix）;%矩阵的SVD分解
x=diag（V）;
x=x.^（-1/2）;
V=diag（x）;
Z=V*U'*Y;
%---------------绘制的球化后的信号图------------------------
figure
plot（Z（1,:
））;
figure
plot（Z（2,:
））;
figure
plot（Z（3,:
））;
%-------------盲源分离程序---------------------------------
p=zeros（3）;%存储每次得到的u1矢量
fori=1:
3
u1=rand（3,1）;
a=sqrt（sum（u1.^2））;%产生初始矢量
u1=u1./a;
n=0;
m=1;
while（m>1e-15）
u2=zeros（3,1）;
forj=1:
100
y=u1'*Z（:
j）;
u2=u2+y^3*Z（:
j）;
end
u2=u2./100;%得到迭代后的u2矢量
w3=zeros（3,1）;
forl=1:
i-1
w3=w3+（u2'*p（:
l））*p（:
l）;
end
u2=u2-w3;%进行矢量的施密特正交化，得到相互正交的向量
d=sqrt（sum（u2.^2））;
u2=u2./d;
m=abs（u1'*u2-1）;
u1=u2;
n=n+1;
end
p（:
i）=u1;
end
q=p'*Z;
%--------------------------绘制分离后的信号--------------------------------------
figure;
plot（q（1,:
））;
figure;
plot（q（2,:
））
figure;
plot（q（3,:
））
运行程序得到的结果为：
解混矩阵p为：
p=
0.00270.9915-0.1302
0.4423-0.1180-0.8891
0.89690.05510.4389
源信号：
混合后的信号：
解混后的信号：
分析：
从上面的结果可以看出，分离出的信号幅度和顺序发生了改变，同源信号几乎一致。
第二种算法：
ICA的共轭下降法
原理：
这种方法以最小化互信息准则作为ICA的判据。
在信息论里，互信息是衡量随机变量间独立性的尺度，用负熵来表示互信息为：
其中J（y）表示y的负熵。
当y中各分量相互独立时，最小化互信息就等价于最大化各分离成分的负熵之和。
负熵可以由
近似。
其中，c为常数，v是具有0均值，单位方差的高斯变量，E表示期望运算，G是非线性非二次的函数，本实验中取
。
所以此问题转化为找出权值矩阵W，使分离出来的估计信号
能使函数
最大。
算法实现：
定义目标函数为：
其中，w为矩阵w的向量。
前面提到的负熵由前面的式子近似的条件是
是0均值单位方差的变量。
对于已经白化的观察信号x，这就相当于限制w的范数为1。
转化为无限制条件的优化问题，根据工程优化的相关知识，可以得到惩罚函数为：
F（w）的梯度为：
其中，
为常数，
为G的导数。
w的计算步骤如下：
1.设定初始点
，置
=
设定迭代次数为n=1000次;
2.若n=1000则停止迭代；否则
，进行正交归一化；
其中步长
为
，（0<
<1）,且
满足条件：
3.
，
4.置k：
=k+1，转2.
算法程序：
%------------------------产生源信号并显示----------------------------------
t=linspace（0,100,100）;
s1=sin（t/2）;
plot（s1）
s2=1/1000*t.^2.*cos（t/4）.^2;
figure;
plot（s2）
s3=exp（-0.2*t）.*cos（pi*t）;
figure;
plot（s3）;
%-------------进行源信号混合---------------------------------------
A=[0.1,0.2,0.7;0.2,0.4,0.4;0.3,0.2,0.5];%混合矩阵
s=[s1;s2;s3];
Y=A*s;
%-----------去均值和球化-------------------------------------------
meanValue=mean（Y'）;
Y=Y-diag（meanValue）*ones（size（Y））;
covarianceMatrix=cov（Y'）;
[U,V,W]=svd（covarianceMatrix）;
x=diag（V）;
x=x.^（-1/2）;
V=diag（x）;
Z=V*U'*Y;%球化后的矩阵
%------------显示球化后的信号-----------------------------
figure;
plot（Z（1,:
））;
figure;
plot（Z（2,:
））;
figure;
plot（Z（3,:
））;
%---------------------盲源分离算法主程序----------------------------
symsyG%定义y和G为符号变量
G=1/4*y^4;
w1=rand（3,1）;
w1=w1./sqrt（sum（w1.^2））;%产生初始的范数为1的随机向量
u=zeros（3,1）;
g=diff（G）;
fori=1:
100
y=w1'*Z（:
i）;
u=u+eval（g）*Z（:
i）;
end
a=0.1;
u=-u./100-a*w1;
d1=-u;%产生初始的搜索方向向量d
p=zeros（3）;
%-----选择初始的步长minb-----------------
fori=1:
3
n=0;
c=a.^[0,1,2,3,4];
e=length（c）;
while（n<=1000）
minb=c
（1）;
w2=w1+minb*d1;
h=0;
fork=1:
100
y=w2'*Z（:
k）;
h=h+eval（G）;
end
minF=-h/100-a/2*（norm（w2）^2-1）;
fork=2:
e
w2=w1+c（k）*d1;
h=0;
forr=1:
100
y=w2'*Z（:
r）;
h=h+eval（G）;
end
h=-h/100-a/2*（norm（w2）^2-1）;
if（hminb=c（k）;
minF=h;
end
end
%-------对迭代求得的向量进行正交归一，进行下一次迭代-------------
w2=w1+minb*d1;%w2为经过迭代求得的向量
w3=zeros（3,1）;
forl=1:
i-1
w3=w3+（w2'*p（:
l））*p（:
l）;
end
w2=w2-w3;%w2进行正交化
w2=w2./sqrt（sum（w2.^2））;
u1=zeros（3,1）;
forj=1:
100
y=w2'*Z（:
j）;
u1=u1+eval（g）*Z（:
j）;
end
u1=-u1./100-a*w2;
t=-norm（u1）^2/（d1'*u）;
b=max（t,0）;
d2=-u1+b*d1;
n=n+1;
d1=d2;
w1=w2;
u=u1;
end
p（:
i）=w1;%把求得满足条件的向量存储在p中
q=w1'*Z
figure;
plot（q）;%绘制解混后的信号
end
运行结果：
求得的解混矩阵为：
p=
0.0072-0.9369-0.3495
0.0778-0.34790.9343
0.99690.0339-0.0704
源信号：
混合球化后的信号：
分离后的信号：
与快速ICA比较，得到的结果比其要好。
ICA的方法很多，这里不意义列举，其应用之一可以做图像分离，当混合图像是几个不相关图像的混合时，可以用ICA进行分离。
升级会员 