www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/tracking_CAUC3_1.m
function [trackResults, channel,ttime]= tracking_CAUC3_1(p,fid, channel, settings,ttime)
% Performs code and carrier tracking for all channels.
% wenguihong amend
%[trackResults, channel] = tracking(fid, channel, settings)
%
% Inputs:
% fid - file identifier of the signal record.
% channel - PRN, carrier frequencies and code phases of all
% satellites to be tracked (prepared by preRum.m from
% acquisition results).
% settings - receiver settings.
% Outputs:
% trackResults - tracking results (structure array). Contains
% in-phase prompt outputs and absolute starting
% positions of spreading codes, together with other
% observation data from the tracking loops. All are
% saved every millisecond.
%--------------------------------------------------------------------------
% SoftGNSS v3.0
%
% Copyright (C) Dennis M. Akos
% Written by Darius Plausinaitis and Dennis M. Akos
% Based on code by DMAkos Oct-1999
%--------------------------------------------------------------------------
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
%USA.
%--------------------------------------------------------------------------
%CVS record:
%$Id: tracking.m,v 1.14.2.32 2007/01/30 09:45:12 dpl Exp $
%% Initialize result structure ============================================
a=inf;
skip=0;
samplesPerCode = round(settings.samplingFreq / ...
(settings.codeFreqBasis / settings.codeLength));
% Channel status
trackResults.status = '-'; % No tracked signal, or lost lock
% The absolute sample in the record of the C/A code start:
trackResults.absoluteSample = zeros(1, settings.msToProcess);
% Freq of the C/A code:
trackResults.codeFreq = inf(1, settings.msToProcess);
% Frequency of the tracked carrier wave:
trackResults.carrFreq = inf(1, settings.msToProcess);
% Outputs from the correlators (In-phase):
trackResults.I_P = zeros(1, settings.msToProcess);
trackResults.I_E = zeros(1, settings.msToProcess);
trackResults.I_L = zeros(1, settings.msToProcess);
% Outputs from the correlators (Quadrature-phase):
trackResults.Q_E = zeros(1, settings.msToProcess);
trackResults.Q_P = zeros(1, settings.msToProcess);
trackResults.Q_L = zeros(1, settings.msToProcess);
% Loop discriminators
trackResults.dllDiscr = inf(1, settings.msToProcess);
trackResults.dllDiscrFilt = inf(1, settings.msToProcess);
trackResults.pllDiscr = inf(1, settings.msToProcess);
trackResults.pllDiscrFilt = inf(1, settings.msToProcess);
%--- Copy initial settings for all channels -------------------------------
trackResults = repmat(trackResults, 1, settings.numberOfChannels);
%% Initialize tracking variables ==========================================
codePeriods = settings.msToProcess; % For GPS one C/A code is one ms
%--- DLL variables --------------------------------------------------------
% Define early-late offset (in chips)
earlyLateSpc = settings.dllCorrelatorSpacing;
% Summation interval
PDIcode = 0.001;
% Calculate filter coefficient values
[tau1code, tau2code] = calcLoopCoef(settings.dllNoiseBandwidth, ...
settings.dllDampingRatio, ...
1.0);
%--- PLL variables --------------------------------------------------------
% Summation interval
PDIcarr = 0.001;
% Calculate filter coefficient values
[tau1carr, tau2carr] = calcLoopCoef(settings.pllNoiseBandwidth, ...
settings.pllDampingRatio, ...
0.25);
% hwb = waitbar(0,'Tracking...');
N=500;
%% Start processing channels ==============================================
for channelNr = 1:settings.numberOfChannels
% Only process if PRN is non zero (acquisition was successful)
if (channel(channelNr).PRN ~= 0)
% Save additional information - each channel's tracked PRN
trackResults(channelNr).PRN = channel(channelNr).PRN;
% Move the starting point of processing. Can be used to start the
% signal processing at any point in the data record (e.g. for long
% records). In addition skip through that data file to start at the
% appropriate sample (corresponding to code phase). Assumes sample
% type is schar (or 1 byte per sample)
fseek(fid, ...
settings.skipNumberOfBytes + channel(channelNr).codePhase-1, ...
'bof');
% Get a vector with the C/A code sampled 1x/chip
caCode = generateCAcode(channel(channelNr).PRN);
% Then make it possible to do early and late versions
caCode = [caCode(1023) caCode caCode(1)];% caCode(1)
%--- Perform various initializations ------------------------------
% define initial code frequency basis of NCO
codeFreq = settings.codeFreqBasis;
% define residual code phase (in chips)
remCodePhase = 0.0;
% define carrier frequency which is used over whole tracking period
carrFreq = channel(channelNr).acquiredFreq;
carrFreqBasis = channel(channelNr).acquiredFreq;
% define residual carrier phase
remCarrPhase = 0.0;
%code tracking loop parameters
oldCodeNco = 0.0;
oldCodeError = 0.0;
%carrier/Costas loop parameters
oldCarrNco = 0.0;
oldCarrError = 0.0;
cPhaseold=channel(channelNr).codePhase;
%=== Process the number of specified code periods =================
for loopCnt = 1:codePeriods
%% GUI update -------------------------------------------------------------
% The GUI is updated every 50ms. This way Matlab GUI is still
% responsive enough. At the same time Matlab is not occupied
% all the time with GUI task.
% if (rem(loopCnt, 50) == 0)
% try
% waitbar(loopCnt/codePeriods, ...
% hwb, ...
% ['Tracking: Ch ', int2str(channelNr), ...
% ' of ', int2str(settings.numberOfChannels), ...
% '; PRN#', int2str(channel(channelNr).PRN), ...
% '; Completed ',int2str(loopCnt), ...
% ' of ', int2str(codePeriods), ' msec']);
% catch
%
% % The progress bar was closed. It is used as a signal
% % to stop, "cancel" processing. Exit.
% disp('Progress bar closed, exiting...');
% return
% end
% end
%% Read next block of data ------------------------------------------------
% Find the size of a "block" or code period in whole samples
% Update the phasestep based on code freq (variable) and
% sampling frequency (fixed)
codePhaseStep = codeFreq / settings.samplingFreq;
% if channelNr==2
% [codePhase,carrFreq]=Wrelax1(p,settings,loopCnt);
% blksize = ceil(codePhase-codePhaseold+6000);
% codePhaseold=codePhase;
% [rawSignal, samplesRead] = fread(fid, blksize, settings.dataType);
% blksize
% trackResults(channelNr).carrFreq(loopCnt) = carrFreq;
% trackResults(channelNr).codeFreq(loopCnt) = codeFreq;
% else
blksize = ceil((settings.codeLength-remCodePhase) / codePhaseStep);
% end
% Read in the appropriate number of samples to process this
% interation
[rawSignal, samplesRead] = fread(fid, blksize, settings.dataType);
%% for bupt if data
rawSignal =rawSignal.';
% If did not read in enough samples, then could be out of
% data - better exit
if (samplesRead ~= blksize)
disp('Not able to read the specified number of samples for tracking, exiting!')
fclose(fid);
return
end
%% 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;
a=find(tcode2>1025);
b=find(tcode2<1);
while tcode2(a)>1025
tcode2(a)=tcode2(a)-1;
a=find(tcode2>1025);
end
while tcode2(b)<1
tcode2(b)=tcode2(b)+1;
b=find(tcode2<1);
end
earlyCode = caCode(tcode2);
% Define index into late code vector
tcode = (remCodePhase+earlyLateSpc) : ...
codePhaseStep : ...
((blksize-1)*codePhaseStep+remCodePhase+earlyLateSpc);
tcode2 = ceil(tcode) + 1;
a=find(tcode2>1025);
b=find(tcode2<1);
while tcode2(a)>1025
tcode2(a)=tcode2(a)-1;
a=find(tcode2>1025);
end
while tcode2(b)<1
tcode2(b)=tcode2(b)+1;
b=find(tcode2<1);
end
lateCode = caCode(tcode2);
% Define index into prompt code vector
tcode = remCodePhase : ...
codePhaseStep : ...
((blksize-1)*codePhaseStep+remCodePhase);
tcode2 = ceil(tcode) + 1;
a=find(tcode2>1025);
b=find(tcode2<1);
while tcode2(a)>1025
tcode2(a)=tcode2(a)-1;
a=find(tcode2>1025);
end
while tcode2(b)<1
tcode2(b)=tcode2(b)+1;
b=find(tcode2<1);
end
promptCode = caCode(tcode2);
remCodePhase = (tcode(blksize) + codePhaseStep) - 1023.0;
%% Generate the carrier frequency to mix the signal to baseband -----------
time = (0:blksize) ./ settings.samplingFreq;
% Get the argument to sin/cos functions
trigarg = ((carrFreq * 2.0 * pi) .* time) + remCarrPhase;
remCarrPhase = rem(trigarg(blksize+1), (2 * pi));
% 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 PLL error and update carrier NCO ----------------------------------
% Implement carrier loop discriminator (phase detector)
carrError = atan(Q_P / I_P) / (2.0 * pi);
% Implement carrier loop filter and generate NCO command
carrNco = oldCarrNco + (tau2carr/tau1carr) * ...
(carrError - oldCarrError) + carrError * (PDIcarr/tau1carr);
oldCarrNco = carrNco;
oldCarrError = carrError;
% Modify carrier freq based on NCO command
carrFreq = carrFreqBasis + carrNco;
trackResults(channelNr).carrFreq(loopCnt) = carrFreq;
%% 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 = settings.codeFreqBasis - codeNco;
trackResults(channelNr).codeFreq(loopCnt) = codeFreq;
%% Record various measures to show in postprocessing ----------------------
% Record sample number (based on 8bit samples)
if channel(channelNr).PRN==13
% if mod((loopCnt-1),N)==0
% [cPhase,cFreq]=Wrelax1(p,settings,cPhaseold,skip,loopCnt);
if loopCnt==1
% [t1_p,fd1]= tracking_CAUC_loop(channel,channelNr,settings,p,loopCnt);
% t1_p=ftell(fid);
% t1=mod(t1_p,samplesPerCode)/samplesPerCode*1e-3;
fd1=carrFreq;
else
% t1=cPhase;
fd1=cFreq;
end
t1_p=ftell(fid);
[cPhase,cFreq]=Wrelax3_jinggu_1(p,settings,loopCnt,fd1);
% [cPhase,cFreq]=Wrelax2(p,settings,loopCnt);
trackResults(channelNr).absoluteSample(loopCnt) = cPhase+samplesPerCode*loopCnt;
% blk=trackResults(channelNr).absoluteSample(loopCnt);
% if loopCnt==1
% trackResults(channelNr).absoluteSample(loopCnt) = cPhase+samplesPerCode;
% else
% trackResults(channelNr).absoluteSample(loopCnt) = round(trackResults(channelNr).absoluteSample(loopCnt-1)+cPhase);
% end
if abs(trackResults(channelNr).absoluteSample(loopCnt)-t1_p)>5
trackResults(channelNr).absoluteSample(loopCnt)=t1_p;
end
trackResults(channelNr).codeFreq(loopCnt) = cFreq;
%%%%%%%%%
% blk = cPhase-cPhaseold+samplesPerCode*N;
% blk1=cPhase-cPhaseold+samplesPerCode;
% if abs(blk-samplesPerCode*N)>60
% blk=blkold;
% end
% % blk=samplesPerCode-cPhase;
% cPhaseold=cPhase;
% blkold=blk;
%%%%%%%%%
loopCnt
if mod(loopCnt,50)==0
clc
end
% fd1
% cFreq
% cFreq-fd1
% t1_p
% blk
% blk-t1_p
% blksize
%%%%%%%%%%%
% if loopCnt==1
% trackResults(channelNr).absoluteSample(loopCnt) = settings.skipNumberOfBytes + channel(channelNr).codePhase-1+blk1;
% else
% trackResults(channelNr).absoluteSample(loopCnt) = trackResults(channelNr).absoluteSample(loopCnt-N)+blk;
% end
%%%%%%%%%%%%
% skip=trackResults(channelNr).absoluteSample(loopCnt);
% skip
% else
% trackResults(channelNr).absoluteSample(loopCnt) = round(trackResults(channelNr).absoluteSample(loopCnt-1)+blk/N);
% if loopCnt~=1
% for ii=1:N-1
% trackResults(channelNr).absoluteSample(loopCnt-N+ii) = trackResults(channelNr).absoluteSample(loopCnt-N-1+ii)+...
% round((trackResults(channelNr).absoluteSample(loopCnt)-trackResults(channelNr).absoluteSample(loopCnt-N))/(N-1));
% end
% end
% end
% if mod((loopCnt-1),10)==0
% % [cPhase,cFreq]=Wrelax1(p,settings,cPhaseold,skip,loopCnt);
% % [cPhase,cFreq]=Wrelax1(p,settings,loopCnt);
% [cPhase,cFreq]=Wrelax2(p,settings,loopCnt);
% blk = cPhase-cPhaseold+samplesPerCode*10;
% blk1=cPhase-cPhaseold+samplesPerCode;
% if abs(blk-samplesPerCode*10)>30
% blk=blkold;
% end
% % blk=samplesPerCode-cPhase;
% cPhaseold=cPhase;
% blkold=blk;
% blk
% blksize
% loopCnt
% if loopCnt==1
% trackResults(channelNr).absoluteSample(loopCnt) = settings.skipNumberOfBytes + channel(channelNr).codePhase-1+blk1;
% else
% trackResults(channelNr).absoluteSample(loopCnt) = trackResults(channelNr).absoluteSample(loopCnt-10)+blk;
% end
%
% else
% trackResults(channelNr).absoluteSample(loopCnt) = trackResults(channelNr).absoluteSample(loopCnt-1)+blk/10;
% end
else
% end
trackResults(channelNr).absoluteSample(loopCnt) = ftell(fid);
end
trackResults(channelNr).dllDiscr(loopCnt) = codeError;
trackResults(channelNr).dllDiscrFilt(loopCnt) = codeNco;
trackResults(channelNr).pllDiscr(loopCnt) = carrError;
trackResults(channelNr).pllDiscrFilt(loopCnt) = carrNco;
trackResults(channelNr).I_E(loopCnt) = I_E;
trackResults(channelNr).I_P(loopCnt) = I_P;
trackResults(channelNr).I_L(loopCnt) = I_L;
trackResults(channelNr).Q_E(loopCnt) = Q_E;
trackResults(channelNr).Q_P(loopCnt) = Q_P;
trackResults(channelNr).Q_L(loopCnt) = Q_L;
end % for loopCnt
if channel(channelNr).PRN==13
ttime(:,p)=trackResults(channelNr).absoluteSample(loopCnt);
end
% If we got so far, this means that the tracking was successful
% Now we only copy status, but it can be update by a lock detector
% if implemented
trackResults(channelNr).status = channel(channelNr).status;
end % if a PRN is assigned
end % for channelNr
save result.mat trackResults channel
% Close the waitbar
% close(hwb)