negentest_01_defl_gauss.m

A Deflationary fastICA using negentropy approximation implementted in Matlab

close all;clear all;clc;

[v1,fs1,bit1] = wavread('v1.wav');
[s3,fs2,bit2] = wavread('s3.wav');
[s1,fs3,bit3] = wavread('s1.wav');

v1 = v1(1:20000,1);
s3 = s3(1:20000,1);
s1 = s1(1:20000,1);

figure;
subplot(3,1,1);plot(v1,'r');grid; % plot mixing 1
subplot(3,1,2);plot(s3,'b');grid; % plot mixing 2
subplot(3,1,3);plot(s1,'g');grid; % plot mixing 2
suptitle('Original Signals');

M1 = 3.61*v1 - 2*s3 + 1.5*s1;            % mixing 1
M2 = 1.73*v1 - 1.41*s3 + 3.2*s1;            % mixing 2
M3 = 2.3*v1 + 3.41*s3 - 1.6*s1;

figure;
subplot(3,1,1);plot(M1,'r');grid; % plot mixing 1
subplot(3,1,2);plot(M2,'b');grid; % plot mixing 2
subplot(3,1,3);plot(M3,'g');grid; % plot mixing 2
suptitle('Mixed Signals');

X = [M1 M2 M3];
% X = [v1 s3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remove the mean and check the data
% [mixedsig, mixedmean] = remmean(mixedsig);
% function [newVectors, meanValue] = remmean(vectors);

mixedsig = X';
% newVectors = zeros (size (mixedsig));
mixedmean = mean (mixedsig')';
mixedsig = mixedsig - mixedmean * ones (1,size (mixedsig, 2));

% mixedsig  = newVectors;
% mixedmean = meanValue;
% clear newVectors meanValue;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate PCA
% [E, D]=pcamat(mixedsig, firstEig, lastEig, interactivePCA, verbose);
% function [E, D] = pcamat(vectors, firstEig, lastEig, s_interactive, ...
%     s_verbose);

[numSamples,vectorSize] = size(X);

% Calculate the covariance matrix.
% if b_verbose, fprintf ('Calculating covariance...\n'); end
% covarianceMatrix = cov(mixedsig', 1)
covarianceMatrix = 1/(numSamples-1) * mixedsig * mixedsig'
% Calculate the eigenvalues and eigenvectors of covariance
% matrix.
[E, D] = eig (covarianceMatrix);

% Sort the eigenvalues - decending.
eigenvalues = flipud(sort(diag(D)));

% ========================================================
% Calculate the whitening
% [whitesig, whiteningMatrix, dewhiteningMatrix] = whitenv ...
% 						     (mixedsig, E, D, verbose);
% function [newVectors, whiteningMatrix, dewhiteningMatrix] = whitenv ...
%     (vectors, E, D, s_verbose);

% Calculate the whitening and dewhitening matrices (these handle
% dimensionality simultaneously).
whiteningMatrix = inv (sqrt (D)) * E';
dewhiteningMatrix = E * sqrt (D);

% Project to the eigenvectors of the covariance matrix.
% Whiten the samples and reduce dimension simultaneously.
% if b_verbose, fprintf ('Whitening...\n'); end
newVectors =  whiteningMatrix * mixedsig;
[numSamples,vectorSize] = size(newVectors);
% whitesig = newVectors;
% X        = whitesig; 
X        = newVectors;
% ========================================================
a1              = 1;
a2              = 1;
max_itr         = 1000;
numFailures     = 0;
% How many times do we try for convergence until we give up.
failureLimit    = 2;
econv           = 0.0001;
numIC           = numSamples;
% B               = [];
B = zeros(numSamples);
itr_fail        = 0;
i               = 1;
olddelta=ones(1,numSamples*numSamples);

while i<=numIC
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Show the progress...
    fprintf('IC %d ', i);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    w_old =rand(numIC,1)-0.5;
%     if i>1
%         w_old = w_old-B*B'*w_old; 
%     end
    w_old = w_old-B*B'*w_old; 
    w_old = w_old/norm(w_old);

    for itr=1:max_itr

% gauss
% This has been split for performance reasons.
        u = X' * w_old;
    	u2=u.^2;
    	ex=exp(-a2 * u2/2);
    	gauss =  u.*ex;
    	dGauss = (1 - a2 * u2) .*ex;
    	w = (X * gauss - sum(dGauss)' * w_old) / vectorSize;

%         if i>1
%             w = w-B*B'*w;
%         end

        % Project the vector into the space orthogonal to the space
        % spanned by the earlier found basis vectors. Note that we can do
        % the projection with matrix B, since the zero entries do not
        % contribute to the projection.
        w = w - B * B' * w;
        w = w / norm(w);
      
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if itr == max_itr
            fprintf('\nComponent number %d did not converge in %d iterations.\n', i, max_itr);
            i = i - 1;
            numFailures = numFailures + 1;
        
            if numFailures > failureLimit
                fprintf('Too many failures to converge (%d). Giving up.\n', numFailures);
            end
            if i == 0
                A=[];
                W=[];
            end
            return;
            
            % numFailures > failurelimit
            break;
        end     
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Show the progress...
        fprintf('.');
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
%         if norm(w - w_old) < econv | norm(w + w_old) < econv
        if min( norm (w - w_old) , norm ( w + w_old) )<econv
            itr_fail = 0;
            
            % Save the vector
            B(:, i) = w;	  
            % Calculate the de-whitened vector.
            A(:,i) = dewhiteningMatrix * w;
            % Calculate ICA filter.
            W(i,:) = w' * whiteningMatrix;
%             W(1:2,i) = w' * whiteningMatrix;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Show the progress...
            fprintf('computed ( %d steps ) \n', itr);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            i = i+1;
%             B = [B w];
            
            break
        end

        w_old = w;
                
    end 
    
end        

% [Dim, NumOfSampl] = size(mixedsig);
% M=numIC;N=numIC;
% delta=reshape(B - eye(3),1,M*N);
% change=delta*delta';
% angle=acos((delta*olddelta')/sqrt((delta*delta')*(olddelta*olddelta')));
% fprintf('****change=%.4f angle=%.1f deg., [N=%d,M=%d] \n',...
% change,180*angle/pi,N,M);

% NumOfSampl = numSamples;
% W = w' * whiteningMatrix;
% W = B' * whiteningMatrix;
Y = W * mixedsig + (W * mixedmean) * ones(1, vectorSize);

figure;
subplot(3,1,1);plot(Y(1,:),'r');grid; % plot mixing 1
subplot(3,1,2);plot(Y(2,:),'b');grid; % plot mixing 2
subplot(3,1,3);plot(Y(3,:),'g');grid; % plot mixing 2
suptitle('Unmixed Signals');
% wavplay(Y(1,:),fs1);
% wavplay(Y(2,:),fs2);
% wavplay(Y(3,:),fs3);

temp_y=Y;           % Y = 2 * 48000
min_y(1,1)=min(temp_y(1,:));
min_y(1,2)=min(temp_y(2,:));
min_y(1,3)=min(temp_y(3,:));

max_y(1,1)=max(temp_y(1,:));
max_y(1,2)=max(temp_y(2,:));
max_y(1,3)=max(temp_y(3,:));

abs_max_y(1,1)=max(abs(min_y(1,1)),abs(max_y(1,1)));
abs_max_y(1,2)=max(abs(min_y(1,2)),abs(max_y(1,2)));
abs_max_y(1,3)=max(abs(min_y(1,3)),abs(max_y(1,3)));

scaling_y(1,:)=Y(1,:)./abs_max_y(1,1);
scaling_y(2,:)=Y(2,:)./abs_max_y(1,2);
scaling_y(3,:)=Y(3,:)./abs_max_y(1,3);

figure;
subplot(3,1,1);plot(scaling_y(1,:),'r');grid;%legend('First','Best');
subplot(3,1,2);plot(scaling_y(2,:),'b');grid;%legend('Second','Best');
subplot(3,1,3);plot(scaling_y(3,:),'g');grid;%legend('Third','Best');
suptitle('SCALED Unmixed Signals');

% wavplay(scaling_y(1,:),fs1);
% wavplay(scaling_y(2,:),fs2);
% wavplay(scaling_y(3,:),fs3);

% figure;
% subplot(311);plot(v1,'r.');grid;hold on;plot(scaling_y(1,:),'bo');shg
% subplot(312);plot(s3,'r.');grid;hold on;plot(scaling_y(2,:),'bo');shg
% subplot(313);plot(s1,'r.');grid;hold on;plot(scaling_y(3,:),'bo');shg

% figure;
% subplot(311);plot(v1' - scaling_y(1,:),'r.');grid;
% subplot(312);plot(s3' - scaling_y(2,:),'r.');grid;
% subplot(313);plot(s1' - scaling_y(3,:),'r.');grid;
% suptitle('DIFFERENCE (Original and Unmixed Signals)');
% 
% wavplay((v1' - scaling_y(1,:))*100,fs1);
% wavplay((s3' - scaling_y(2,:)*-1)*100,fs2);
% wavplay((s1' - scaling_y(3,:)*-1)*100,fs3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [icasig] = fastica (X');
% [icasig] = fastica (X', 'approach', 'symm', 'g', 'tanh');
% [icasig] = fastica (X', 'approach', 'defl', 'g', 'tanh');
% figure;
% subplot(3,1,1);plot(icasig(1,:),'r');grid; % plot mixing 1
% subplot(3,1,2);plot(icasig(2,:),'b');grid; % plot mixing 2
% subplot(3,1,3);plot(icasig(3,:),'g');grid; % plot mixing 2
% 
% temp_y=icasig;           % Y = 2 * 48000
% min_y(1,1)=min(temp_y(1,:));
% min_y(1,2)=min(temp_y(2,:));
% min_y(1,3)=min(temp_y(3,:));
% 
% max_y(1,1)=max(temp_y(1,:));
% max_y(1,2)=max(temp_y(2,:));
% max_y(1,3)=max(temp_y(3,:));
% 
% abs_max_y(1,1)=max(abs(min_y(1,1)),abs(max_y(1,1)));
% abs_max_y(1,2)=max(abs(min_y(1,2)),abs(max_y(1,2)));
% abs_max_y(1,3)=max(abs(min_y(1,3)),abs(max_y(1,3)));
% 
% scaling_y(1,:)=icasig(1,:)./abs_max_y(1,1);
% scaling_y(2,:)=icasig(2,:)./abs_max_y(1,2);
% scaling_y(3,:)=icasig(3,:)./abs_max_y(1,3);
% 
% figure;
% subplot(3,1,1);plot(scaling_y(1,:),'r');grid;%legend('First','Best');
% subplot(3,1,2);plot(scaling_y(2,:),'b');grid;%legend('Second','Best');
% subplot(3,1,3);plot(scaling_y(3,:),'g');grid;%legend('Third','Best');
% suptitle('SCALED Unmixed Signals');
% 
% wavplay(scaling_y(1,:),fs3);
% wavplay(scaling_y(2,:),fs1);
% wavplay(scaling_y(3,:),fs2);