matlab_GPS信号产生、捕获、追踪全套程序

function codeError= test(PRN,data) format long; % Outputs: % trackResults - tracking results (structure array). Contains % in-phase prompt outputs and absolute spreading % code's starting positions, together with other % observation data from the tracking loops. All are % saved every millisecond. %% Initialize result structure ============================================ chiprate = 1.023e6; %码速率 IF = 9.548e6; %载波频率 samplingFreq = 38.192e6; %采样频率 codeLength = 1023; %码长 msToProcess = 50; %[ms] % Outputs from the correlators (In-phase): trackResults.I_P = zeros(1, msToProcess); trackResults.I_E = zeros(1, msToProcess); trackResults.I_L = zeros(1, msToProcess); % Outputs from the correlators (Quadrature-phase): trackResults.Q_E = zeros(1, msToProcess); trackResults.Q_P = zeros(1, msToProcess); trackResults.Q_L = zeros(1, msToProcess); % Loop discriminators trackResults.dllDiscr = 0.4.*ones(1, msToProcess); trackResults.dllDiscrFilt = inf(1, msToProcess); %% Initialize tracking variables ========================================== codePeriods = msToProcess; % For GPS one C/A code is one ms %--- DLL variables -------------------------------------------------------- % Define early-late offset (in chips) earlyLateSpc = 0.5; % Summation interval PDIcode = 100; % Code tracking loop parameters settings.dllDampingRatio = 0.7; settings.dllNoiseBandwidth = 2; %[Hz] settings.dllCorrelatorSpacing = 0.5; %[chips] % Calculate filter coefficient values [tau1code, tau2code] = calcLoopCoef(settings.dllNoiseBandwidth, ... settings.dllDampingRatio, ... 1.0); %% Start processing channels ============================================== % Get a vector with the C/A code sampled 1x/chip caCode = cacode(PRN); % Then make it possible to do early and late versions caCode = [caCode(1023) caCode caCode(1)]; % define initial code frequency basis of NCO codeFreq = chiprate; % define residual code phase (in chips) remCodePhase = 0.5; %code tracking loop parameters oldCodeNco = 0.0; oldCodeError = 0.0; %=== Process the number of specified code periods ================= for loopCnt = 1:codePeriods %% Read next block of data ------------------------------------------------ codePhaseStep = codeFreq / samplingFreq; %1.023e6/38.192e6 blksize = ceil((codeLength) / (1.023e6/38.192e6)); %38192个数 % Read in the appropriate number of samples to process this % interation rawSignal = data; %transpose vector %% Set up all the code phase tracking information ------------------------- % Define index into early code vector tcode = (remCodePhase-earlyLateSpc) : ... codePhaseStep : ... ((blksize-1)*codePhaseStep+remCodePhase-earlyLateSpc); tcode2 = ceil(tcode) + 1; earlyCode = caCode(tcode2); %为什么加1呢？我验证了，会使得下标达到1024,因为这不是cacode， %而是caCode，定义变成了caCode = [caCode(1023) caCode caCode(1)] % Define index into late code vector tcode = (remCodePhase+earlyLateSpc) : ... codePhaseStep : ... ((blksize-1)*codePhaseStep+remCodePhase+earlyLateSpc); tcode2 = ceil(tcode) + 1; lateCode = caCode(tcode2); % Define index into prompt code vector tcode = remCodePhase : ... codePhaseStep : ... ((blksize-1)*codePhaseStep+remCodePhase); tcode2 = ceil(tcode) + 1; promptCode = caCode(tcode2); %% Generate the carrier frequency to mix the signal to baseband ----------- time = (0:blksize-1) ./ samplingFreq; % Get the argument to sin/cos functions trigarg = ((IF * 2.0 * pi) .* time) ; % Finally compute the signal to mix the collected data to bandband carrCos = cos(trigarg(1:blksize)); carrSin = sin(trigarg(1:blksize)); %% Generate the six standard accumulated values --------------------------- % First mix to baseband qBasebandSignal = carrCos .* rawSignal; iBasebandSignal = carrSin .* rawSignal; % Now get early, late, and prompt values for each I_E = sum(earlyCode .* iBasebandSignal); Q_E = sum(earlyCode .* qBasebandSignal); I_P = sum(promptCode .* iBasebandSignal); Q_P = sum(promptCode .* qBasebandSignal); I_L = sum(lateCode .* iBasebandSignal); Q_L = sum(lateCode .* qBasebandSignal); %% Find DLL error and update code NCO ------------------------------------- codeError = (sqrt(I_E * I_E + Q_E * Q_E) - sqrt(I_L * I_L + Q_L * Q_L)) / ... (sqrt(I_E * I_E + Q_E * Q_E) + sqrt(I_L * I_L + Q_L * Q_L)); % Implement code loop filter and generate NCO command codeNco = oldCodeNco + (tau2code/tau1code) * ... (codeError - oldCodeError) + codeError * (PDIcode/tau1code); oldCodeNco = codeNco; oldCodeError = codeError; % Modify code freq based on NCO command codeFreq = chiprate - codeNco; trackResults.codeFreq(loopCnt) = codeFreq; %% Record various measures to show in postprocessing ---------------------- % Record sample number (based on 8bit samples) trackResults.dllDiscr(loopCnt) = codeError; trackResults.dllDiscrFilt(loopCnt) = codeNco; trackResults.I_E(loopCnt) = I_E; trackResults.I_P(loopCnt) = I_P; trackResults.I_L(loopCnt) = I_L; trackResults.Q_E(loopCnt) = Q_E; trackResults.Q_P(loopCnt) = Q_P; trackResults.Q_L(loopCnt) = Q_L; end % for loopCnt figure(8); plot(trackResults.dllDiscr);