高斯信号.docx
《高斯信号.docx》由会员分享，可在线阅读，更多相关《高斯信号.docx（13页珍藏版）》请在冰点文库上搜索。
高斯信号
亚高斯信号，高斯信号，超高斯信号
一个信号的高斯性是通过其峭度定义的。
在信号x的均值为零的条件下，其峭度定义如下：
kurt（x）=E{x^4}-3[E{x^2}]^2
<0次高斯信号
kurt（x）=0 高斯信号
>0超高斯信号
当我们拿到任意信号x的一个样本后,可通过如下的计算求其峭度，进而判断高斯性：
假设x是1*N的行向量
x=x-mean（x）*ones（1,N）%去均值
KurtX=mean（x.^4）-3*（mean（x.^2））^2%求峭度
均匀分布的信号是次高斯信号，拉普拉斯分布的信号是超高斯信号。
语音信号是超高斯信号。
根据中心极限定理的意义，N个不同分布信号的联合分布有高斯化的趋势，所以信号的非高斯性是盲分离一个很好的优化判据。
另外，在脑电信号处理中，由于一般脑电信号、肌电伪差及工频干扰等多属于亚高斯信号,而心电伪差和眨眼伪差等则多属于超高斯信号。
基本原理是独立元分析（ICA）,它被看作是主元分析（PCA）的扩展。
其差别在于,主元分析只利用了观察数据的二阶统计量且要求各分量正交,而独立元分析利用了观察数据的高阶统计量且要求各分量统计独立。
对于高斯分布的信号,二阶统计量足以描述其特性,但是对于通信系统中典型的通信信号,其分布通常是欠高斯的,所以二阶统计量不足以描述其特性,必须用高阶统计量描述其特性。
系统的基本方程是：
解混的判据是使辅助输出
的信息熵
极大。
%读入混合前的原始图片并显示
I1=imread（'input1.jpg'）;
I2=imread（'input2.jpg'）;
I3=imread（'input3.jpg'）;
I1=rgb2gray（I1）;
I2=rgb2gray（I2）;
I3=rgb2gray（I3）;
subplot（4,3,1）,imshow（I1）,title（'输入图像1'）,
subplot（4,3,2）,imshow（I2）,title（'输入图像2'）,
subplot（4,3,3）,imshow（I3）,title（'输入图像3'）,
%计算图片数据的维数
[x1,y1,z1]=size（I1）;
[x2,y2,z2]=size（I2）;
[x3,y3,z3]=size（I3）;
%将其重新排列为一维行向量并组成矩阵
s1=[reshape（I1,[1,x1*y1*z1]）];
s2=[reshape（I2,[1,x2*y2*z2]）];
s3=[reshape（I3,[1,x3*y3*z3]）];
S_all=[s1;s2;s3];%图片个数即为变量数，图片的像素数即为采样数
%因此S_all是一个变量个数＊采样个数的矩阵
S=double（S_all）;%将图像数据转换为双精度格式
Sweight=rand（size（S_all,1））;%取一随机矩阵，作为信号混合的权矩阵
MixedS=Sweight*S;%得到三个图像的混合信号矩阵
%将混合矩阵重新排列为原始的图片矩阵并输出
ms1=reshape（MixedS（1,:
）,[x1,y1,z1]）;
ms2=reshape（MixedS（2,:
）,[x2,y2,z2]）;
ms3=reshape（MixedS（3,:
）,[x3,y3,z3]）;
I1_mixed=uint8（round（ms1））;
I2_mixed=uint8（round（ms2））;
I3_mixed=uint8（round（ms3））;
subplot（4,3,4）,imshow（I1_mixed）,title（'混合图像1'）,
subplot（4,3,5）,imshow（I2_mixed）,title（'混合图像2'）,
subplot（4,3,6）,imshow（I3_mixed）,title（'混合图像3'）,
MixedS_bak=MixedS;%将混合后的数据备份，以便在恢复时直接调用
%%%%%%%%%%%%%%%%%%%%%%%%%%标准化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_mean=zeros（3,1）;
fori=1:
3
MixedS_mean（i）=mean（MixedS（i,:
））;
end%计算MixedS的均值与方差
fori=1:
3
forj=1:
size（MixedS,2）
MixedS（i,j）=MixedS（i,j）-MixedS_mean（i）;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%白化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_cov=cov（MixedS'）;%cov为求协方差的函数
[E,D]=eig（MixedS_cov）;%对图片矩阵的协方差函数进行特征值分解
Q=inv（sqrt（D））*（E）';%Q为白化矩阵
MixedS_white=Q*MixedS;%MixedS_white为白化后的图片矩阵
IsI=cov（MixedS_white'）;%IsI应为单位阵
%%%%%%%%%%%%%%%%%%%%%%%% FASTICA算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=MixedS_white;%以下算法将对X进行操作
[VariableNum,SampleNum]=size（X）;
numofIC=VariableNum;%在此应用中，独立元个数等于变量个数
B=zeros（numofIC,VariableNum）;%初始化列向量b的寄存矩阵,B=[b1,b2,...,bd]
forr=1:
numofIC%迭代求取每一个独立元
i=1;maxIterationsNum=150;%设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数）
b=2*（rand（numofIC,1）-.5）;%随机设置b初值
b=b/norm（b）;%对b标准化
whilei=maxIterationsNum+1
ifi==maxIterationsNum%循环结束处理
fprintf（'\n第%d分量在%d次迭代内并不收敛。
',r,maxIterationsNum）;
break;
end
bOld=b;%初始化前一步b的寄存器
u=1;
t=X'*b;
g=t.^3;
dg=3*t.^2;
b=（（1-u）*t'*g*b+u*X*g）/SampleNum-mean（dg）*b;
%核心公式，参见理论部分公式2.52
b=b-B*B'*b;%对b正交化
b=b/norm（b）;
ifabs（abs（b'*bOld）-1）<1e-9%如果收敛，则
B（:
r）=b;%保存所得向量b
break;
end
i=i+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%数据复原并构图%%%%%%%%%%%%%%%%%%%%%%%%%%
ICAedS=B'*Q*MixedS_bak;%参见理论部分公式2.55
ICAedS_bak=ICAedS;
ICAedS=abs（55*ICAedS）;
%t1=ICAedS（1,:
）;
%t2=ICAedS（2,:
）;
%t3=ICAedS（3,:
）;
%
%fori=1:
size（t1,2）
%t1（1,i）=t1（1,i）+abs（min（t1））;
%t1（1,i）=t1（1,i）*256/（max（t1）-min（t1））;
%t2（1,i）=t2（1,i）+abs（min（t2））;
%t2（1,i）=t2（1,i）*256/（max（t2）-min（t2））;
%t3（1,i）=t3（1,i）+abs（min（t3））;
%t3（1,i）=t3（1,i）*256/（max（t3）-min（t3））;
%
%end
%
%ICAedS（1,:
）=t1;
%ICAedS（2,:
）=t2;
%ICAedS（3,:
）=t3;
%将计算后的混合矩阵重新排列为图片矩阵并输出
is1=reshape（ICAedS（1,:
）,[x1,y1,z1]）;
is2=reshape（ICAedS（2,:
）,[x2,y2,z2]）;
is3=reshape（ICAedS（3,:
）,[x3,y3,z3]）;
I1_icaed=uint8（round（is1））;
I2_icaed=uint8（round（is2））;
I3_icaed=uint8（round（is3））;
subplot（4,3,7）,imshow（I1_icaed）,title（'ICA解混图像1'）,
subplot（4,3,8）,imshow（I2_icaed）,title（'ICA解混图像2'）,
subplot（4,3,9）,imshow（I3_icaed）,title（'ICA解混图像3'）,
%%%%%%%%%%%%%%%%%%%%%%%%%%PCA计算并构图%%%%%%%%%%%%%%%%%%%%%%%%%%
[V,D]=eig（MixedS_cov）;
Vtmp=zeros（size（V,1）,1）;
forj=1:
2
fori=1:
2
ifD（i,i）tmp=D（i,i）;Vtmp=V（:
i）;
D（i,i）=D（i+1,i+1）;V（:
i）=V（:
i+1）;
D（i+1,i+1）=tmp;V（:
i+1）=Vtmp;
end
end
end
t1=（MixedS'*V（:
1））';
t2=（MixedS'*V（:
2））';
t3=（MixedS'*V（:
3））';
fori=1:
size（t1,2）%亮度调整
t1（1,i）=t1（1,i）+abs（min（t1））;
t1（1,i）=t1（1,i）*256/（max（t1）-min（t1））;
t2（1,i）=t2（1,i）+abs（min（t2））;
t2（1,i）=t2（1,i）*256/（max（t2）-min（t2））;
t3（1,i）=t3（1,i）+abs（min（t3））;
t3（1,i）=t3（1,i）*256/（max（t3）-min（t3））;
end
%%%%%%%%%%%%%%%%%%%%%%%构图%%%%%%%%%%%%%%%%%%%%%%
ts1=reshape（t1,[x1,y1,z1]）;
ts2=reshape（t2,[x2,y2,z2]）;
ts3=reshape（t3,[x3,y3,z3]）;
T1_icaed=uint8（round（ts1））;
T2_icaed=uint8（round（ts2））;
T3_icaed=uint8（round（ts3））;
subplot（4,3,10）,imshow（T1_icaed）,title（'PCA解混图像1'）,
subplot（4,3,11）,imshow（T2_icaed）,title（'PCA解混图像2'）,
subplot（4,3,12）,imshow（T3_icaed）,title（'PCA解混图像3'）,
%%%%%%%%%%%%%%%%%%%%%%%%%%输出一维矢量图%%%%%%%%%%%%%%%%%%%%%%%%%%
%subplot（3,3,1）,plot（S（1,:
））,title（'输入图像1一维矢量图'）
%subplot（3,3,2）,plot（S（2,:
））,title（'输入图像2一维矢量图'）
%subplot（3,3,3）,plot（S（3,:
））,title（'输入图像3一维矢量图'）
%subplot（3,3,4）,plot（ICAedS（1,:
））,title（'ICA解混图像1一维矢量图'）
%subplot（3,3,5）,plot（ICAedS（2,:
））,title（'ICA解混图像2一维矢量图'）
%subplot（3,3,6）,plot（ICAedS（3,:
））,title（'ICA解混图像3一维矢量图'）
%subplot（3,3,7）,plot（t1）,title（'PCA解混图像1一维矢量图'）
%subplot（3,3,8）,plot（t2）,title（'PCA解混图像2一维矢量图'）
%subplot（3,3,9）,plot（t3）,title（'PCA解混图像3一维矢量图'）
im=imread（imgname）;
im=rgb2gray（im）;
imnoise_1=imnoise（im,'speckle',1）;
imnoise_1=uint8（imnoise_1）;
g2=size（imnoise_1）;
figure
（2）imshow（imnoise_1）;
title（'Noiseimage'）;
imnoise_2=double（imnoise_1）;
imnoise_3=imnoise_2+var;
imnoise_3=log（imnoise_3）;
immoise_3=double（imnoise_3）;
imnoise_4=exp（imnoise_3）;
imnoise_4=imnoise_4-var;
imnoise_5=uint8（imnoise_4）;
figure（5）imshow（imnoise_5）title（'recoveryspeckleimage'）;
升级会员 