www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/testCAUC_DCT.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 = 8; %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:1
%% 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
filename1=['E:\realDN_' int2str(p) '.dat'];
settings.fileName =filename1;
settings.dataType ='double';
% Intermediate, sampling and code frequencies
settings.IF = 1.25e6; %[Hz]
settings.samplingFreq = 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 = 2.5;
%% Tracking loops settings =======================================
% Code tracking loop parameters
settings.dllDampingRatio = 0.7;
settings.dllNoiseBandwidth = 2; %[Hz]
settings.dllCorrelatorSpacing = 0.5; %[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, 280*samplesPerCode, settings.dataType);
% data = data.';
% fclose(fid);
%
% % bb=fir1(128,[0.331,0.531]);
% % freqz(bb,1,512)
% % y2 = filter(bb,1,data);
%
%
% 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.';
% a1=dct(data);
% a1(abs(a1)<100)=0;
% K = idct(a1);
% plot(abs(fft(K-data)))
% % data=K;
% data=data-K;
% ma=max(abs(real(data)));
% data=round(real(data)*128/ma);
%--- Do the acquisition -------------------------------------------
disp ('Acquiring satellites...');
tic
acqResults = acquisition(data, settings);
toc
% acqResults.codePhase(1)=4315;
plotAcquisition1(acqResults);
end
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% phase(p)=acqResults.codePhase(13);
%% 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 = preRun(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,ttime] = tracking_CAUC(fid, channel, settings,p,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 = postNavigation(trackResults, settings);
% navSolutions = postNavigation2(p,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;
if a(1)==39
laterror=lat(:,3)*31;
else
laterror=(60-lat(:,3))*31;
end
if b(1)==50
longerror=cos(39/180*pi)*long(:,3)*31;
else
longerror=cos(39/180*pi)*(60-long(:,3))*31;
end
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);
d1=sqrt(laterror.^2+longerror.^2);
end
mse=sqrt(sum(d.^2)/p);
mse1=sqrt(sum(d1.^2)/p);
% save r3.mat lat long hgt laterror longerror
% ti=toc