figure8success.m

老师给的一个范本

zclc;clear all;
a=1/2;
M=4;
n=100;        %
pi=3.14159;      
id=sqrt(-1);
T=1/32000;
sample_rate=8;
fs=sample_rate/T;
n_t=[-3 3];
tao=0;
theta=0;
%-------------------------------------
%码元产生fax   ran
times=100;     % 100 times
v=0:320:4800;% 
Tv=v*T;
snr=2;
sigma=sqrt(1/(2*10.^(snr/10)));      % sqrt(N0/2)
ve1=zeros(times,length(v));  
ve2=zeros(times,length(v));  
ve3=zeros(times,length(v));  
Tavve1=zeros(1,length(v));
Tavve2=zeros(1,length(v));
Tavve3=zeros(1,length(v));
% var1=zeros(times,length(v));
% var2=zeros(times,length(v));
% var3=zeros(times,length(v));
% Tavvar1=zeros(1,length(v));
% Tavvar2=zeros(1,length(v));
% Tavvar3=zeros(1,length(v));
% ==============
for t=2:length(v)   % v(0)=0
  for k=1:times
  ran=randn(1,n);                  %随机产生Bpsk的信号比特1或-1
   c=zeros(1,n);
   for i=1:n
      if ran(1,i)>0
         c(1,i)=1;
      else
        c(1,i)=-1;
      end   
   end;
%  calculate ck\
ck=zeros(1,length(c));
for i=1:(length(c)-1)     %  
     switch c(i)
      case -1,
         switch c(i+1)
             case -1,%[0, 0] case
               phase=pi/4;
           case 1,%[0,1] case
                phase=3*pi/4;
         end  
      case 1,
        switch c(i+1)
           case -1,%[1,0] case
               phase=7*pi/4;
           case 1,%[1,1] case
               phase=5*pi/4;
        end
     end
   ck(i)=exp(id*phase);      %(3.2.14)
end
% interper into zeros 
ck_interp=zeros(1,length(ck)*sample_rate);
for i=1:length(ck)
      ck_interp((i-1)*sample_rate+1)=ck(i);
end     
% rcosfir
[g sample_time]=rcosfir(a,n_t,sample_rate,T);   % the same concequence
cg=conv(ck_interp,g);
r=zeros(1,length(cg));
y1=zeros(1,length(cg));
for i=1:length(cg)
    r(i)=cg(i)*exp(id*(2*pi*v(t)*(i/fs+tao)+theta));  %
end
% add noise
y1=r+sigma*(randn(1,length(r))+id*randn(1,length(r)));
% choose point from y1
b=25;
y=zeros(1,n);  
for i=1:length(y)
      y(i)=y1(sample_rate*(i-1)+b);
end
% calculate z
z=zeros(1,length(y));
for i=1:length(y)
    z(i)=y(i)*conj(ck(i));
end
%  may method
L0=32;
gama=zeros(1,(L0-1));
for i=1:L0-1
     gama(i)=3/2*L0/(L0^2-1)*(1-((2*i-L0)/L0)^2);
end
for i=2:L0
    ve1((t-1)*times+k)=ve1((t-1)*times+k)+gama(i-1)*angle(z(i)*conj(z(i-1)));
 end
ve1((t-1)*times+k)=ve1((t-1)*times+k)/(2*pi*T)*T;%
% var1((t-1)*times+k)=((ve1((t-1)*times+k)-v(t))/v(t))^2;
Tavve1(t)=Tavve1(t)+ve1((t-1)*times+k)/times;
% Tavvar1(t)=Tavvar1(t)+var1((t-1)*times+k)/times;
%  fitz method
N=16;
R=zeros(1,L0-1);
for m=1:L0-1
    for  i=m+1:L0        %   equal to k=(m+1):L0
      R(m)=R(m)+1/(L0-m)*z(i)*conj(z(i-m));     %???(3.2.30)
    end
end
for m=1:N
   ve2((t-1)*times+k)=ve2((t-1)*times+k)+angle(R(m));
end
ve2((t-1)*times+k)=ve2((t-1)*times+k)/(pi*N*(N+1)*T)*T;%
% var2((t-1)*times+k)=((ve2((t-1)*times+k)-v(t))/v(t))^2;
Tavve2(t)=Tavve2(t)+ve2((t-1)*times+k)/times;
% Tavvar2(t)=Tavvar2(t)+var2((t-1)*times+k)/times;
% L&R method
for m=1:N
    ve3((t-1)*times+k)=ve3((t-1)*times+k)+R(m);
end
ve3((t-1)*times+k)=angle(ve3((t-1)*times+k))/(pi*(N+1)*T)*T;%
% var3((t-1)*times+k)=((ve3((t-1)*times+k)-v(t))/v(t))^2;
Tavve3(t)=Tavve3(t)+ve3((t-1)*times+k)/times;
% Tavvar3(t)=Tavvar3(t)+var3((t-1)*times+k)/times;
end
end
%  drawn pictures
figure(1)
plot(Tv,Tavve1,'-');
hold on
plot(Tv,Tavve2,'-');
hold on
plot(Tv,Tavve3,'-');
hold off
xlabel('Normalized frequency, Tv'),ylabel('Average Estimates E(Tv)');
%  =============
Tv1=-1*fliplr(Tv);
Tv2=[Tv1 Tv(2:end)];
Tavve11=-1*fliplr(Tavve1);
Tavve12=[Tavve11 Tavve1(2:end)];
Tavve21=-1*fliplr(Tavve2);
Tavve22=[Tavve21 Tavve2(2:end)];
Tavve31=-1*fliplr(Tavve3);
Tavve32=[Tavve31 Tavve3(2:end)];
figure(2)
plot(Tv2,Tavve12);
hold on
plot(Tv2,Tavve22);
hold on
plot(Tv2,Tavve32);
hold off
xlabel('        Normalized frequency, Tv'),ylabel('         Average Estimates E(Tv)');
gtext('Kay');
hold on
gtext('L&R');
hold on
gtext('Fitz');
hold off