www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/testCAUC3.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 = 36000;%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=2:2 %% 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 = 1.25e6; %[Hz] settings.samplingFreq = 7.5e6; %[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 = 2; %[Hz] settings.dllCorrelatorSpacing = 0.1; %[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); % [carrFreq,codePhase]=Wrelax3(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_CAUC3_1(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