demo.m

完整的独立分量分析(ICA)程序包,可以用来进行盲源分离

% demo --- using FastICA to separating 3 mixed images
%PCA-ICA方法实现SAR图像中的相干斑去除
%FastICA算法：收敛速度快但稳定性不够好



clear;clc;close all;

% 获取三个大小为256×256的双精度(0-1)图象
%-------------------------- case 1 begin ---------------------------------------
% 第一幅图象：
info = imfinfo('HH.jpg');
Width = info.Width;
Height = info.Height;
I1 = imread('HH.jpg'); 
I1 = im2double(I1); 
figure; imshow(I1,[]); title('原始图象1');


% 第二幅图象：
I2 = imread('HV.jpg');
I2 = im2double(I2); 
figure; imshow(I2,[]);title('原始图象2');


% 第三幅图象：
I3 = imread('VV.jpg'); 
I3 = im2double(I3); 
figure; imshow(I3,[]); title('原始图象3');
% 第四幅图象：
I4 = imread('hav.jpg'); 
I4 = im2double(I3); 
% 第五幅图象：
I5 = imread('hsv.jpg'); 
I5 = im2double(I3); 
%----------------------------- case 1 end ---------------------------------------

% 把图象按行读取，存贮在行向量中，并求其峭度值
imv1 = m1'; kur1=kurtosis(imv1) 
imv2 = m2'; kur2=kurtosis(imv2) 
imv3 = m3'; kur3=kurtosis(imv3) 
imv4 = m4'; kur4=kurtosis(imv4) 

% 构造原始图象数据矩阵，其中每行为一个图象数组
S = [imv1;imv2;imv3;imv4];

imv1=eqim(S(1,:));%归一化
imv2=eqim(S(2,:));
imv3=eqim(S(3,:));
imv4=eqim(S(4,:));



Seqim=[imv1;imv2;imv3;imv4];

[S, mixedmean] = remmean(Seqim);%去均值

covarianceMatrix0=cov(S');          % 产生协方差阵
[E,D]=eig(covarianceMatrix0);         % 求出协方差的特征向量和特征值
% col1=E(:,1);
col2=E(:,2);
col3=E(:,3);
col4=E(:,4);
E1=[col2';col3';col4']';
vec=[D(2,2),D(3,3),D(4,4)];
D1=diag(vec);


% 对数据进行白化预处理

[Z, V] = whiten(S, E1, D1);




% 采用快速固定点算法分离图象
[Y,B] = fastICA(Z,'nonlinfunc','tanh','OverValue',0.000001);
   B
   W=B'*V
   Y= W *S+ (W * mixedmean) * ones(1, 20848);

[tempY,m0]=remmean(Y);


re1 = Y(1,:); 
re2 = Y(2,:); 
re3 = Y(3,:); 
re1=(re1-min(re1))/(max(re1)-min(re1));
re2=(re2-min(re2))/(max(re2)-min(re2));
re3=(re3-min(re3))/(max(re3)-min(re3));
re1=(re1-0.5)*2;
re2=(re2-0.5)*2;
re3=(re3-0.5)*2;
wavwrite(re1,'re1.wav');
wavwrite(re2,'re2.wav');
wavwrite(re3,'re3.wav');

% figure; imshow(im1,[]); 
% figure; imshow(im2,[]); 
% figure; imshow(im3,[]); 
%  
% imwrite(im1,'IC1.jpg');
% imwrite(im2,'IC2.jpg');
% imwrite(im3,'IC3.jpg');
YOut=logeqim(Y);
%[re1CR,re1ENL]=ContrastRatio((YOut(1,:))');
%[re2CR,re2ENL]=ContrastRatio((YOut(2,:))');
fprintf('-------------------------The perfomance index of images after ICA--------------------------------\n');
fprintf('img1\t\t\timg2\n');
%fprintf('-------------------------Contrast ratio of images------------------------\n');
%fprintf('%f\t\t\t%f\n',re1CR,re2CR);
%fprintf('-------------------------ENL of images------------------------\n');
%fprintf('%f\t\t\t%f\n',re1ENL,re2ENL);
k1=kurtosis(YOut(1,:))
k2=kurtosis(YOut(2,:))
k3=kurtosis(YOut(3,:))