%自然梯度算法     觀測信號>=源信號
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