pastd.asv

这是盲信号的代码 都已经通过编译了 做这方面的同仁可以参考一下 我觉得蛮惯用的

function SINR=pastd()
%
%2005.12.6
load r;
[dim,N]=size(r);
mixedsig=r(:,1:50);
[Us,Ds,Un,Dn] = subspace(mixedsig);
K=size(Ds,1);
D=[diag(Ds);diag(Dn)];
U=[Us Un];
%D=flipud(D);
beta=0.99;
sigma=mean(diag(Dn));
for t=1:1000
    x=r(:,t);
    for i=1:K
        y=U(:,i)'*x;    
        D(i)=beta*D(i)+abs(y)^2;
        Us(:,i)=U(:,i)+(x-U(:,i)*y)*y/D(i);
        x=x-U(:,i)*y;
    end
    sigma=beta*sigma+norm(x)/(N-K);
    for k=1:K
        a=0;
        b=1;
        for i=k+1:N
            a=a+D(i)/(N-k);
            b=b*D(i);
        end
        afa(k)=a/(b^(1/(N-k)));
        AIC(k)=(N-k)*ln(afa(k))/(1-beta)+k*(2*N-k);    
    end
    [m,Kt]=min(AIC);
    K=Kt+1;

    
    Us=U(:,1:K);
    Ds=diag(D(1:K));
    afa=s1'*Us*(Ds^-1)*Us'*s1;
    w=(Us*(Ds^-1)*Us'*s1)/afa;
    
    n=(A2*w'*s2)^2+(A3*w'*s3)^2+(A4*w'*s4)^2+(A5*w'*s5)^2+(A6*w'*s6)^2+w'*0.01*w;
    SINR(i)=((w'*s1)^2)/n;
end

i=1:1000;
%plot(i,10*log10(SINR));
%%semilogy(i,SINR);
%xlabel('迭代次数');
%ylabel('信号干扰噪声比SINR(dB)');