www.pudn.com > fastfixedpoint.rar > fastfixedpoint.m

%经典快速固定点算法（基于峭度，梯度下降，超高斯） 
 
clear all; 
 
% 初始化独立源数M和数据长度N 
M=3;%4; 
sn=3;       %周期数 
sl=128;     %每个周期128点 
N=sl*sn;     
 
 
load EPsignal1r; 
% 产生独立信源 
signal1=zeros(1,N);         %清零 
for i=0:sn-1 
    signal1(i*sl+1:(i+1)*sl)=EPsignal1;            %将信号周期延拓，模拟EP信号 
end 
signal2=((rand(1,N)<0.5)*2-1).*log(rand(1,N));     %模拟超高斯信号  
signal3=randn(1,N);         %产生白高斯信号 
%signal4=rand(1,N);         %均匀分布模拟亚高斯信号 
%源的均值要求为0，去直流分量 
signal1=signal1-mean(signal1);        %去均值 
signal1=signal1/max(abs(signal1));    %归一化 
signal2=signal2-mean(signal2); 
signal3=signal3-mean(signal3); 
%signal4=signal4-mean(signal4); 
signal=[signal1;signal2;signal3];%;signal4]; 
 
 
% 初始化 
v=zeros(M,N);    % 混合后阵矩 v=As 
x=zeros(M,N);    % 白化后阵矩 x=Mv=MAs 
y=zeros(M,N);    % 输出信号 
w=zeros(M,1);    % 解混向量，列 
                 % 初值均为零 
A=rand(M);       % 混合矩阵，随机产生 
%wold=rand(M,1);    % 存放k-1时刻的解混向量，用随机数初始化权向量，放在迭代部分也可以 
Wopt=zeros(M,M);    % 权值收敛后的最优权值矩阵 
B=zeros(M);         % 由已经求得的解混向量构成的矩阵。为了每次提取出不同的独立源，需要对解混向量进行正交化。 
flag=zeros(1,M);       % 判断收敛的标志。收敛，flag=1；不收敛，flag=0。行向量 
stepsize=N*ones(1,M);  % 收敛步长初始化，每个源的初始步长均为最大值即总点数 
 
 
% 信号混合 
v=A*signal; 
 
 
% 白化  
Ev=zeros(M); 
eigvalue=zeros(M);    %特征值 
eigvector=zeros(M);   %每一列对应一个特征向量 
WhiteM=zeros(M);      %白化矩阵 
for i=1:N 
    Ev=Ev+v(:,i)*v(:,i)'; 
end 
Ev=Ev/N ;             %近似协防差矩阵 E[xx']，平均近似期望 
[eigvector,eigvalue]=eig(Ev); 
for j=1:M 
    eigvalue(j,j)=power(eigvalue(j,j),-0.5);%此时eigvalue为对角阵,对角元为特征值开方 
end 
WhiteM=eigvector*eigvalue*eigvector'; 
x=WhiteM*v; 
 
 
% 权值迭代，产生输出 
meanvalue=zeros(M,1);          % 迭代算法中要用到,用来纪录,初始为0 
for si=1:M                     %用sn不好==si 表示源 
    wold=rand(M,1);            %列，我觉得把他放到循环之内比较好? 
    if  si>1 
        wold=wold-B*B'*wold;   %投影运算，为满足正交化条件,每个分量之间保持独立 
    end 
    for i=1:N 
        if flag(1,si)==0 
            y(si,i)=wold'*x(:,i); 
            for j=1:N 
                meanvalue=meanvalue+x(:,j)*power(wold'*x(:,j),3); 
            end 
            meanvalue=meanvalue/N; 
            w=meanvalue-3*wold;       %迭代公式实现 
            meanvalue=zeros(M,1);     %此步清零很重要            
            if si>1 
                w=w-B*B'*w; 
            end 
            w=w/norm(w);              %归一化,除以范数 
            if abs(abs(w'*wold)-1)<1e-6 
                flag(1,si)=1; 
                stepsize(1,si)=i;     %步长是根据收敛确定的 
                Wopt(:,si)=w; 
                B(:,si)=w;            %w是矩阵相应的一列，B矩阵用来投影运算 
            else 
                wold=w;  
            end 
        else 
            y(si,i)=Wopt(:,si)'*x(:,i); 
        end 
    end 
end  
 
 
% 画图 
figure(1);                   %源信号 
subplot(M,3,1); 
plot(signal1); 
subplot(M,3,4); 
plot(signal2); 
subplot(M,3,7); 
plot(signal3); 
%subplot(M,3,10); 
%plot(signal4); 
 
subplot(M,3,2);              %混合信号 
plot(v(1,:)); 
subplot(M,3,5); 
plot(v(2,:)); 
subplot(M,3,8); 
plot(v(3,:)); 
%subplot(M,3,11); 
%plot(v(4,:)); 
 
subplot(M,3,3);             %输出信号 
%plot(x(1,:)); 
plot(y(1,:)); 
subplot(M,3,6); 
%plot(x(2,:)); 
plot(y(2,:)); 
subplot(M,3,9); 
%plot(x(3,:)); 
plot(y(3,:)); 
%subplot(M,3,12); 
%plot(x(4,:)); 
%plot(y(4,:)); 
 
 
 
%stepsize=[5	 8	2] 
 
 
%A=0.11992	0.84077	0.71796 
%  0.32012	0.70461	0.13022 
%  0.96523	0.14933	0.68027 
 
 
%Wopt=0.62105	-0.48073	0.61903 
%     0.77039	0.51973	   -0.36929 
%     -0.1442	0.70625	    0.69313