showicawav.m

这是个关于ICA的MATLAB的源程序

%ShowICAwav
% Mixture two sound and then separat them
close all
clear all

% Read the sound data from wave file,x-data,f-sampling rate, b-data precision
[x1,f1,b1]=wavread('so1.wav');
[x2,f2,b2]=wavread('so3.wav');
x1=x1/std(x1);
x2=x2/std(x2);
    % play the sound
    input('Anykey to play sound 1 :');
    wavplay(x1,f1);
    input('Anykey to play sound 2 :');
    wavplay(x2,f2);
X=[x1';x2'];
N=length(x1);

clear x1 x2

% Cancel mean, original data
        figure(1)
        subplot(4,1,1)
            plot(X(1,:),'r','linewidth',3)
        subplot(4,1,2)
            plot(X(2,:),'b-.','linewidth',3)
        subplot(2,3,4)
            [H,A]=hist(X(1,:));
            bar(A,H,'r')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
        subplot(2,3,5)
            [H,A]=hist(X(2,:));
            bar(A,H,'b')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
        subplot(2,3,6)
            plot(X(1,:),X(2,:),'*','linewidth',3)
            axis([min(X(1,:))-0.2,max(X(1,:))+0.2,min(X(2,:))-0.2,max(X(2,:))+0.2])
            
% Mixture data
A=[0.7,0.3;0.4,0.6];
X=A*X;
        figure(2)
        subplot(4,1,1)
            plot(X(1,:),'r','linewidth',3)
        subplot(4,1,2)
            plot(X(2,:),'b-.','linewidth',3)
        subplot(2,3,4)
            [H,A]=hist(X(1,:));
            bar(A,H,'r')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
        subplot(2,3,5)
            [H,A]=hist(X(2,:));
            bar(A,H,'b')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10]);
        subplot(2,3,6)
            plot(X(1,:),X(2,:),'*','linewidth',3)
            axis([min(X(1,:))-0.1,max(X(1,:))+0.1,min(X(2,:))-0.1,max(X(2,:))+0.1])
                % play the sound
                input('Anykey to play sound 1 :');
                wavplay(X(1,:),f1);
                input('Anykey to play sound 2 :');
                wavplay(X(2,:),f2);

% Cancel and Whitening

    X(1,:)=X(1,:)-mean(X(1,:));
    X(2,:)=X(2,:)-mean(X(2,:));
    CXX=X*X'/N;
    [D,L]=eig(CXX);
    X=sqrt(inv(L))*D*X;
        figure(3)
        subplot(4,1,1)
            plot(X(1,:),'r','linewidth',3)
        subplot(4,1,2)
            plot(X(2,:),'b-.','linewidth',3)
        subplot(2,3,4)
            [H,A]=hist(X(1,:));
            bar(A,H,'r')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
        subplot(2,3,5)
            [H,A]=hist(X(2,:));
            bar(A,H,'b')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
        subplot(2,3,6)
            plot(X(1,:),X(2,:),'*','linewidth',3)
            axis([min(X(1,:))-0.2,max(X(1,:))+0.2,min(X(2,:))-0.2,max(X(2,:))+0.2])
% FastIAC
    w1=rand(1,2);
    w1=w1/sqrt(sum(w1.^2));
    for k=1:10,
        k
        y1=w1*X;
        y1=y1-mean(y1);
        k4=sum(y1.^4)/length(y1);
        k2=sum(y1.^2)/length(y1);
        kt1(k)=k4/k2-3;
        dw=[0;0];
        for n=1:N
            dw=dw+X(:,n)*(w1*X(:,n))^3;
        end
        w1=sum(dw')/length(dw)-3*w1;
        w1=w1/sqrt(sum(w1.^2));
    end
    
    w2=rand(1,2);
    w2=w2/sqrt(sum(w2.^2));
    for k=1:50
        w2=w2-(w2*w1')*w1;
    end
    y2=w2*X;
    
    y1=(y1-mean(y1))/sqrt(sum(y1.^2)/length(y1));
    y2=(y2-mean(y2))/sqrt(sum(y2.^2)/length(y2));
        % play the sound
    input('Anykey to play sound 1 :');
    wavplay(y1,f1);
    input('Anykey to play sound 2 :');
    wavplay(y2,f2);

    
        figure(4)
        subplot(4,1,1)
            plot(y1,'r','linewidth',3)
        subplot(4,1,2)
            plot(y2,'b-.','linewidth',3)
        subplot(2,3,4)
            [H,A]=hist(y1);
            bar(A,H,'r')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
         subplot(2,3,5)
            [H,A]=hist(y2);
            bar(A,H,'b')
            dA=(A(2)-A(1))/1.9;
            axis([min(A)-dA,max(A)+dA,0,max(H)+10])
         subplot(2,3,6)
            plot(y1,y2,'*','linewidth',3)
            axis([min(y1)-0.2,max(y1)+0.2,min(y2)-0.2,max(y2)+0.2])
            
        figure(5)
        hold on
            plot(kt1,'r','linewidth',3)
            
       W=[w1;w2]