aic.m

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

function AIC()
%
%
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=10;
A6=10;

r=A1*r1+A2*r2+A3*r3+A4*r4+A5*r5+A6*r6;
SNR=20;
r=awgn(r,SNR);
N=31;
mixedsig=r(:,1:50);
[Us,Ds,Un,Dn] = subspace(mixedsig);
K=size(Ds,1);
D=[diag(Ds);diag(Dn)];
U=[Un Us];
beta=0.95;
for t=1:2000;
    c=cov(r(:,t:t+50-1)',1);
    [E,D]=eig(c);
    D=diag(D);
    D=flipud(D);
%    K=sum(D>0.5);
    for k=1:K
        a=0;
        b=1;
        for i=k+1:N
            a=a+D(i);
            b=b*D(i);
        end 
        temp=N-k;
         
            
        a=a/temp;
        b=b^(1/temp);
        
        afa(k)=a/b;
        AIC(k)=temp*log(afa(k))/(1-beta)+0.5*k*(2*N-k)*log(1-beta);    
    end
    [m,Kt]=min(AIC);
    K=Kt+1;
    dim(t)=Kt;
end
plot(dim);
figure(2)
plot(D);