这是子空间算法中的pastd算法，但结果不对，高手们帮忙啊，本人不胜感激
function SINR = pastd()
%加载数据
fbconnection=[0 1 0 0 1]; % gold序列
mseq=m_sequence(fbconnection);
fbconnection1=[0 0 1 0 1];
fbconnection2=[0 1 1 1 1];
goldseq=gold_seq(fbconnection1,fbconnection2);
N=2^length(fbconnection)-1;
ind1=find(mseq==0);
mseq(ind1)=-1;
ind2=find(goldseq==0);
goldseq(ind2)=-1;
bit_number=10000; % 用户数据
K=31;
r=zeros(K,bit_number);
for i=1:K
r(i,:)=randsrc(1,bit_number);
end
goldseq = 1/sqrt(31)*goldseq;%扩频码归一化
s1 = goldseq(1,:)';
s2 = goldseq(2,:)';
s3 = goldseq(3,:)';
s4 = goldseq(4,:)';
s5 = goldseq(5,:)';
s6 = goldseq(6,:)';
A2 = 10;
A3 = 10;
A4 = 10;
A5 = 10;
A6 = 10;
%[dim, N] = size(r)
N = 31;
mixedsig = r(:,1:50);
[Us,Ds,Un,Dn] = subspace(mixedsig);
K = size(Ds,1);
D = [diag(Ds);diag(Dn)];
U = [Us Un];
beta = 0.9;
sigma = mean(diag(Dn));
for t = 51:1050
x = r(:,t);
for i = 1:K
y = U(:,i)'*x;%y=u'*x
D(i) = beta*D(i)+abs(y)^2;
e = x-U(:,i)*y;
Us(:,i) = U(:,i)+e*(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);
b = b*D(i);
end
a = a/(N-k);
afa(k) = a/(b^(1/(N-k)));
AIC(k) = (N-k)*log(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(t-50) = ((w'*s1)^2)/n;
end
i = 1:1000;
plot(i,10*log10(SINR));
xlable('迭代次数');
ylable('信干比SINR(dB)');
grid on