aquisition.m

matlab仿真的GPS信号捕获程序

clear;
%指标是signalpower=-185dBW，此时可以进行室内定位；
%指标是SNR=-52dB，此时noisepower=-133dBW
Dopler_bound = 10; %原来是20，为加快运行速度，减少频率槽数
Freq_bin = 500;
Fc = 1250000;
Length_of_data = 5000;
Length_of_fft = 4096;
NUM = 5;%50010;
Signalpower = -177 %I路信号的功率（dBW）
Noisepower = -143 %噪声的复功率（dBW）
Threshold = 1e5;
Times_of_ncoh = 40;

%变量初始化
successcount=0;
%在卫星从地平线升起或降落时，用户速度900m/s时,径向速度差不多(+ -)10K，

fid_ca=fopen('CA_5000.txt','r');
CA_samp=fscanf(fid_ca,'%d');%5000点
fclose(fid_ca);
CA_samp=CA_samp';

%本地伪码FFT:CA_fft
CA_local=Average(Length_of_data, Length_of_fft, CA_samp);
%plot(CA_local)
CA_fft=fft(CA_local);
%plot(imag(CA_fft))

Idata=ReadIn( ); %把文本Idata全部读到Idata数组中
Qdata=ReadIn( ); %  
Idata=Idata';
Qdata=Qdata';
Idata(1:NUM)=[]; %删去前NUM个元素
Qdata(1:NUM)=[];

AMP=sqrt(2*10^(Signalpower/10))   %输入量为幅度系数，单位相当于（V），%AMP=1e-9，此时signal power=-183dBW
DIV_AWGN= 10^(Noisepower/10) %输入量为高斯白噪声的方差。DIV_AWGN：在均值means为0时，高斯白噪声的方差即噪声的功率单位为（W）

for k=1:200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %变量初始化
    acq_tag = 0;
    coh_count = 1;
    ncoh_count = 1;
    highdop_index = 1;%初始值
    lowdop_index = 0;
    %初始化矩阵并清零
    R = zeros(20,(2*Dopler_bound+1)*Length_of_fft);
    sum_cohe_accum = zeros(Times_of_ncoh,(2*Dopler_bound+1)*Length_of_fft);%相干累积和清零
    sum_ncohe_accum = zeros(1,(2*Dopler_bound+1)*Length_of_fft);%非相干累积和清零
    ncohe_accum = zeros(1,(2*Dopler_bound+1)*Length_of_fft);
   
    [gr,gi]=MyAwgn(62*5000*20-NUM,DIV_AWGN);%产生随机高斯白噪声：
    Ich = Idata*AMP+gr;
    Qch = Qdata*AMP+gi;
        
    for n=1:Times_of_ncoh %把I和Q转变成二维数组Isn和Qsn
        Isn(n,1:1e5) = Ich(1+(n-1)*1e5 : n*1e5);
        Qsn(n,1:1e5) = Qch(1+(n-1)*1e5 : n*1e5);
    end;
    for ncoh_count=1:Times_of_ncoh%%%进行非相干累积
        while coh_count<=20 %%进行相干累积
            while (highdop_index <= Dopler_bound) | (lowdop_index>=(-1)*Dopler_bound) %一次5000点搜索
                if  lowdop_index >= Dopler_bound*(-1)
                    cos_localcarr = LocalCarrier(lowdop_index, Freq_bin, Fc);
                    sin_localcarr = LocalCarrier(lowdop_index, Freq_bin, Fc);
                    I=Isn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*cos_localcarr-Qsn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*sin_localcarr;
                    Q=Isn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*sin_localcarr+Qsn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*cos_localcarr;
                    I_short = Average( Length_of_data, Length_of_fft, I);
                    Q_short = Average( Length_of_data, Length_of_fft, Q);
                    Rev_fft = fft(I_short+i*Q_short); %接收来的信号做fft
                    R(coh_count,1+(lowdop_index+Dopler_bound)*Length_of_fft : (lowdop_index+Dopler_bound+1)*Length_of_fft) = ifft(Rev_fft.*conj(CA_fft)); %相关值R，R其实是一个矩阵
