qiaoduu.m

基于ICA的信号分离程序

clc
clear
v=0.5:0.0001:0.6;
sig=[];
sig(1,:)=sin(300*v)+cos(700*v);
sig(2,:)=sawtooth(800*v);                        %1000hz取样
sig(3,:)=square(1000*v);
sig(4,:)=sin(450*v);
for t=1:4
sig(t,:)=sig(t,:)/std(sig(t,:));
end
s1=sig(1,:);
s2=sig(2,:);
s3=sig(3,:);
s4=sig(4,:);
subplot(4,4,1)
plot(v,s1)
xlabel('v')  
ylabel('s1')
subplot(4,4,2)
plot(v,s2)
xlabel('v')
ylabel('s2')
subplot(4,4,3)
plot(v,s3)
xlabel('v')
ylabel('s3')
subplot(4,4,4)
plot(v,s4)
xlabel('v')
ylabel('s4')
Aorig=rand(size(sig,1));
omixedsig=Aorig*sig;
x1=omixedsig(1,:);
x2=omixedsig(2,:);
x3=omixedsig(3,:);
x4=omixedsig(4,:);
figure
subplot(4,4,1);
plot(v,x1)
xlabel('v')
ylabel('x1')
subplot(4,4,2);
plot(v,x2)
xlabel('v')
ylabel('x2')
subplot(4,4,3);
plot(v,x3)
xlabel('v')
ylabel('x3')
subplot(4,4,4);
plot(v,x4)
xlabel('v')
ylabel('x4')
mixedsig=zeros(size(omixedsig));
mixedmean=mean(omixedsig')';
mixedsig=omixedsig-mixedmean*ones(1,size(omixedsig,2));

covMatric=cov(mixedsig',1);
[E,D]=eig(covMatric);
eigenvalue=flipud(sort(diag(D)));
whiteningMatric=inv(sqrt(D))*E';

%用于白化数据的白化矩阵
dewhieningMatric=E*sqrt(D);
whitsig=whiteningMatric*mixedsig;
figure
subplot(4,4,1);
plot(v,whitsig(1,:))
subplot(4,4,2);
plot(v,whitsig(2,:))
subplot(4,4,3);
plot(v,whitsig(3,:))
subplot(4,4,4);
plot(v,whitsig(4,:))
X=whitsig;
[vecterSize,numSamples]=size(X);
A=zeros(vecterSize);
B=zeros(vecterSize);
W=zeros(vecterSize);
D=zeros(vecterSize);
num=vecterSize;
for i=1:num
    w=rand(4,1);   
    w=w/norm(w);
    j=1;   
    while(j<2500) 
        w=w-B*B'*w;     
        w=w/norm(w);                 
        y=w'*X;             
        c=X*(y.^3)'/numSamples;        
        w1=c-3*w;   
        w1=w1/norm(w1);            
        w=w1;
        j=j+1;
    end
    j      
    A(:,i)=dewhieningMatric*w
    W(i,:)=w'*whiteningMatric
    B(:,i)=w    
end
   icasig=W*mixedsig+(W*mixedmean)*ones(1,numSamples);
    y1=icasig(1,:);
    y2=icasig(2,:);
    y3=icasig(3,:);
    y4=icasig(4,:);
    figure
    subplot(4,4,1);
    plot(v,y1)
    xlabel('v')
    ylabel('y1')
    subplot(4,4,2);
    plot(v,y2)
    xlabel('v')
    ylabel('y2')
    subplot(4,4,3);
    plot(v,y3)
    xlabel('v')
    ylabel('y3')
    subplot(4,4,4);
    plot(v,y4)
    xlabel('v')
    ylabel('y4')