function SINR = opast()
%
%2005.12.5
%加載數據
load goldseq;
load bs;
%gold碼歸一化
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;

N=length(s1);
w=s1;
%u=10^-4;


r1=goldseq(1,:)'*b(1,:);
r2=goldseq(2,:)'*b(2,:);
r3=goldseq(3,:)'*b(3,:);
r4=goldseq(4,:)'*b(4,:);
r5=goldseq(5,:)'*b(5,:);
r6=goldseq(6,:)'*b(6,:);

r=r1+A2*r2+A3*r3+A4*r4+A5*r5+A6*r6;
%r=r1+r2+r3+r4+r5+r6;
SNR=20;
r=awgn(r,SNR);

[Us,Ds,Un,Dn] = subspace(r(:,1:50));
k=size(Ds,1);
W=Us;
Z=inv(Ds);
beta=0.98;

%[Us,Ds,Un,Dn] = subspace(r);

for i=1:1000
    
    y=W'*r(:,i);
    q=(1/beta)*Z*y;
    gama=1/(1+y'*q);
    p=gama*(r(:,i)-W*y);
    tao=(1/norm(q))*(1/sqrt(1+norm(p)*norm(q))-1);
    f=tao*W*q+(1+tao*norm(q))*p;
    Z=(1/beta)*Z-gama*q*q';
    W=W+f*q';
    %w=W(:,1);
    D=diag(diag(inv(Z)));
    %w=(1/(s1'*(W*inv(D)*W')*s1))*(W*inv(D)*W')*s1;
    w=W*inv(D)*W'*s1;
    %性能測量
    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;
%    error(i)=norm(W*inv(W'*W)*W'-Us*inv(Us'*Us)*Us');
end
i=1:1000;
plot(i,10*log10(SINR));
xlabel('迭代次數');
ylabel('信號干擾噪聲比SINR(dB)');
grid on;