%                    R_mod(coh_count,1+(lowdop_index+Dopler_bound)*Length_of_fft : (lowdop_index+Dopler_bound+1)*Length_of_fft) = abs(R(coh_count,1+(lowdop_index+Dopler_bound)*Length_of_fft : (lowdop_index+Dopler_bound+1)*Length_of_fft));
                end;
                if  highdop_index <= Dopler_bound
                    cos_localcarr = LocalCarrier(highdop_index, Freq_bin, Fc);
                    sin_localcarr = LocalCarrier(highdop_index, Freq_bin, Fc);
                    I=Isn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*cos_localcarr-Qsn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*sin_localcarr;
                    Q=Isn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*sin_localcarr+Qsn(ncoh_count,1+(coh_count-1)*Length_of_data : coh_count*Length_of_data ).*cos_localcarr;
                    I_short = Average( Length_of_data, Length_of_fft, I);
                    Q_short = Average( Length_of_data, Length_of_fft, Q);
                    Rev_fft = fft(I_short+i*Q_short); %接收来的信号做fft
                    R(coh_count,1+(highdop_index+Dopler_bound)*Length_of_fft : (highdop_index+Dopler_bound+1)*Length_of_fft) = ifft(Rev_fft.*conj(CA_fft));
    %                 R_mod(coh_count,1+(highdop_index+Dopler_bound)*Length_of_fft : (highdop_index+Dopler_bound+1)*Length_of_fft) = abs(R(coh_count,1+(highdop_index+Dopler_bound)*Length_of_fft : (highdop_index+Dopler_bound+1)*Length_of_fft));
                end;
                    lowdop_index = lowdop_index-1;
                    highdop_index = highdop_index+1;
            end; %结束时域搜索频率，求得单次相关值

    %共进行2×Dopler_bound+1次，即 R 实际是 一个21×4096型矩阵
    %再进行相干累积20次，即 R 实际是 一个20×21×4096型矩阵
    %         [peak,locate]=max(R_mod(coh_count, : ));

    %       disp(Threshold);
    %       if peak>Threshold  %门限值要科学慎重计算，调试，得到。
    %          acq_tag=1;
    %          freqbin_index = ceil(locate/Length_of_fft); %向上取整，求出行号即多普勒频槽号
    %          freqbin_index = freqbin_index-Dopler_bound-1;
    %          CA_phase = mod(locate,Length_of_fft);       %取模，求出列号即码相位偏移号
    %          CA_phase = Unaverage(Length_of_fft, Length_of_data, CA_phase);
    %          CA_phase = CA_phase-1;
    %          disp(coh_count);
    %          disp(freqbin_index);
    %          disp(CA_phase);
    %          break;

    %       else
               sum_cohe_accum(ncoh_count, : ) = sum_cohe_accum(ncoh_count, : ) + R(coh_count, : );
               coh_count=coh_count+1;
%               if coh_count == 2 & ncoh_count == 1
%                   for p=1:21
%                    picture(p,:) = R(coh_count, 1+(p-1)*4096 : 4096*p );
%                   end;
%                   mesh(picture);
%               end;%把一次相关结果转化成picture二维矩阵，画三维图
    %       end;
    
        end; %结束20ms的20次相干累积
%         if  ncoh_count == 1
%             for p=1:21
%                 picture(p,:) = sum_cohe_accum(1, 1+(p-1)*4096 : 4096*p );
%             end;
%             mesh(picture);
%         end;%把20次相干累积结果转化成picture二维矩阵，画三维图
            ncohe_accum( : ) = (abs(sum_cohe_accum(ncoh_count, : )/20)).^2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            sum_ncohe_accum( : ) = sum_ncohe_accum( : ) + ncohe_accum( : );%求非相干累积和
%             
    end;%结束40次非相干累积
    for p=1:21
        picture(p, : ) = sum_ncohe_accum( 1+(p-1)*4096 : 4096*p );
    end;
    mesh(picture);
    hold on;
    [peak, locate] = max(sum_ncohe_accum);
    freqbin_index = ceil(locate/Length_of_fft); %向上取整，求出行号即多普勒频槽号
    freqbin_index = freqbin_index-Dopler_bound-1;%freqbin_index取值范围从 -Dopler_bound，到 +Dopler_bound
    CA_phase = mod(locate,Length_of_fft);       %取模，求出列号即码相位偏移号
    CA_phase = Unaverage(Length_of_fft, Length_of_data, CA_phase);
    CA_phase = CA_phase-1;
    disp(peak);
    disp(freqbin_index);
    disp(CA_phase);
    if (freqbin_index==2) & (CA_phase==(Length_of_data-mod(NUM,Length_of_data))|CA_phase==(Length_of_data-mod(NUM,Length_of_data)-1)|CA_phase==(Length_of_data-mod(NUM,Length_of_data)-2)|CA_phase==(Length_of_data-mod(NUM,Length_of_data)+1))
        successcount = successcount+1;
        acq_tag = 1;
    end;
end;
    
    disp(successcount);

%     fclose(fid1);%捕获结束，关掉文本指针
%     fclose(fid2);%捕获结束，关掉文本指针