www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/testCAUC1.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% used to test the CAUC GPS IF data
% by Haitao Liu
% OK
% the sample data must be convert by "cauc_convert_gps_if_data.m"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by Haitao Liu
% 2009.06.06
% used to test the IF data
%%------------------------------------------------------------------------
clear all
clc;
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath include % The software receiver functions
addpath geoFunctions % Position calculation related functions
ttime=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Processing settings ===========================================
% Number of milliseconds to be processed used 36000 + any transients (see
% below - in Nav parameters) to ensure nav subframes are provided
settings.msToProcess = 19;%36100; %[ms]
% Number of channels to be used for signal processing
settings.numberOfChannels = 4; %8
% 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). fseek
% function is used to move the file read point, therefore advance is byte
% based only.
settings.skipNumberOfBytes =0;
for p=1:20
%% Raw signal file name and other parameter ======================
% This is a "default" name of the data file (signal record) to be used in
% the post-processing mode
% file3=['result_f' int2str(p) '.mat'];
% file4=['result_c' int2str(p) '.mat'];
% fid3=fopen(file3,'r');
% fid4=fopen(file4,'r');
% acqResults.carrFreq=fread(fid3,'single');
% acqResults.codePhase=fread(fid4,'single');
% acqResults.carrFreq=round(acqResults.carrFreq).';
% acqResults.codePhase=round(acqResults.codePhase).';
% % acqResults.peakMetric=zeros(1,32);
% fclose(fid3);
% fclose(fid4);
filename1=['g:\real_' int2str(p) '.dat'];
settings.fileName =filename1;
settings.dataType ='int8';
% Intermediate, sampling and code frequencies
settings.IF = 4.309e6; %[Hz]
settings.samplingFreq = 20e6; %[Hz]
settings.codeFreqBasis = 1.023e6; %[Hz]
% Define number of chips in a code period
settings.codeLength = 1023;
%% Acquisition settings ==========================================
% Skips acquisition in the script postProcessing.m if set to 1
settings.skipAcquisition = 0;
% List of satellites to look for. Some satellites can be excluded to speed
% up acquisition
settings.acqSatelliteList = 1:32; %[PRN4 numbers]
% Band around IF to search for satellite signal. Depends on max Doppler
settings.acqSearchBand = 10; %[kHz]
% Threshold for the signal presence decision rule
settings.acqThreshold = 1.8;
%% Tracking loops settings =======================================
% Code tracking loop parameters
settings.dllDampingRatio = 0.7;
settings.dllNoiseBandwidth = 10; %[Hz]
settings.dllCorrelatorSpacing = 0.2; %[chips]
% Carrier tracking loop parameters
settings.pllDampingRatio = 0.7;
settings.pllNoiseBandwidth = 25; %[Hz]
%% Navigation solution settings ==================================
% Period for calculating pseudoranges and position
settings.navSolPeriod = 500; %[ms]
% Elevation mask to exclude signals from satellites at low elevation
settings.elevationMask = 10; %[degrees 0 - 90]
% Enable/dissable use of tropospheric correction
settings.useTropCorr = 1; % 0 - Off
% 1 - On
% True position of the antenna in UTM system (if known). Otherwise enter
% all NaN's and mean position will be used as a reference .
settings.truePosition.E = nan;
settings.truePosition.N = nan;
settings.truePosition.U = nan;
%% Plot settings =================================================
% Enable/disable plotting of the tracking results for each channel
settings.plotTracking = 1; % 0 - Off
% 1 - On
%% Constants =====================================================
settings.c = 299792458; % The speed of light, [m/s]
settings.startOffset = 68.802; %[ms] Initial sign. travel time
%%------------------------------------------------------------------------
[fid, message] = fopen(settings.fileName ,'r');
if (fid > 0)
% 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).
fseek(fid, settings.skipNumberOfBytes, 'bof');
% Find number of samples per spreading code
samplesPerCode = round(settings.samplingFreq/(settings.codeFreqBasis / settings.codeLength));
% Read 11ms of signal
[data, count] = fread(fid, 11*samplesPerCode, settings.dataType);
data = data.';
fclose(fid);
if (count < 11*samplesPerCode)
% The file is to short
error('Could not read enough data from the data file.');
end
%----------------------------------------------------------------------
% display the IF data
% plotspec(data-mean(data),1/settings.samplingFreq);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% acquisition the satelliate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initialization ================================================
disp ('Starting acquisiting ...');
[fid, message] = fopen(settings.fileName ,'r','b');
% If success, then process the data
if (fid > 0)
% Move the starting point of processing. Can be used to start the
% signal processing at any point in the data record (e.g. good for long
% records or for signal processing in blocks).
fseek(fid, settings.skipNumberOfBytes, 'bof');
%% Acquisition ===================================================
% Do acquisition if it is not disabled in settings or if the variable
% acqResults does not exist.
if (settings.skipAcquisition == 0)
% Find number of samples per spreading code
samplesPerCode = round(settings.samplingFreq /(settings.codeFreqBasis / settings.codeLength));
% Read data for acquisition. 11ms of signal are needed for the fine
% frequency estimation
data = fread(fid, 11*samplesPerCode, settings.dataType);
data = data.';
% data=round(data*2^24);
%--- Do the acquisition -------------------------------------------
disp ('Acquiring satellites...');
acqResults = acquisition(data, settings);
% plotAcquisition1(acqResults);
end
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[carrFreq,codePhase]=Wrelax(p,settings);
acqResults.carrFreq(13)=carrFreq;
acqResults.codePhase(13)=codePhase;
% acqResults=Wrelax(p,settings);
%% Initialize channels and prepare for the run ============================
% Start further processing only if a GNSS signal was acquired (the
% field FREQUENCY will be set to 0 for all not acquired signals)
if (any(acqResults.carrFreq))
channel = preRun1(acqResults, settings);
showChannelStatus(channel, settings);
else
% No satellites to track, exit
disp('No GNSS signals detected, signal processing finished.');
trackResults = [];
return;
end
%% Track the signal =======================================================
startTime = now;
disp ([' Tracking started at ', datestr(startTime)]);
% Process all channels for given data block
% [trackResults, channel] = tracking_CAUC(fid, channel, settings);
[trackResults, channel,ttime] = tracking_CAUC1(p,fid, channel, settings,ttime);
% Close the data file
fclose(fid);
disp([' Tracking is over (elapsed time ', ...
datestr(now - startTime, 13), ')'])
% Auto save the acquisition & tracking results to a file to allow
% running the positioning solution afterwards.
disp(' Saving Acq & Tracking results to file "trackingResults.mat"')
save('trackingResults', ...
'trackResults', 'settings', 'acqResults', 'channel');
%% Calculate navigation solutions =========================================
disp(' Calculating navigation solutions...');
% navSolutions = postNavigation1(trackResults, settings, p);
navSolutions = postNavigation(trackResults, settings);
disp(' Processing is complete for this data block');
%% Plot all results ===================================================
disp ('Ploting results...');
if settings.plotTracking
% plotTracking(1:settings.numberOfChannels, trackResults, settings);
end
% plotNavigation(navSolutions, settings);
disp('Post processing of the signal is over.');
load temp.mat
lat(p,:)=a;
long(p,:)=b;
hgt(p,:)=c;
laterror=lat(:,3)*31;
longerror=cos(39/180*pi)*long(:,3)*31;
x1=(39+0/60.0+0/3600.0)/180.0*pi;
y1=(50+0/60.0+0/3600.0)/180.0*pi;
x2=(lat(p,1)+lat(p,2)/60.0+lat(p,3)/3600.0)/180.0*pi;
y2=(long(p,1)+long(p,2)/60.0+long(p,3)/3600.0)/180.0*pi;
AB=sqrt(2*(1-cos(abs(x2-x1))+cos(x1)*cos(x2)-cos(x1)*cos(x2)*cos(abs(y2-y1))));
d(p,:)=2*6371e3*asin(AB/2);
end
% save r2.mat lat long hgt laterror longerror
% ti=toc