1、高斯信号亚高斯信号，高斯信号，超高斯信号一个信号的高斯性是通过其峭度定义的。在信号x的均值为零的条件下，其峭度定义如下：kurt(x)=Ex4-3Ex22 0 超高斯信号当我们拿到任意信号x的一个样本后,可通过如下的计算求其峭度，进而判断高斯性：假设x是1*N的行向量x=x-mean(x)*ones(1,N) %去均值KurtX=mean(x.4)-3*(mean(x.2)2 %求峭度 均匀分布的信号是次高斯信号，拉普拉斯分布的信号是超高斯信号。语音信号是超高斯信号。根据中心极限定理的意义，N个不同分布信号的联合分布有高斯化的趋势，所以信号的非高斯性是盲分离一个很好的优化判据。另外，在脑电信号
2、处理中，由于一般脑电信号、肌电伪差及工频干扰等多属于亚高斯信号,而心电伪差和眨眼伪差等则多属于超高斯信号。基本原理是独立元分析(ICA),它被看作是主元分析(PCA)的扩展。其差别在于,主元分析只利用了观察数据的二阶统计量且要求各分量正交,而独立元分析利用了观察数据的高阶统计量且要求各分量统计独立。对于高斯分布的信号,二阶统计量足以描述其特性,但是对于通信系统中典型的通信信号,其分布通常是欠高斯的,所以二阶统计 量不足以描述其特性,必须用高阶统计量描述其特性。 系统的基本方程是： 解混的判据是使辅助输出的信息熵极大。% 读入混合前的原始图片并显示 I1=imread (input1.jpg);
3、 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); % 将其重新排列为一维行向量
4、并组成矩阵 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
5、(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),t
6、itle(混合图像3), MixedS_bak=MixedS; % 将混合后的数据备份，以便在恢复时直接调用 % 标准化 % MixedS_mean=zeros(3,1); for i=1:3 MixedS_mean(i)=mean(MixedS(i,:); end % 计算MixedS的均值与方差 for i=1:3 for j=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); %
7、对图片矩阵的协方差函数进行特征值分解 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
8、 for r=1:numofIC % 迭代求取每一个独立元 i=1;maxIterationsNum=150; % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数） b=2*(rand(numofIC,1)-.5); % 随机设置b初值 b=b/norm(b); % 对b标准化 while i=maxIterationsNum+1 if i = maxIterationsNum % 循环结束处理 fprintf(n第%d分量在%d次迭代内并不收敛。, r,maxIterationsNum); break; end bOld=b; % 初始化前一步b的寄存器 u=1; t=X*b;
9、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); if abs(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,:
10、); % t2=ICAedS(2,:); % t3=ICAedS(3,:); % % for i=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,:)=t
11、1; % 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解混图像
12、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); for j=1:2 for i=1:2 if D(i,i)D(i+1,i+1) 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 en
13、d t1=(MixedS*V(:,1); t2=(MixedS*V(:,2); t3=(MixedS*V(:,3); for i=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
14、=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
15、(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(ICA
16、edS(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(Noise image);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(recovery speckle image);
