www.pudn.com > ICAbss.rar > ica_ng.m

%自然梯度算法     观测信号>=源信号 
clear 
k=4000;                                                   %数据点                
fs=10000;                                                 %采样频率 
for t=1:k 
    s(1,t)=sign(cos(2*pi*155*t/fs));                      %符号信号 
    s(2,t)=sin(2*pi*800*t/fs);                            %高频正弦信号 
    s(3,t)=sin(2*pi*90*t/fs);                             %低频正弦信号 
    s(4,t)=sin(2*pi*9*t/fs)*sin(2*pi*300*t/fs);           %幅值调制信号 
end 
a=1;                                                      %生成[-a,a]的均匀分布随机噪声 
s(5,:)=a-2*a*rand(1,k);                                    
 
figure(1)                                                 %源信号图 
for n=1:5 
subplot(5,1,n);plot(s(n,:)); 
end 
 
A=a-2*a*rand(5,5);                                        %混合矩阵,[-a,a]的均匀分布 
x=A*s;                                                    %观测信号 
 
figure(2)                                                 %观察信号图 
for n=1:5 
    subplot(5,1,n);plot(x(n,:)); 
end 
 
for n=1:5                                                %观测信号零均值处理 
    ave=mean(x(5,:)); 
    x(5,:)=x(5,:)-ave; 
end 
 
I=eye(5,5);                                               %生成单位矩阵 
w1=0.5*eye(5,5);                                       %初始化W1 
for t=1:k 
    y(:,t)=w1*x(:,t); %迭代  
    for l=1:5  
         FIy(l)=y(l,t)^3;   %非线性函数 
    end 
    w1=w1+0.005*(I-FIy'*y(:,t)')*w1; 
     
    c=w1*A;                                               %性能矩阵 
    for p=1:5                                             %计算串音误差 
        max1(p)=abs(c(p,1));                            
        for q=1:5                                 
            if max1(p)<=abs(c(p,q)) 
                max1(p)=abs(c(p,q)); 
            else max1(p)=max1(p); 
            end 
        end 
    end 
    s2=0;    
    for p=1:5 
        s1=0; 
        for q=1:5 
            s1=s1+abs(c(p,q))/max1(p); 
        end 
        s2=s2+abs(s1-1); 
    end 
     
     
    for q=1:5 
        max2(q)=abs(c(1,q)); 
        for p=1:5 
            if max2(q)<=abs(c(p,q)) 
                max2(q)=abs(c(p,q)); 
            else max2(q)=max2(q); 
            end 
        end 
    end 
     
    s4=0; 
    for q=1:5 
        s3=0; 
        for p=1:5 
            s3=s3+abs(c(p,q))/max2(q); 
        end 
        s4=s4+abs(s3-1); 
    end 
    e(t)=s2+s4;  
  
  
end 
 
             
     
y=w1*x;                                               %估计源信号 
 
figure(3)                                             %估计源信号图 
for n=1:5 
    subplot(5,1,n);plot(y(n,:)); 
end 
 
figure(4)                                             %串音误差图 
plot(e) 
 
%FIy=3*(y(:,t).^11)/4+15*(y(:,t).^9)/4-14*(y(:,t).^7)/3-19*(y(:,t).^5)/4+29*(y(:,t).^3)/4