pastd.m

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

function SINR=pastd()
%
%2005.12.6
load bs;
load goldseq;
goldseq = 1/sqrt(31)*goldseq;
%b(6,1:1000)=zeros(1,1000);
s1=goldseq(1,:)';
s2=goldseq(2,:)';
s3=goldseq(3,:)';
s4=goldseq(4,:)';
s5=goldseq(5,:)';
s6=goldseq(6,:)';

r1=s1*b(1,:);
r2=s2*b(2,:);
r3=s3*b(3,:);
r4=s4*b(4,:);
r5=s5*b(5,:);
r6=s6*b(6,:);
A1=1;
A2=10;
A3=10;
A4=10;
A5=0;
A6=0;

r=A1*r1+A2*r2+A3*r3+A4*r4+A5*r5+A6*r6;
SNR=20;
r=awgn(r,SNR);
%[dim,N]=size(r);
N=31;
u=10^-4;
mixedsig=r(:,1:50);
[Us,Ds,Un,Dn] = subspace(mixedsig);
K=size(Ds,1);
D=[diag(Ds);diag(Dn)];
U=[Us Un];
beta=0.97;

sigma=mean(diag(Dn));
for t=51:2050
    x=r(:,t);
    for i=1:K
        y=U(:,i)'*x;    
        D(i)=beta*D(i)+u*abs(y)^2;
        e=x-U(:,i)*y;
        U(:,i)=U(:,i)+u*e*(y/D(i));
        %Us(:,i)=U(:,i)+(x-U(:,i)*y)*y/D(i);
        x=x-U(:,i)*y;
    end   
    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(t-50)=((w'*s1)^2)/n;
end

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