blinding_timing.m

来自「超宽带（UWB）系统的盲同步仿真程序」· M 代码 · 共 136 行

%%%%%%%%%%  The simulation of blind UWB timing with dirty template  %%%%%%%%%%%  
%%%%%%%%%%%%%%%%%%%%%  start  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

%%%%%%%%%%%%%%%%%%%%%%  Parameters Setting  %%%%%%%%%%%%%%%%%%%%%%%%%%%

Tp=1;                                    %  Duration of a monocycle
Tf=100;                                  %  Duration of a frame
Nf=32;                                   %  The number of impulse per symbol
tm=0.288;                                %  The selection of tm here is to ensure the waveform is FCC-comply                                      
M=8;                                    %  The size of the observation interval
Num=20;                                  %  The number of data group for MSE
Sam_Rate=10;                             %  The number of sampling during a monocycle period---Tp
Tfsam=Tf*Sam_Rate;                       %  The number of sampling during a frame---Tf
Tssam=Tf*Sam_Rate*Nf;                    %  The number of sampling during a symbol---Ts
Nmultipath=10;                           %  The number of multipath
E=14073;                                 %  Determined by t3
%%%%%%%%%%%  Generates the waveform of the monocycle P(t)  %%%%%%%%%%%%%

k=floor(Sam_Rate/2);
e=mod(Sam_Rate,2);
t=linspace(0,Tp/2,k);    
s=(1-4.*pi.*((t./tm).^2)).*exp(-2.*pi.*((t./tm).^2));
if e
    for i=1:k
        Pt(k+1+i)=s(i);
        Pt(k+1)=1;
        Pt(k+1-i)=s(i);
    end
else
    for i=1:k
        Pt(k+i)=s(i);
        Pt(k+1-i)=s(i);
    end
end
en=Pt*Pt';                        
A=(1/en)^0.5;
Pt=A.*Pt;                    %  Energy normalizition
clear k e t s A i en 
%plot(Pt)

%%%%%%%%%%%%%  Genetates the transmitted signal  %%%%%%%%%%%%%%%%%%%%%

Timing_data=zeros(Num,M+2);
Timing_data=randint(Num,M+2);
seq=zeros(1,Tssam*(M+2));
PTt1=zeros(1,Tssam*(M+2)+Sam_Rate-1);
PTt=zeros(Num,Tssam*(M+2));   
for k=1:Num
    for i=1:M+2
        for j=1:Nf
            index=(i-1)*Tssam+(j-1)*Tfsam+1;        
            seq(index)=1-2*Timing_data(k,i);
        end
    end
    PTt1=conv(seq,Pt);
    PTt(k,:)=PTt1(1:Tssam*(M+2));
end
clear Pt seq PTt1 index k i j Timing_data

%%%%%%%%%%%%%%%%%%%%%%%  Multipath Channel  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Reference to J. R. Foerster, "The effects of multipath interference on 
%%%%%%%%%  the performance of UWB systems in an indoor wireless channel"  

%  Generates the attenuation of all the users, following Rayleigh distribution
atten=zeros(1,Nmultipath);
delta=-log(0.001)/Nmultipath;    
nc=0;
for L=0:Nmultipath-1
    nc=nc+exp(-delta*L);
end
oumiga0=1/nc;    
for L=0:Nmultipath-1
    atten(L+1)=raylrnd(sqrt(oumiga0*exp(-delta*L)/2));
end
atten=sort(atten,'descend');     %  Sorting is to approve the performance of Rake receiver

%  introduce propagation attenuation and multipath delay

Rect=zeros(Num,Tssam*(M+2)+(Nmultipath-1)*Sam_Rate);
for k=1:Num
    for L=0:Nmultipath-1
        Rect(k,:)=Rect(k,:)+atten(L+1)*[zeros(1,L*Sam_Rate),PTt(k,:),zeros(1,(Nmultipath-1-L)*Sam_Rate)];
    end
    PTt(k,:)=Rect(k,1:Tssam*(M+2));
   
end
clear Rect atten delta L nc oumiga en A 
                   

%%%%%%%%%%%%%%%%%%%%%%%%%%%  acquisition of timing  %%%%%%%%%%%%%%%%%%%%%%%%
            

SNR=0:2:20;
for ind=1:length(SNR)    
    Ne=zeros(1,Num);    
    for k=1:Num
        Rt=zeros(1,Tssam*(M+2));
        Rt=awgn(PTt(k,:),SNR(ind));              %  The received signal
        %subplot(2,1,1)
        %plot(PTt(1,:))
        %subplot(2,1,2)
        %plot(Rt)
        Xt(1:Tssam*(M+2)-E)=Rt(E+1:Tssam*(M+2));  
        clear Rt
        Rxx=zeros(1,Nf);
        for n=0:Nf-1
            for m=0:M/2-1
                Rx=abs(Xt((n*Tfsam+2*m*Tssam+1):(n*Tfsam+(2*m+1)*Tssam))*...
                    Xt((n*Tfsam+(2*m+1)*Tssam+1):(n*Tfsam+(2*m+2)*Tssam))');
                Rxx(n+1)=Rxx(n+1)+Rx;
            end
            Rxx(n+1)=Rxx(n+1)/M*2;
        end
        mRxx=max(Rxx);                        %  The maximum value of Rxx
        for n=0:Nf-1
            if Rxx(n+1)==mRxx
                Ne(k)=n;
            end
        end
        clear mRxx Rxx Xt n m Rx
    end
    MSE(ind)=(sum((Ne-17073).^2))/Num/Nf^2  ;   %  The MSE of Ne
    clear Ne
end