ica_ng.m

此代码只要是自然梯度算法（naturalness_gradident_arithmetic

%自然梯度算法     观测信号>=源信号
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