acquisition_nc_gps_sim.m

GPS acquisition 使用noncoherent 方式积分

%  function [] = acquisition( )
% The function perform an acquisition process for  satellites
% yes = 1 if the satellite signal is acquired
% yes = 0 if the satellite signal is NOT acquired
% f = carrier frequency
% ph = code phase (samples)
% K = number of sequential correlations to be accumlated
% Threshold = threshold (typically 6 if K = 10)

clc
clear all;
% close all;
%integration time
n=20;
%satellite prn vector
sat_id_vector =[2 4 6 8 10 12 15 24 26 27 29 129 137 127];
% sat_id_vector =[127 128];
%sat_id_vector =[27];
% carrier center frequency
fc =4.1304e6;
%doppler frequency shift range
fc_lo=fc-7e3;
fc_hi=fc+7e3;
%doppler frequency shift step
step = 100;
%sampling frequency
fs = 16.368e6;
%acquistion threshold
Threshold = 10;
%chip sample
csample = fs/1.023e6*100;

[fname.name , fname.pathname]=uigetfile({'*.sim','All Files(*.*)'},'LOAD left rawif');
length_1ms = fs/1e3;

n_num=size(n,2);
data=zeros(length_1ms,1);
freq = fc_lo:step:fc_hi;

nn = 0:length_1ms*n-1;

plot_col = ceil(size(sat_id_vector,2)/3);
if  fname.name ~= 0

    %Get number of satellite
    sat_num = size(sat_id_vector,2);
    %Open Raw_IF file
    fid = fopen([fname.pathname fname.name],'rb');
    IF_data  = (fread(fid, fs/1e3*n, 'int8'))*2+1;

    
    % time 
    figure(31);
    time = IF_data(1:length_1ms,1);
    hist(IF_data);
   
    % Plot sprectrum
    figure(12);
    spectrum = abs(fft(IF_data(1:length_1ms,1))).^2;
    plot(spectrum);
    figure(13);
    sp = 10*log10(spectrum(2:length_1ms/2)/mean(spectrum(2:10))) - 144;
    plot(1000:1000:fs/2-1000,sp);
    xlabel('Frequency (MHz)');
    ylabel('Power Sprectal Density (dBW/Hz)');
    title('IF Signal Spectrum')
    axis([1000 fs/2 -inf inf])

    % The noise power per second
    P_n = mean(spectrum(2:100))/1000;

    for sat_i = 1:sat_num
        %Select satellite PRN
        sat_id = sat_id_vector(sat_i);

        C = zeros(length(fc_lo:step:fc_hi),length_1ms);
        Correlator_I = zeros(length(fc_lo:step:fc_hi),length_1ms);
        Correlator_Q = zeros(length(fc_lo:step:fc_hi),length_1ms);
        bar1 = waitbar(0,'RUN -- Acquisition');
        counterl=size(freq,2)*2;

        tic
        % K sequential correlation is accumulated
        %         % Generate C/A code
        %         CA = cacode(sat_id, 1.023e6, fs, length_1ms);
        %         % DFT of C/A code
        %         F_CA = conj(fft(CA));

        %         for nonch_i =1:n(n_num)
        %             data = IF_data(((nonch_i-1)*length_1ms + 1):(nonch_i*length_1ms),1);
        %             i=0;
        %             for fcc = fc_lo:step:fc_hi
        %                 i = i+1;
        %                 % Correlation
        % %                 I = cos(2*pi*fcc/fs*nn).*data';
        % %                 Q = sin(2*pi*fcc/fs*nn).*data';
        % %                 Xk = fft(I+ j*Q);
        %                 Xk = fft(exp(j*2*pi*fcc/fs*nn).*data');
        %                 corcodeIQ = Xk.* F_CA;
        %                 Cor_Result = abs(ifft(corcodeIQ));
        %                 C(i,:) = C(i,:) + Cor_Result(1,1:length_1ms);
        %                 waitbar((i+(nonch_i-1)*size(freq,2))/(size(freq,2)*n));
        %             end
        %
        %         end

        % Generate C/A code
        CA = cacode(sat_id, 1.023e6, fs, length_1ms);
        % DFT of C/A code
        F_CA = conj(fft(CA));

        i=0;
        for fcc = fc_lo:step:fc_hi
            i = i+1;
            mixed_signal =exp(j*2*pi*fcc/fs*nn).*IF_data';
            for nonch_i =1:n(n_num)
                Xk = fft(mixed_signal(1,((nonch_i-1)*length_1ms + 1):(nonch_i*length_1ms)));
                corcodeIQ = Xk.* F_CA;
                Cor_Result = abs(ifft(corcodeIQ));
                C(i,:) = C(i,:) + Cor_Result(1,1:length_1ms);
            end
            waitbar(i/(size(freq,2)));
        end

        toc

        [m,ph] = find(C == max(max(C)));
        %Calculate signal power
        P_s = 0;
        if m > size(freq,2) - 2
            up_bound = size(freq,2);
        else
            up_bound = m + 2;
        end

        if m < 3
            low_bound = 1;
        else
            low_bound = m - 2;
        end

        for bin_i = low_bound : up_bound
            if C(bin_i,ph) > 12.5 * P_n
                P_s = P_s + C(bin_i,ph);
            end
        end

        C_N0 = 10*log10(P_s/P_n);

        yes = C_N0 > Threshold;

        if  yes
            %normalize
            C = C./P_n/1e4;
            axispower = max(max(max(C)));
            %find frequency
            f = fc_lo + (m-1)*step;
            % find code phase
            code_phase = ph;

            %             figure(100*sat_id_vector(sat_i) + n(n_num));
            %             subplot(311);
            if ph < csample
                D = C(m,length_1ms-(csample-ph):length_1ms);
                D = [D C(m,1:ph+csample)];
                CM = C(:,length_1ms-(csample-ph):length_1ms);
                CM = [CM C(:,1:ph+csample)];
            elseif ph > length_1ms - csample
                D = C(m,ph-csample:length_1ms);
                D = [D C(m,1:csample-(length_1ms-ph))];
                CM = C(:,ph-csample:length_1ms);
                CM = [CM C(:,1:csample-(length_1ms-ph))];
            else
                D = C(m,ph-csample:ph+csample);
                CM = C(:,ph-csample:ph+csample);
            end
            figure(14);
            subplot(plot_col,3,sat_i);
            xaxis=-csample/16:1/16:csample/16;
            yaxis=((fc_lo-fc):step:(fc_hi-fc));
            mesh(xaxis,yaxis,CM);
            hold on;
            mesh(csample/16:1/16:(csample+1)/16,yaxis,CM(:,csample:csample+1));
            hold on;
            if m == 1
                mesh(xaxis,(fc_hi-fc):(fc_hi-fc)+1,CM(m:m+1,:));
            else
                mesh(xaxis,(fc_hi-fc):(fc_hi-fc)+1,CM(m-1:m,:));
            end
            minipower = min(min(min(CM)));
            axis([-csample/16 csample/16+1 (fc_lo-fc) (fc_hi-fc)+1 minipower axispower]);
            title(['PRN = ' num2str(sat_id_vector(1,sat_i)) ' C/N_0 = ' num2str(C_N0) 'dB Doppler = ' num2str(f - fc) 'Hz' ] );
            xlabel('Code phase (chip)');
            ylabel('Doppler (Hz)');
            zlabel('Magnitude');

        else
            f = 0;
            code_phase = 0;
        end
        close(bar1);
    end
    % check threshold
    fclose(fid);
end