diff --git a/BDS/B1C/include/acquisition.m b/BDS/B1C/include/acquisition.m index b259f2f..5334d4f 100644 --- a/BDS/B1C/include/acquisition.m +++ b/BDS/B1C/include/acquisition.m @@ -168,7 +168,15 @@ %--- Input signal power for GLRT statistic calculation -------------------- sigPower = sqrt(var(sig10PlusXms(1:samplesXmsLen)) * samplesXmsLen); - + %--- Generate carrier wave frequency grid ----------------------- +initFreq = settings.IF - settings.acqSearchBand; +% Generate local sine and cosine +sigCarr = exp(1i*initFreq* phasePoints); +% "Remove carrier" from the signal +I1 = real(sigCarr .* sig10PlusXms); +Q1 = imag(sigCarr .* sig10PlusXms); +% Convert the baseband signal to frequency domain +IQfreqDom = fft(I1 + 1i*Q1); % Perform search for all listed PRN numbers ... fprintf('('); for PRN = settings.acqSatelliteList @@ -190,27 +198,17 @@ zeros(1,len10PlusXms - samplesXmsLen)]; PilotPriFreqDom = conj(fft(localPilot)); end - %--- Make the correlation for whole frequency band (for all freq. bins) for frqBinIndex = 1:numberOfFrqBins - %--- Generate carrier wave frequency grid ----------------------- - frqBins(frqBinIndex) = settings.IF - settings.acqSearchBand + ... - settings.acqStep * (frqBinIndex - 1); - % Generate local sine and cosine - sigCarr = exp(1i*frqBins(frqBinIndex) * phasePoints); - % "Remove carrier" from the signal - I1 = real(sigCarr .* sig10PlusXms); - Q1 = imag(sigCarr .* sig10PlusXms); - % Convert the baseband signal to frequency domain - IQfreqDom = fft(I1 + 1i*Q1); + IQfreqDomShifted = circshift(IQfreqDom, frqBinIndex - 1); % Multiplication in the frequency domain (correlation in time domain) - convCodeIQ1 = IQfreqDom .* DataPriFreqDom; + convCodeIQ1 = IQfreqDomShifted .* DataPriFreqDom; % Perform inverse DFT and store correlation results results(frqBinIndex, :) = abs(ifft(convCodeIQ1)); % Pilot signal acquisition if (settings.pilotACQflag == 1) - convCodeIQ1 = IQfreqDom .* PilotPriFreqDom; + convCodeIQ1 = IQfreqDomShifted .* PilotPriFreqDom; % Non-coherent combining of data and pilot results results(frqBinIndex, :) = (results(frqBinIndex, :)* sqrt(11)+ ... abs(ifft(convCodeIQ1))* sqrt(29) )/ sqrt(40); @@ -224,6 +222,7 @@ % Find the correlation peak and the carrier frequency [~, frequencyBinIndex] = max(max(results, [], 2)); + selFreq = initFreq + (frequencyBinIndex -1)*settings.acqStep; % Find code phase of the same correlation peak [peakSize, codePhase] = max(max(results)); % Store GLRT statistic @@ -253,8 +252,7 @@ %--- Search different frequency bins ------------------------------ for FineBinIndex = 1 : NumOfFineBins % Carrier frequencies of the fine frequency bins - FineFrqBins(FineBinIndex) = frqBins(frequencyBinIndex) -... - settings.acqStep + fineStep * (FineBinIndex - 1); + FineFrqBins(FineBinIndex) = selFreq - settings.acqStep + fineStep * (FineBinIndex - 1); % Generate local sine and cosine sigCarr10cm = exp(1i*FineFrqBins(FineBinIndex) * finePhasePoints); FineResult(FineBinIndex) = abs(sum(xCarrier .* sigCarr10cm)); diff --git a/BDS/B1C/include/acquisitionCS.m b/BDS/B1C/include/acquisitionCS.m deleted file mode 100644 index fdb0fef..0000000 --- a/BDS/B1C/include/acquisitionCS.m +++ /dev/null @@ -1,307 +0,0 @@ -function acqResults = acquisitionCS(longSignal, settings) -%Function performs cold start acquisition on the collected "data". It -%searches for GPS signals of all satellites, which are listed in field -%"acqSatelliteList" in the settings structure. Function saves code phase -%and frequency of the detected signals in the "acqResults" structure. -% -%acqResults = acquisition(longSignal, settings) -% -% Inputs: -% longSignal - 20 ms of raw IF signal from the front-end.. -% settings - Receiver settings. Provides information about -% sampling and intermediate frequencies and other -% parameters including the list of the satellites to -% be acquired. -% Outputs: -% acqResults - Function saves code phases and frequencies of the -% detected signals in the "acqResults" structure. The -% field "carrFreq" is set to 0 if the signal is not -% detected for the given PRN number. - -%-------------------------------------------------------------------------- -% CU Multi-GNSS SDR -% (C) Developed for BDS B1C SDR by Yafeng Li, Nagaraj C. Shivaramaiah -% and Dennis M. Akos. -% Based on the original framework for GPS C/A SDR by Darius Plausinaitis, -% Peter Rinder, Nicolaj Bertelsen and Dennis M. Akos -% -% Reference: Li, Y., Shivaramaiah, N.C. & Akos, D.M. Design and -% implementation of an open-source BDS-3 B1C/B2a SDR receiver. -% GPS Solut (2019) 23: 60. https://doi.org/10.1007/s10291-019-0853-z -%-------------------------------------------------------------------------- -%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: acquisition.m,v 1.1.2.12 2006/08/14 12:08:03 dpl Exp $ - - -%% Condition input signal to speed up acquisition =================== - -% If input IF signal freq. is too high, a resampling strategy is applied -% to speed up the acquisition. This is user selectable. -if (settings.samplingFreq > settings.resamplingThreshold && ... - settings.resamplingflag == 1) - - %--- Filiter out signal power outside the main lobe of CM code -------- - fs = settings.samplingFreq; - IF = settings.IF; - % Bandwidth of CM mian lobe - BW = 9e6; - % Filter parameter - w1 = (IF)-BW/2; - w2 = (IF)+BW/2; - wp = [w1*2/fs-0.002 w2*2/fs+0.002]; - % Filter coefficients - b = fir1(700,wp); - % Filter operation - longSignal = filtfilt(b,1,longSignal); - - % --- Find resample frequency ----------------------------------------- - % Refer to bandpass sampling theorem(Yi-Ran Sun,Generalized Bandpass - % Sampling Receivers for Software Defined Radio) - - % Upper boundary frequency of the bandpass IF signal - fu = settings.IF + BW/2; - - % Lower freq. of the acceptable sampling Freq. range - n = floor(fu/BW); - if (n<1) - n = 1; - end - lowerFreq = 2*fu/n; - - % Lower boundary frequency of the bandpass IF signal - fl = settings.IF - BW/2; - - % Upper boundary frequency of the acceptable sampling Freq. range - if(n>1) - upperFreq = 2*fl/(n-1); - else - upperFreq = lowerFreq; - end - - % Save orignal Freq. for later use - oldFreq = settings.samplingFreq; - - % Take the center of the acceptable sampling Freq. range as - % resampling frequency. As settings are used to generate local - % CM code samples, so assign the resampling freq. to settings. - % This can not change the settings.samplingFreq outside this - % acquisition function. - settings.samplingFreq = ceil((lowerFreq + upperFreq)/2); - - %--- Downsample input IF signal --------------------------------------- - % Signal length after resampling - signalLen = floor((length(longSignal)-1) /oldFreq * settings.samplingFreq); - % Resampled signal index - index = ceil((0:signalLen-1)/settings.samplingFreq * oldFreq); - index(1) = 1; - % Resampled signal - longSignal = longSignal(index); - - % Foe latter use - oldIF = settings.IF; - - % Resampling is equivalent to down-converting the original IF by integer - % times of resampling freq.. So the IF after resampling is equivalent to: - settings.IF = rem(settings.IF,settings.samplingFreq); - -end % resampling input IF signals - -%% Acquisition initialization ======================================= -%--- Varaibles for coarse acquisition ------------------------------------- -% Number of samples per B1C code -samplesPerCode = round(settings.samplingFreq / ... - (settings.codeFreqBasis / settings.codeLength)); -% Number of samples of settings.acqCohT ms B1C spreading codes -samplesXmsLen= round(samplesPerCode /10*settings.acqCohT); -% Number of samples of 10 plus acqCohT(X) ms B1C spreading codes used to do -% zero padding correlation for coarse acquisition -len10PlusXms = round(samplesPerCode /10*(10+settings.acqCohT)); -% Cut 10 plus X cm input signal to do correlate -sig10PlusXms = longSignal(1:len10PlusXms); -% Find sampling period -ts = 1 / settings.samplingFreq; -% Find phase points of the local carrier wave -phasePoints = (0 : (len10PlusXms-1)) * 2 * pi * ts; - -% Number of the frequency bins for the given acquisition band -numberOfFrqBins = round(settings.acqSearchBand * 2 / settings.acqStep) + 1; -% Search results of all frequency bins and code shifts (for one satellite) -results = zeros(numberOfFrqBins, len10PlusXms); -% Carrier frequencies of the frequency bins -frqBins = zeros(1, numberOfFrqBins); - -%--- Initialize acqResults and related variables -------------------------- -% Carrier frequencies of detected signals -acqResults.carrFreq = zeros(1, max(settings.acqSatelliteList)); -% PRN code phases of detected signals -acqResults.codePhase = zeros(1, max(settings.acqSatelliteList)); -% Correlation peak ratios of the detected signals -acqResults.peakMetric = zeros(1, max(settings.acqSatelliteList)); - -%--- Varaibles for fine acquisition --------------------------------------- -% Number of the frequency bins for fine acquisition: use -% settings.acqStep*2 fine acquisition band, and 25Hz step -fineStep = 25; -NumOfFineBins = round(settings.acqStep/25)*2 + 1; -% Carrier frequencies of the frequency bins -FineFrqBins = zeros(1, NumOfFineBins); -% Non-coherent results of all frequency bins -FineResult = zeros(1,NumOfFineBins); -% Phase points of the local carrier wave -finePhasePoints = (0 : samplesPerCode-1) * 2 * pi * ts; - -%--- Input signal power for GLRT statistic calculation -------------------- -sigPower = sqrt(var(sig10PlusXms(1:samplesXmsLen)) * samplesXmsLen); - %--- Generate carrier wave frequency grid ----------------------- -initFreq = settings.IF - settings.acqSearchBand; -% Generate local sine and cosine -sigCarr = exp(1i*initFreq* phasePoints); -% "Remove carrier" from the signal -I1 = real(sigCarr .* sig10PlusXms); -Q1 = imag(sigCarr .* sig10PlusXms); -% Convert the baseband signal to frequency domain -IQfreqDom = fft(I1 + 1i*Q1); -% Perform search for all listed PRN numbers ... -fprintf('('); -for PRN = settings.acqSatelliteList - %% Coarse acquisition =========================================== - - % Generate B1C data codes and sample them according to the sampling freq. - DataPriTable = makeDataTable(settings,PRN); - % generate local code duplicate to do correlate - localData = [DataPriTable(1:samplesXmsLen) ... - zeros(1,len10PlusXms - samplesXmsLen)]; - - % Perform DFT of B1C data code - DataPriFreqDom = conj(fft(localData)); - - % Pilot signal acquisition - if (settings.pilotACQflag == 1) - PilotPriTable = makePilotTable(settings,PRN); - localPilot = [PilotPriTable(1:samplesXmsLen) ... - zeros(1,len10PlusXms - samplesXmsLen)]; - PilotPriFreqDom = conj(fft(localPilot)); - end - %--- Make the correlation for whole frequency band (for all freq. bins) - for frqBinIndex = 1:numberOfFrqBins - IQfreqDomShifted = circshift(IQfreqDom, frqBinIndex - 1); - % Multiplication in the frequency domain (correlation in time domain) - convCodeIQ1 = IQfreqDomShifted .* DataPriFreqDom; - % Perform inverse DFT and store correlation results - results(frqBinIndex, :) = abs(ifft(convCodeIQ1)); - - % Pilot signal acquisition - if (settings.pilotACQflag == 1) - convCodeIQ1 = IQfreqDomShifted .* PilotPriFreqDom; - % Non-coherent combining of data and pilot results - results(frqBinIndex, :) = (results(frqBinIndex, :)* sqrt(11)+ ... - abs(ifft(convCodeIQ1))* sqrt(29) )/ sqrt(40); - end - - end % frqBinIndex = 1:numberOfFrqBins - - %% Look for correlation peaks in the results ==================== - % Find the highest peak and compare it to the second highest peak - % The second peak is chosen not closer than 1 chip to the highest peak - - % Find the correlation peak and the carrier frequency - [~, frequencyBinIndex] = max(max(results, [], 2)); - selFreq = initFreq + (frequencyBinIndex -1)*settings.acqStep; - % Find code phase of the same correlation peak - [peakSize, codePhase] = max(max(results)); - % Store GLRT statistic - acqResults.peakMetric(PRN) = peakSize/sigPower; - - % To prevent index from exceeding matrix dimensions in the fine - % qacquisition, move to previous code start position. - if (codePhase+samplesPerCode-1) > length(longSignal) - codePhase = codePhase - samplesPerCode; - end - - % If the result is above threshold, then there is a signal ... - if acqResults.peakMetric(PRN) > settings.acqThreshold - %% Fine resolution frequency search ========================= - % Indicate PRN number of the detected signal - fprintf('%02d ', PRN); - % Cut 10ms input signal - signal0DC = longSignal(codePhase:(codePhase + samplesPerCode-1)); - % Remove B1C code modulation from the original signal - xCarrier = signal0DC .* DataPriTable; - - % Pilot signal acquisition - if (settings.pilotACQflag == 1) - xCarrierPilot = signal0DC .* PilotPriTable; - end - - %--- Search different frequency bins ------------------------------ - for FineBinIndex = 1 : NumOfFineBins - % Carrier frequencies of the fine frequency bins - FineFrqBins(FineBinIndex) = selFreq - settings.acqStep + fineStep * (FineBinIndex - 1); - % Generate local sine and cosine - sigCarr10cm = exp(1i*FineFrqBins(FineBinIndex) * finePhasePoints); - FineResult(FineBinIndex) = abs(sum(xCarrier .* sigCarr10cm)); - - %--- Pilot signal acquisition - if (settings.pilotACQflag == 1) - FineResult(FineBinIndex) = (FineResult(FineBinIndex)*11 + ... - abs(sum(xCarrierPilot .* sigCarr10cm))*29)/40; - end - end - - % Find the fine carrier freq. ------------------------------------- - % Corresponding to the largest noncoherent power - [~,maxFinBin] = max(FineResult); - acqResults.carrFreq(PRN) = FineFrqBins(maxFinBin); - - %signal found, if IF =0 just change to 1 Hz to allow processing - if(acqResults.carrFreq(PRN) == 0) - acqResults.carrFreq(PRN) = 1; - end - acqResults.codePhase(PRN) = codePhase; - - %% Downsampling recovery =================================== - % Find acquisition results corresponding to orignal sampling freq - if (exist('oldFreq', 'var') && settings.resamplingflag == 1) - % Find code phase - acqResults.codePhase(PRN) = floor((codePhase - 1)/ ... - settings.samplingFreq * oldFreq)+1; - - % Doppler frequency - if (settings.IF >= settings.samplingFreq/2) - % In this condition, the FFT computed freq. is symmetric - % with the true frequemcy about half of the sampling - % frequency, so we have the following: - IF_temp = settings.samplingFreq - settings.IF; - doppler = IF_temp - acqResults.carrFreq(PRN); - else - doppler = acqResults.carrFreq(PRN) - settings.IF; - end - % Carrier freq. corresponding to orignal sampling freq - acqResults.carrFreq(PRN) = doppler + oldIF; - end - - else - %--- No signal with this PRN -------------------------------------- - fprintf('. '); - end % if (peakSize/secondPeakSize) > settings.acqThreshold - -end % for PRN = satelliteList - -%=== Acquisition is over ================================================== -fprintf(')\n'); diff --git a/BDS/B1C/include/postProcessing.m b/BDS/B1C/include/postProcessing.m index 0ad8863..7e84802 100644 --- a/BDS/B1C/include/postProcessing.m +++ b/BDS/B1C/include/postProcessing.m @@ -42,7 +42,8 @@ %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); % Normal acquisition - acqResults = acquisitionCS(data, settings); + acqResults = acquisition(data, settings); + save("acqResults") end %% Initialize channels and prepare for the run ============================ @@ -72,6 +73,7 @@ % Tracking with Matlab correlator [trkResults, ~] = WB_tracking(fid, channel, settings); end +save("trkResults") % Close the data file fclose(fid); @@ -82,10 +84,16 @@ %% Calculate navigation solutions ========================================= disp(' Calculating navigation solutions...'); [navResults, ~] = postNavigation(trkResults, settings); +save("navResults") disp(' Processing is complete for this data block'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults, settings); +end + if settings.plotTracking plotTracking(1:settings.numberOfChannels, trkResults, settings); end @@ -95,6 +103,7 @@ end disp('Post processing of the signal is over.'); + else % Error while opening the data file. error('Unable to read file %s: %s.', settings.fileName, message); diff --git a/BDS/B1C/init.m b/BDS/B1C/init.m index a138ad5..2636077 100644 --- a/BDS/B1C/init.m +++ b/BDS/B1C/init.m @@ -1,20 +1,16 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- -% CU Multi-GNSS SDR -% (C) Developed for BDS B1C SDR by Yafeng Li, Nagaraj C. Shivaramaiah -% and Dennis M. Akos. -% Based on the original framework for GPS C/A SDR by Darius Plausinaitis, -% Peter Rinder, Nicolaj Bertelsen and Dennis M. Akos -% -% Reference: Li, Y., Shivaramaiah, N.C. & Akos, D.M. Design and -% implementation of an open-source BDS-3 B1C/B2a SDR receiver. -% GPS Solut (2019) 23: 60. https://doi.org/10.1007/s10291-019-0853-z +% CU Multi-GNSS SDR +% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -34,15 +30,15 @@ % $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ %% Clean up the environment first ========================================= -% clear; close all; clc; +% clear all; close all; clc; format ('compact'); format ('long', 'g'); %--- Include folders with functions --------------------------------------- -addpath include % The software receiver functions -addpath Common % Common functions between differnt SDR receivers -% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers %% Print startup ========================================================== % fprintf(['\n',... @@ -58,15 +54,14 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try fprintf('Probing data (%s)...\n', settings.fileName) - probeData(settings); + probeData(settings); catch % There was an error, print it and exit - errStruct = lasterror; %#ok + errStruct = lasterror; disp(errStruct.message); disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') return; @@ -75,10 +70,10 @@ disp(' Raw IF data plotted ') disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') disp(' '); -gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit: '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); if (gnssStart == 1) disp(' '); - %start things rolling... + % start things rolling... postProcessing end diff --git a/BDS/B1C/initSettings.m b/BDS/B1C/initSettings.m index 2e94702..09150e7 100644 --- a/BDS/B1C/initSettings.m +++ b/BDS/B1C/initSettings.m @@ -43,6 +43,7 @@ % $Id: initSettings.m,v 1.9.2.31 2006/08/18 11:41:57 dpl Exp $ %% Raw signal file name and related parameter ======================= +settings.fileName = '/media/gnss/Ext2TB/data/L1/B200Ex_1.bin'; % Data type used to store one sample settings.dataType = 'schar'; % File Types @@ -51,7 +52,7 @@ settings.fileType = 2; %Intermediate, sampling, code freqency and code length settings.IF = -20e3; %[Hz] -settings.samplingFreq = 5e6; %[Hz] +settings.samplingFreq = 18e6; %[Hz] % Front-end bandwidth for calculation of weighting factor in wideband (WB) % tracking mode (this is required only for WB tracking mode) settings.FEBW = 27e6; % [Hz] @@ -62,7 +63,7 @@ settings.msToProcess = 40000; %[ms] % List of satellites to look for. Some satellites can be excluded to speed % up acquisition. In May 2020, BeiDou3 includes PRNs with 19:30,32:46,59,60. -settings.acqSatelliteList = [19:30,32:46,59,60]; +settings.acqSatelliteList = [1:69]; % Enable use of pilot channel to do acquisition settings.pilotACQflag = 1; % 0 - Off; 1 - On % Pilot channel tracking mode setting @@ -126,11 +127,11 @@ settings.truePosition.N = nan; settings.truePosition.U = nan; -%% Plot settings ==================================================== +%% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 0; % 0 - Off; 1 - On -settings.plotAcquisition = 0; -settings.plotNavigation = 0; +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ======================================================== settings.c = 299792458; % The speed of light, [m/s] diff --git a/BDS/B1I/include/makeCaTableDMA.m b/BDS/B1I/include/makeCaTableDMA.m new file mode 100644 index 0000000..76800d6 --- /dev/null +++ b/BDS/B1I/include/makeCaTableDMA.m @@ -0,0 +1,82 @@ +function caCodesTable = makeCaTableDMA2B(settings, Ncodes) +%Function generates CA codes for all 32 satellites based on the settings +%provided in the structure "settings". The codes are digitized at the +%sampling frequency specified in the settings structure. +%One row in the "caCodesTable" is one C/A code. The row number is the PRN +%number of the C/A code. +% +%caCodesTable = makeCaTable(settings) +% +% Inputs: +% settings - receiver settings +% Outputs: +% caCodesTable - an array of arrays (matrix) containing C/A codes +% for all satellite PRN-s + +%-------------------------------------------------------------------------- +% CU Multi-GNSS SDR +% (C) Updated by Daehee Won, Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos +%-------------------------------------------------------------------------- +%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: makeCaTable.m,v 1.1.2.6 2006/08/14 11:38:22 dpl Exp $ + +%-------------------------------------------------------------------------- +% Modified for Beidou by Daehee Won +% Final update: Sep. 26, 2013 + +%--- Find number of samples per spreading code ---------------------------- +samplesPerCode = round(settings.samplingFreq / ... + (settings.codeFreqBasis / settings.codeLength)); + +%--- Prepare the output matrix to speed up function ----------------------- +caCodesTable = zeros(58, samplesPerCode); % Beidou has 58 PRNs + +%--- Find time constants -------------------------------------------------- +ts = 1/settings.samplingFreq; % Sampling period in sec +tc = 1/settings.codeFreqBasis; % C/A chip period in sec + +samplesPerCode = round(settings.samplingFreq / ... + (settings.codeFreqBasis / (Ncodes*settings.codeLength))); +caCodesTable = zeros(58, samplesPerCode); % Beidou has 58 PRNs + +%=== For all satellite PRN-s ... +for PRN = 1:58 % Beidou has 58 PRNs + %--- Generate CA code for given PRN ----------------------------------- + caCode = generateCAcode53(PRN); + caCode = [caCode caCode]; + + %=== Digitizing ======================================================= + + %--- Make index array to read C/A code values ------------------------- + % The length of the index array depends on the sampling frequency - + % number of samples per millisecond (because one C/A code period is one + % millisecond). + codeValueIndex = ceil((ts * (1:samplesPerCode)) / tc); + + %--- Correct the last index (due to number rounding issues) ----------- + codeValueIndex(end) = (Ncodes*2046); % Beidou chip lenth is 2046 + + %--- Make the digitized version of the C/A code ----------------------- + % The "upsampled" code is made by selecting values form the CA code + % chip array (caCode) for the time instances of each sample. + caCodesTable(PRN, :) = caCode(codeValueIndex); + +end % for PRN = 1:58 diff --git a/BDS/B1I/include/postProcessing.m b/BDS/B1I/include/postProcessing.m index 53d85bc..27371f3 100644 --- a/BDS/B1I/include/postProcessing.m +++ b/BDS/B1I/include/postProcessing.m @@ -1,3 +1,57 @@ +% Script postProcessing.m processes the raw signal from the specified data +% file (in settings) operating on blocks of 37 seconds of data. +% +% First it runs acquisition code identifying the satellites in the file, +% then the code and carrier for each of the satellites are tracked, storing +% the 1msec accumulations. After processing all satellites in the 37 sec +% data block, then postNavigation is called. It calculates pseudoranges +% and attempts a position solutions. At the end plots are made for that +% block of data. +% +%-------------------------------------------------------------------------- +% CU Multi-GNSS SDR +% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos +%-------------------------------------------------------------------------- +%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. +%-------------------------------------------------------------------------- + +% THE SCRIPT "RECIPE" +% +% The purpose of this script is to combine all parts of the software +% receiver. +% +% 1.1) Open the data file for the processing and seek to desired point. +% +% 2.1) Acquire satellites +% +% 3.1) Initialize channels (preRun.m). +% 3.2) Pass the channel structure and the file identifier to the tracking +% function. It will read and process the data. The tracking results are +% stored in the trkResults structure. The results can be accessed this +% way (the results are stored each millisecond): +% trkResults(channelNumber).XXX(fromMillisecond : toMillisecond), where +% XXX is a field name of the result (e.g. I_P, codePhase etc.) +% +% 4) Pass tracking results to the navigation solution function. It will +% decode navigation messages, find satellite positions, measure +% pseudoranges and find receiver position. +% +% 5) Plot the results %% Initialization ========================================================= disp ('Starting processing...'); @@ -52,7 +106,7 @@ %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); acqResults = acquisition(data, settings); - save(acqPath, 'acqResults') + save('acqResults') end %% Initialize channels and prepare for the run ============================ @@ -75,6 +129,7 @@ % Process all channels for given data block [trkResults, ~] = tracking(fid, channel, settings); + save("trkResults") % Close the data file fclose(fid); @@ -82,10 +137,16 @@ %% Calculate navigation solutions ========================================= disp(' Calculating navigation solutions...'); [navResults, ~] = postNavigation(trkResults, settings); + save("navResults") disp(' Processing is complete for this data block'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults, settings); +end + if settings.plotTracking plotTracking(1:settings.numberOfChannels, trkResults, settings); end @@ -94,6 +155,7 @@ plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); + else % Error while opening the data file. error('Unable to read file %s: %s.', settings.fileName, message); diff --git a/BDS/B1I/init.m b/BDS/B1I/init.m index 3d76fcf..2636077 100644 --- a/BDS/B1I/init.m +++ b/BDS/B1I/init.m @@ -1,9 +1,11 @@ + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR -% (C) Updated by Daehee Won, Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -28,27 +30,27 @@ % $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ %% Clean up the environment first ========================================= -clear; close all; clc; +% clear all; close all; clc; format ('compact'); format ('long', 'g'); %--- Include folders with functions --------------------------------------- -addpath include % The software receiver functions -addpath ('Common') % Common functions between differnt SDR receivers -addpath ('IF_Data_Set') % IF data sets for each SDR receivers +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers %% Print startup ========================================================== -fprintf(['\n',... - 'Welcome to: softGNSS\n\n', ... - 'An open source GNSS SDR software project initiated by:\n\n', ... - ' Danish GPS Center/Aalborg University\n\n', ... - 'The code was improved by GNSS Laboratory/University of Colorado.\n\n',... - 'The software receiver softGNSS comes with ABSOLUTELY NO WARRANTY;\n',... - 'for details please read license details in the file license.txt. This\n',... - 'is free software, and you are welcome to redistribute it under\n',... - 'the terms described in the license.\n\n']); -fprintf(' -------------------------------\n\n'); +% fprintf(['\n',... +% 'Welcome to: softGNSS\n\n', ... +% 'An open source GNSS SDR software project initiated by:\n\n', ... +% ' Danish GPS Center/Aalborg University\n\n', ... +% 'The code was improved by GNSS Laboratory/University of Colorado.\n\n',... +% 'The software receiver softGNSS comes with ABSOLUTELY NO WARRANTY;\n',... +% 'for details please read license details in the file license.txt. This\n',... +% 'is free software, and you are welcome to redistribute it under\n',... +% 'the terms described in the license.\n\n']); +% fprintf(' -------------------------------\n\n'); %% Initialize constants, settings ========================================= settings = initSettings(); @@ -72,7 +74,6 @@ if (gnssStart == 1) disp(' '); - %start things rolling... + % start things rolling... postProcessing end - diff --git a/BDS/B1I/initSettings.m b/BDS/B1I/initSettings.m index f61c7d5..dbd8726 100644 --- a/BDS/B1I/initSettings.m +++ b/BDS/B1I/initSettings.m @@ -38,6 +38,8 @@ % $Id: initSettings.m,v 1.9.2.31 2006/08/18 11:41:57 dpl Exp $ %% Processing settings ==================================================== + +settings.fileName = '/media/gnss/Ext2TB/data/L1I/B200Ex_1.bin'; % Number of milliseconds to be processed used 36000 + any transients (see % below - in Nav parameters) to ensure nav subframes are provided settings.msToProcess = 61000; %61000; %[ms] @@ -66,8 +68,8 @@ settings.fileType = 2; % Intermediate, sampling and code frequencies -settings.IF = -44800; %-58000; % -40000, -44800 -settings.samplingFreq = 22e6; %53e6; %[Hz] +settings.IF = -20e3; %-58000; % -40000, -44800 +settings.samplingFreq = 18e6; %53e6; %[Hz] settings.codeFreqBasis = 2.046e6; %[Hz] Beidou has 2046 chip length. % Define number of chips in a code period @@ -122,9 +124,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 1; % 0 - Off - % 1 - On - +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] diff --git a/BDS/B2a/include/postProcessing.m b/BDS/B2a/include/postProcessing.m index 782adff..37dd1e6 100644 --- a/BDS/B2a/include/postProcessing.m +++ b/BDS/B2a/include/postProcessing.m @@ -96,6 +96,7 @@ %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); acqResults = acquisition(data, settings); + save("acqResults") end %% Initialize channels and prepare for the run ============================ @@ -119,6 +120,7 @@ % Process all channels for given data block [trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") % Close the data file fclose(fid); @@ -129,17 +131,23 @@ %% Calculate navigation solutions ========================================= disp(' Calculating navigation solutions...'); [navResults,~] = postNavigation(trkResults, settings); +save("navResults") disp(' Processing is complete for this data block'); disp('Post processing of the signal is over.'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking -plotTracking(1:settings.numberOfChannels, trkResults, settings); + plotTracking(1:settings.numberOfChannels, trkResults, settings); end if settings.plotNavigation -plotNavigation(navResults, settings); + plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); diff --git a/BDS/B2a/init.m b/BDS/B2a/init.m index a401554..2636077 100644 --- a/BDS/B2a/init.m +++ b/BDS/B2a/init.m @@ -1,10 +1,11 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -29,15 +30,15 @@ % $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ %% Clean up the environment first ========================================= -% clear; close all; clc; +% clear all; close all; clc; format ('compact'); format ('long', 'g'); %--- Include folders with functions --------------------------------------- -addpath include % The software receiver functions -addpath Common % Common functions between differnt SDR receivers -% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers %% Print startup ========================================================== % fprintf(['\n',... @@ -53,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -61,19 +61,19 @@ probeData(settings); catch % There was an error, print it and exit - errStruct = lasterror; %#ok + errStruct = lasterror; disp(errStruct.message); - disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') + disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit: '); -% -% if (gnssStart == 1) -% disp(' '); -% %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end diff --git a/BDS/B2a/initSettings.m b/BDS/B2a/initSettings.m index fa56e4b..2566bc4 100644 --- a/BDS/B2a/initSettings.m +++ b/BDS/B2a/initSettings.m @@ -51,7 +51,7 @@ %% 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 -settings.fileName = 'Beidou_B2a_IF_signal.bin'; +settings.fileName = '/media/gnss/Ext2TB/data/L5/B200Ex_3.bin'; % Data type used to store one sample settings.dataType = 'schar'; % File Types @@ -59,8 +59,8 @@ %2 - 8 bit I/Q samples I0,Q0,I1,Q1,I2,Q2,... settings.fileType = 2; % Intermediate, sampling -settings.IF = -0.05e6; % [Hz] -settings.samplingFreq = 20e6; % [Hz] +settings.IF = -20e3; % [Hz] +settings.samplingFreq = 18e6; % [Hz] %% Code parameter setting % Define number of chips in a code period and code frequencies of B2a settings.codeLength = 10230; % [chip] @@ -112,10 +112,11 @@ settings.truePosition.N = nan; settings.truePosition.U = nan; -%% Plot settings ==================================================== +%% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 1; % 0 - Off - % 1 - On +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ======================================================== settings.c = 299792458; % The speed of light, [m/s] settings.startOffset = 68.802; % [ms] Initial sign. travel time diff --git a/BDS/B3I/include/postProcessing.m b/BDS/B3I/include/postProcessing.m index c4fa8da..692cc34 100644 --- a/BDS/B3I/include/postProcessing.m +++ b/BDS/B3I/include/postProcessing.m @@ -97,7 +97,7 @@ %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); acqResults = acquisition(data, settings); - save(acqPath, "acqResults") + save("acqResults") end %% Initialize channels and prepare for the run ============================ @@ -121,6 +121,7 @@ % Process all channels for given data block [trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") % Close the data file fclose(fid); disp([' Tracking is over (elapsed time ', ... @@ -128,17 +129,24 @@ %% Calculate navigation solutions ========================================= disp(' Calculating navigation solutions...'); [navResults, ~] = postNavigation(trkResults, settings); +save("navResults") disp(' Processing is complete for this data block'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults, settings); +end + if settings.plotTracking -plotTracking(1:settings.numberOfChannels, trkResults, settings); + plotTracking(1:settings.numberOfChannels, trkResults, settings); end if settings.plotNavigation -plotNavigation(navResults, settings); + plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); + else % Error while opening the data file. error('Unable to read file %s: %s.', settings.fileName, message); diff --git a/BDS/B3I/init.m b/BDS/B3I/init.m index 6e4694e..2636077 100644 --- a/BDS/B3I/init.m +++ b/BDS/B3I/init.m @@ -1,11 +1,11 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- -% CU Multi-GNSS SDR -% (C) Developed for BDS B3I SDR by Yafeng Li, Nagaraj C. Shivaramaiah -% and Dennis M. Akos. -% Based on the original framework for GPS C/A SDR by Darius Plausinaitis, -% Peter Rinder, Nicolaj Bertelsen and Dennis M. Akos +% CU Multi-GNSS SDR +% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -30,15 +30,15 @@ % $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ %% Clean up the environment first ========================================= -% clear; close all; clc; +% clear all; close all; clc; format ('compact'); format ('long', 'g'); %--- Include folders with functions --------------------------------------- -addpath include % The software receiver functions -addpath Common % Common functions between differnt SDR receivers -% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers %% Print startup ========================================================== % fprintf(['\n',... @@ -54,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -68,13 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); -% %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end diff --git a/BDS/B3I/initSettings.m b/BDS/B3I/initSettings.m index 40e582e..1995999 100644 --- a/BDS/B3I/initSettings.m +++ b/BDS/B3I/initSettings.m @@ -38,6 +38,7 @@ % $Id: initSettings.m,v 1.9.2.31 2006/08/18 11:41:57 dpl Exp $ %% Processing settings ==================================================== +settings.fileName = '/media/gnss/Ext2TB/data/B3I/B200Ex_1.bin'; % Number of milliseconds to be processed used 36000 + any transients (see % below - in Nav parameters) to ensure necessary nav subframes are provided settings.msToProcess = 45000; %[ms] @@ -50,9 +51,9 @@ settings.skipNumberOfBytes = 0; %% Raw signal file name and other parameter =============================== +settings.fileName = '/media/gnss/Ext2TB/data/B3/B200Ex_1.bin'; % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = 'Beidou_B3I_IF_signal.bin'; % Data type used to store one sample settings.dataType = 'schar'; % File Types @@ -60,8 +61,8 @@ %2 - 8 bit I/Q samples I0,Q0,I1,Q1,I2,Q2,... settings.fileType = 2; % Intermediate and sampling frequencies -settings.IF = -0.52e6; % [Hz] -settings.samplingFreq = 10e6; % [Hz] +settings.IF = -20e3; % [Hz] +settings.samplingFreq = 18e6; % [Hz] %% Code parameter setting ================================================= % Define number of chips in a code period and code frequencies @@ -110,10 +111,11 @@ settings.truePosition.N = nan; settings.truePosition.U = nan; -%% Plot settings ==================================================== +%% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 1; % 0 - Off - % 1 - On +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ======================================================== % The speed of light settings.c = 299792458; %[m/s] diff --git a/GAL/GAL_E1C/include/postProcessing.m b/GAL/GAL_E1C/include/postProcessing.m index 08d31f8..01779f7 100644 --- a/GAL/GAL_E1C/include/postProcessing.m +++ b/GAL/GAL_E1C/include/postProcessing.m @@ -7,12 +7,12 @@ % data block, then postNavigation is called. It calculates pseudoranges % and attempts a position solutions. At the end plots are made for that % block of data. + %-------------------------------------------------------------------------- -% CU Multi-GNSS SDR -% (C) Developed for BDS B3I SDR by Yafeng Li, Nagaraj C. Shivaramaiah -% and Dennis M. Akos. -% Based on the original framework for GPS C/A SDR by Darius Plausinaitis, -% Peter Rinder, Nicolaj Bertelsen and Dennis M. Akos +% CU Multi-GNSS SDR +% (C) Updated by Jakob Almqvist, Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- %This program is free software; you can redistribute it and/or %modify it under the terms of the GNU General Public License @@ -52,6 +52,7 @@ % pseudoranges and find receiver position. % % 5) Plot the results. + %% Initialization ========================================================= disp ('Starting processing...'); @@ -82,8 +83,9 @@ samplesPerCode = round(settings.samplingFreq / ... (settings.codeFreqBasis / settings.codeLength)); - % At least 27ms of signal are needed for fine frequency estimation - codeLen = max(27,settings.acqNonCohTime+2); + + % At least 42ms of signal are needed for fine frequency estimation + codeLen = max(42,settings.acqNonCohTime+2); % Read data for acquisition. data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; @@ -96,7 +98,8 @@ %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); acqResults = acquisition(data, settings); - end + save("acqResults") +end %% Initialize channels and prepare for the run ============================ @@ -119,25 +122,33 @@ % Process all channels for given data block [trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") % Close the data file fclose(fid); - disp([' Tracking is over (elapsed time ', ... - datestr(now - startTime, 13), ')']) + datestr(now - startTime, 13), ')']) %% Calculate navigation solutions ========================================= disp(' Calculating navigation solutions...'); + [navResults, ~] = postNavigation(trkResults, settings); -disp(' Processing is complete for this data block'); +save("navResults") +disp(' Processing is complete for this data block'); +disp('Post processing of the signal is over.'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking -plotTracking(1:settings.numberOfChannels, trkResults, settings); + plotTracking(1:settings.numberOfChannels, trkResults, settings); end if settings.plotNavigation -plotNavigation(navResults, settings); + plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); @@ -145,3 +156,4 @@ % Error while opening the data file. error('Unable to read file %s: %s.', settings.fileName, message); end % if (fid > 0) + diff --git a/GAL/GAL_E1C/init.m b/GAL/GAL_E1C/init.m index 46552e0..2636077 100644 --- a/GAL/GAL_E1C/init.m +++ b/GAL/GAL_E1C/init.m @@ -1,10 +1,11 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -29,15 +30,15 @@ % $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ %% Clean up the environment first ========================================= -% clear; close all; clc; +% clear all; close all; clc; format ('compact'); format ('long', 'g'); %--- Include folders with functions --------------------------------------- -addpath include % The software receiver functions -addpath Common % Common functions between differnt SDR receivers -% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers %% Print startup ========================================================== % fprintf(['\n',... @@ -53,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -67,14 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); - %start things rolling... - postProcessing -% end +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); +if (gnssStart == 1) + disp(' '); + % start things rolling... + postProcessing +end diff --git a/GAL/GAL_E1C/initSettings.m b/GAL/GAL_E1C/initSettings.m index ff59d4f..b9a3b17 100644 --- a/GAL/GAL_E1C/initSettings.m +++ b/GAL/GAL_E1C/initSettings.m @@ -55,7 +55,7 @@ % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = '/media/gnss/Ext2TB/JMBF/L1/Lab16Nov_L1schar_Fs5e6_IF20e3.bin'; +settings.fileName = '/media/gnss/Ext2TB/data/L5/B200Ex_2.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -66,8 +66,8 @@ settings.fileType = 2; %Intermediate, sampling, code and L1 frequencies -settings.IF = 20e3; %[Hz] -settings.samplingFreq = 5e6; %[Hz] +settings.IF = -20e3; %[Hz] +settings.samplingFreq = 18e6; %[Hz] settings.codeFreqBasis = 1.023e6; % [Hz] % Define number of chips in a code period @@ -125,10 +125,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 0; % 0 - Off - % 1 - On -settings.plotNavigation = 0; -settings.plotAcquisition = 0; +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] diff --git a/GAL/GAL_E5a/include/NAVdecodingv2.m b/GAL/GAL_E5a/include/NAVdecodingv2.m new file mode 100644 index 0000000..d99a8c9 --- /dev/null +++ b/GAL/GAL_E5a/include/NAVdecodingv2.m @@ -0,0 +1,178 @@ +function [eph, firstSubFrame] = NAVdecodingv2(I_P,settings) + + % findPreambles finds the first preamble occurrence in the bit stream of + % each channel. The preamble is verified by check of the spacing between + % preambles (6sec) and parity checking of the first two words in a + % subframe. At the same time function returns list of channels, that are in + % tracking state and with valid preambles in the nav data stream. + % + %[eph, subFrameStart,SOW] = CNAVdecoding(I_P_InputBits) + % + % Inputs: + % I_P_InputBits - output from the tracking function + % + % Outputs: + % firstSubFrame - Starting positions of the first message in the + % input bit stream I_P_InputBits in each channel. + % The position is CNAV bit(20ms before convolutional decoding) + % count since start of tracking. Corresponding value will + % be set to inf if no valid preambles were detected in + % the channel. + % SOW - Time Of Week (SOW) of the first message(in seconds). + % Corresponding value will be set to inf if no valid preambles + % were detected in the channel. + % eph - SV ephemeris. + + %-------------------------------------------------------------------------- + % CU Multi-GNSS SDR + % (C) Written by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos + %-------------------------------------------------------------------------- + % + %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: findPreambles.m,v 1.1.2.10 2017/01/19 21:13:22 dpl Exp $ + + + %% Initialize ephemeris structute --------------------------------------- + % Preamble search can be delayed to a later point in the tracking results + % to avoid noise due to tracking loop transients + searchStartOffset = 0; + % TOW of the first message + TOW = inf; + %Creates a cyclic redundancy code (CRC) detector System object + % crcDet = comm.CRCDetector([24 23 18 17 14 11 10 7 6 5 4 3 1 0]); + crcDet = comm.CRCDetector([24 23 18 17 14 11 10 7 6 5 4 3 1 0]); + + % Convert CNAV-producing convolutional code polynomials to trellis description + % Note that the difference from GPS is that the second branch G2 is + % inverted at the end (see ICD) + trellis = poly2trellis(7,[171 ~133]); + % Viterbi traceback depth for vitdec(function) + tblen = 35; + % Secondary code is "842E9" + % (Consider each bit from the nav message is a full secondary code.) + secondCode = 1-2*[1 0 0 0 0 1 0 0 0 0 1 0 1 1 1 0 1 0 0 1]; + % Sync pattern + anti_sync_bits = [-1 1 -1 -1 1 -1 -1 -1 1 1 1 1]; + % Define eph structure and initialize validity flag to zero. + eph = eph_structure_init(); + firstSubFrame = 0; + eph.flag = 0; + TOWflg = 1; + wordCount = []; + %% Find preamble to decode page ------------------------------------------- + % "Upsample" the preamble - make 20 vales per one bit. The preamble must be + % found with precision of a sample. + % This does: + % [ [1 0 1 1 0 1 1 1 0 0 0 0]*-1; + % [1 0 1 1 0 1 1 1 0 0 0 0]* 1 ]; + preamble_ms = kron(anti_sync_bits, secondCode); + + % Use the prompt correlator as the symbol stream, skip start of record if + % set in initSettings to avoid tracking loop transients + symbols = I_P(1 + searchStartOffset : end); + + % Now threshold the output and convert it to -1 and +1 + symbols(symbols > 0) = 1; + symbols(symbols <= 0) = -1; + + % Correlate tracking output with the preamble + tlmXcorrResult = xcorr(symbols, preamble_ms); + xcorrLength = (length(tlmXcorrResult) + 1) /2; + + + % Find at what index/ms the preambles start + index = find(round(abs(tlmXcorrResult(xcorrLength : xcorrLength * 2 - 1))) >= 239.99)'; + % Check each index is 10e3 a part from at least another index + newIndex = []; + for idx=1:size(index,1) + if find(abs(index-index(idx))==10e3) + newIndex = [newIndex;index(idx)]; + end + end + index = newIndex; + if isempty(index) + fprintf(' Sync pattern not found in navigation data') + return + end + % Save the starting point + firstSubFrame = index(1); + %% Start decoding current subframe, page by page. ------------------------- + %Indicator the requisite messages are all decoded + wordType = 0; + validCounter = 0; + % Analyze detected preamble like patterns ================================ + for ind = 1:size(index) % For each occurrence + % Check there are enough bits. + if (settings.msToProcess - index(ind)) < 10e3 + fprintf(' Not enough symbols to demodulate navData!', ind) + break + end + + % Select a whole page from the current subframe (1000 symbols, or 500 nav bits) + page = symbols(index(ind):(index(ind) + 10e3 - 1)); + % Convert from symbols to navigation bits + navBits = (secondCode*reshape(page, 20, 500)) > 0; + %--- Correct polarity of the all data bits according to preamble bits + sync_bits = [1 0 1 1 0 1 1 1 0 0 0 0]; + if navBits(1:12)*sync_bits'~=6 + navBits = not(navBits); + end + % Remove preamble from data + navBits = navBits(13:end); + + %--- De-interleave page parts ----------------------------------------- + symMat = reshape(navBits,61,8)'; + pageSym = reshape(symMat,1,[])'; + + %--- Remove convolutional encoding from page parts -------------------- + decBits = vitdec(pageSym,trellis,tblen,'trunc','hard'); + + %--- Retrieve full page -------------------------------------------- + page = decBits(1:238)'; + + %--- Check the CRC ---------------------------------------------------- + [~,frmError] = step(crcDet,page'); + if (frmError) +%--- Decode the message type ------------------------------------------ + fprintf('\nChecksum error.\n') + return + end + + %% Retrieve ephemeris data from current page -------------------------- + navWord = dec2bin(page)'; + [eph, validWord, wordType]= ephemerisv2(navWord, eph); + % Correct TOW if could be retrieved + if TOWflg && ~isempty(eph.TOW) + % Correct TOW to time for first page part + TOW = eph.TOW - (ind-1)*10; + TOWflg = 0; + end + % Increase counter only if a different page has been read + if ~any(wordType == wordCount) + wordCount = [wordCount, wordType]; + validCounter = validCounter + validWord; + end + % Check if enough navigation pages have been decoded + if validCounter ==4 + break + end + end % + eph.TOW = TOW; + eph.flag = floor(validCounter/4); +end diff --git a/GAL/GAL_E5a/include/acquisition.m b/GAL/GAL_E5a/include/acquisition.m index a0a4800..50517c8 100644 --- a/GAL/GAL_E5a/include/acquisition.m +++ b/GAL/GAL_E5a/include/acquisition.m @@ -298,4 +298,4 @@ end % for PRN = satelliteList %=== Acquisition is over ================================================== -fprintf(')\n'); +fprintf(')\n'); \ No newline at end of file diff --git a/GAL/GAL_E5a/include/ephemerisv2.m b/GAL/GAL_E5a/include/ephemerisv2.m new file mode 100644 index 0000000..cb0dda0 --- /dev/null +++ b/GAL/GAL_E5a/include/ephemerisv2.m @@ -0,0 +1,96 @@ +function [eph, validWord, wordType] = ephemerisv2(navWord, eph) + +%-------------------------------------------------------------------------- +% CU Multi-GNSS SDR +% (C) Written by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +%-------------------------------------------------------------------------- + +%--- Convert navigation data word to binary --------------------------- + +%--- Decode the message type ------------------------------------------ +wordType = bin2dec(navWord(1:6)); + +% Galileo Pi definition +galPi = 3.1415926535898; +% +validWord = 0; +%--- Decode sub-frame based on the message type ----------------------- +switch wordType + case 1 %--- SVID , Clock correction, SISA, Ionospheric correction, + % BGD, GST, Signal health and Data validity status + eph.SVID = bin2dec(navWord(7:12)); % 6 bits + eph.IODnav1 = bin2dec(navWord(13:22)); % 10 bits + eph.t_oc = bin2dec(navWord(23:36)) * 60; % 14 bits in s + eph.a_f0 = twosComp2dec(navWord(37:67)) * 2^(-34); % 31 bit in s + eph.a_f1 = twosComp2dec(navWord(68:88)) * 2^(-46); % 21 bit in s/s + eph.a_f2 = twosComp2dec(navWord(89:94)) * 2^(-59); % 6 bit in s/(s^2) +% eph.SISA_E1_E5a = bin2dec(navWord(95:102)); % 8 bit eph.SISA = navWord(121:128) % content not currently defined + eph.a_i0 = bin2dec(navWord(103:113)) * 2^(-2); % 11 bit + eph.a_i1 = twosComp2dec(navWord(114:124)) * 2^(-8); % 11 bit + eph.a_i2 = twosComp2dec(navWord(125:138)) * 2^(-15); % 14 bit + eph.iono_SF1 = bin2dec(navWord(139)); % 1 bit + eph.iono_SF2 = bin2dec(navWord(140)); % 1 bit + eph.iono_SF3 = bin2dec(navWord(141)); % 1 bit + eph.iono_SF4 = bin2dec(navWord(142)); % 1 bit + eph.iono_SF5 = bin2dec(navWord(143)); % 1 bit + eph.BGD_E1E5a = twosComp2dec(navWord(144:153)) * 2^(-32); % 10 bit + eph.E5a_HS = bin2dec(navWord(154:155)); % 2 bit E5a signal Health Status (0-OK) + eph.WN = bin2dec(navWord(156:167)); % 12 bit + eph.TOW = bin2dec(navWord(168:187)); % 20 bit + eph.E5a_DVS = bin2dec(navWord(188)); % 1 bit E5a Data validity status (0-valid) + + validWord = 1; + + case 2 %--- Ephemeris (1/3) and GST ------------------------------------- + eph.IODnav2 = bin2dec(navWord(7:16)); % 10 bit + eph.M_0 = twosComp2dec(navWord(17:48)) * 2^(-31) * galPi; % 32 bit in rad + eph.OmegaDot = twosComp2dec(navWord(49:72)) * 2^(-43) * galPi; % 24 bit in rad + eph.e = bin2dec(navWord(73:104)) * 2^(-33); % 32 bit + eph.sqrtA = bin2dec(navWord(105:136)) * 2^(-19); % 32 bit in m^0.2 + eph.Omega_0 = twosComp2dec(navWord(137:168)) * 2^(-31) * galPi; % 32 bit in rad + eph.iDot = twosComp2dec(navWord(169:182)) * 2^(-43) * galPi; % 14 bit in rad + + validWord = 1; + + case 3 %--- Ephemeris (2/3) and GST ------------------------------- + eph.IODnav3 = bin2dec(navWord(7:16)); % 10 bit + eph.i_0 = twosComp2dec(navWord(17:48)) * 2^(-31) * galPi; % 32 bit in rad + eph.omega = twosComp2dec(navWord(49:80)) * 2^(-31) * galPi; % 32 bit in rad + eph.deltan = twosComp2dec(navWord(81:96)) * 2^(-43) * galPi; % 16 bit in rad + eph.CUC = twosComp2dec(navWord(97:112)) * 2^(-29); % 16 bit in rad + eph.CUS = twosComp2dec(navWord(113:128)) * 2^(-29); % 16 bit in rad + eph.CRC = twosComp2dec(navWord(129:144)) * 2^(-5); % 16 bit in m + eph.CRS = twosComp2dec(navWord(145:160)) * 2^(-5); % 16 bit in m + eph.t_oe = bin2dec(navWord(161:174)) * 60; % 14 bit in s + + validWord = 1; + + case 4 %--- Ephemeris (3/3), GST-UTC conversion, GST-GPS conversion and TO W + eph.IODnav4 = bin2dec(navWord(7:16)); % 10 bit + eph.CIC = twosComp2dec(navWord(17:32)) * 2^(-29); % 16 bit in rad + eph.CIS = twosComp2dec(navWord(33:48)) * 2^(-29); % 16 bit in rad + eph.A0 = twosComp2dec(navWord(49:80)) * 2^(-30); % 32 bit + eph.A1 = twosComp2dec(navWord(81:104)) * 2^(-50); % 24 bit + eph.delt_LS = twosComp2dec(navWord(105:112)); % 8 bit + eph.t_ot = bin2dec(navWord(113:120)) * 3600; % 8 bit + eph.WN_ot = bin2dec(navWord(121:128)); % 8 bit + eph.WN_LSF = bin2dec(navWord(129:136)); % 8 bit + eph.DN = bin2dec(navWord(137:139)); % 3 bit + eph.delt_LSF = twosComp2dec(navWord(140:147)); % 8 bit + eph.t_og = bin2dec(navWord(148:155)) * 3600; % 8 bit + eph.A0_G = twosComp2dec(navWord(156:171)) * 2^(-35); % 16 bit + eph.A1_G = twosComp2dec(navWord(172:183)) * 2^(-51); % 12 bit + eph.WN_og = bin2dec(navWord(184:189)); % 6 bit + + validWord = 1; + + %case 5 %--- Almanac ------------------------------------- + + %case 6 %--- Almanac ------------------------------------- + + +end % switch word type + +end % all pages + + diff --git a/GAL/GAL_E5a/include/postProcessing.m b/GAL/GAL_E5a/include/postProcessing.m index 81e5d53..e4c2cff 100644 --- a/GAL/GAL_E5a/include/postProcessing.m +++ b/GAL/GAL_E5a/include/postProcessing.m @@ -7,125 +7,12 @@ % data block, then postNavigation is called. It calculates pseudoranges % and attempts a position solutions. At the end plots are made for that % block of data. -%-------------------------------------------------------------------------- -% CU Multi-GNSS SDR -% (C) Developed for BDS B3I SDR by Yafeng Li, Nagaraj C. Shivaramaiah -% and Dennis M. Akos. -% Based on the original framework for GPS C/A SDR by Darius Plausinaitis, -% Peter Rinder, Nicolaj Bertelsen and Dennis M. Akos -%-------------------------------------------------------------------------- -%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. -%-------------------------------------------------------------------------- -% THE SCRIPT "RECIPE" -% -% The purpose of this script is to combine all parts of the software -% receiver. -% -% 1.1) Open the data file for the processing and seek to desired point. -% -% 2.1) Acquire satellites -% -% 3.1) Initialize channels (preRun.m). -% 3.2) Pass the channel structure and the file identifier to the tracking -% function. It will read and process the data. The tracking results are -% stored in the trkResults structure. The results can be accessed this -% way (the results are stored each millisecond): -% trkResults(channelNumber).XXX(fromMillisecond : toMillisecond), where -% XXX is a field name of the result (e.g. I_P, codePhase etc.) -% -% 4) Pass tracking results to the navigation solution function. It will -% decode navigation messages, find satellite positions, measure -% pseudoranges and find receiver position. -% -% 5) Plot the results. - -%% Initialization ========================================================= -disp ('Starting processing...'); - -[fid, message] = fopen(settings.fileName, 'rb'); - -%Initialize the multiplier to adjust for the data type -if (settings.fileType==1) -dataAdaptCoeff=1; -else -dataAdaptCoeff=2; -end - -%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, dataAdaptCoeff*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) || ~exist('acqResults', 'var')) - - % Find number of samples per spreading code - samplesPerCode = round(settings.samplingFreq / ... - (settings.codeFreqBasis / settings.codeLength)); - - % Read data for acquisition. At least 102ms of signal are needed - % for the fine frequency estimation - codeLen = max(102,settings.acqNonCohTime+2); - % are needed for the fine frequency estimation - data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; - - if (dataAdaptCoeff == 2) - data1 = data(1:2:end); - data2 = data(2:2:end); - data = data1 + 1i .* data2; - end - - %--- Do the acquisition ------------------------------------------- - disp (' Acquiring satellites...'); - acqResults = acquisition(data, settings); -end -%% 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.'); - trkResults = []; - navResults = []; - return;% Script postProcessing.m processes the raw signal from the specified data -% file (in settings) operating on blocks of 37 seconds of data. -% -% First it runs acquisition code identifying the satellites in the file, -% then the code and carrier for each of the satellites are tracked, storing -% the 1msec accumulations. After processing all satellites in the 37 sec -% data block, then postNavigation is called. It calculates pseudoranges -% and attempts a position solutions. At the end plots are made for that -% block of data. %-------------------------------------------------------------------------- -% CU Multi-GNSS SDR -% (C) Developed for BDS B3I SDR by Yafeng Li, Nagaraj C. Shivaramaiah -% and Dennis M. Akos. -% Based on the original framework for GPS C/A SDR by Darius Plausinaitis, -% Peter Rinder, Nicolaj Bertelsen and Dennis M. Akos +% CU Multi-GNSS SDR +% (C) Updated by Jakob Almqvist, Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- %This program is free software; you can redistribute it and/or %modify it under the terms of the GNU General Public License @@ -196,22 +83,24 @@ samplesPerCode = round(settings.samplingFreq / ... (settings.codeFreqBasis / settings.codeLength)); - % Read data for acquisition. At least 102ms of signal are needed - % for the fine frequency estimation + + % At least 42ms of signal are needed for fine frequency estimation codeLen = max(102,settings.acqNonCohTime+2); - % are needed for the fine frequency estimation + % Read data for acquisition. data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; - if (dataAdaptCoeff == 2) - data1 = data(1:2:end); - data2 = data(2:2:end); - data = data1 + 1i .* data2; + if (dataAdaptCoeff==2) + data1=data(1:2:end); + data2=data(2:2:end); + data=data1 + 1i .* data2; end %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); acqResults = acquisition(data, settings); + save("acqResults") end + %% Initialize channels and prepare for the run ============================ % Start further processing only if a GNSS signal was acquired (the @@ -233,57 +122,27 @@ % Process all channels for given data block [trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") % Close the data file fclose(fid); - disp([' Tracking is over (elapsed time ', ... - datestr(now - startTime, 13), ')']) - + datestr(now - startTime, 13), ')']) %% Calculate navigation solutions ========================================= disp(' Calculating navigation solutions...'); + [navResults, ~] = postNavigation(trkResults, settings); -disp(' Processing is complete for this data block'); +save("navResults") +disp(' Processing is complete for this data block'); +disp('Post processing of the signal is over.'); %% Plot all results =================================================== disp (' Ploting results...'); -if settings.plotTracking - plotTracking(1:settings.numberOfChannels, trkResults, settings); -end - -if settings.plotNavigation - plotNavigation(navResults, settings); -end -disp('Post processing of the signal is over.'); -else -% Error while opening the data file. -error('Unable to read file %s: %s.', settings.fileName, message); - -end % if (fid > 0) - +if settings.plotAcquisition + plotAcquisition(acqResults); end -%% Track the signal ======================================================= -startTime = now; -disp ([' Tracking started at ', datestr(startTime)]); - -% Process all channels for given data block -[trkResults, ~] = tracking(fid, channel, settings); -% Close the data file -fclose(fid); - -disp([' Tracking is over (elapsed time ', ... - datestr(now - startTime, 13), ')']) - - -%% Calculate navigation solutions ========================================= -disp(' Calculating navigation solutions...'); -[navResults, ~] = postNavigation(trkResults, settings); -disp(' Processing is complete for this data block'); - -%% Plot all results =================================================== -disp (' Ploting results...'); if settings.plotTracking plotTracking(1:settings.numberOfChannels, trkResults, settings); end @@ -292,8 +151,9 @@ plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); + else % Error while opening the data file. error('Unable to read file %s: %s.', settings.fileName, message); - end % if (fid > 0) + diff --git a/GAL/GAL_E5a/init.m b/GAL/GAL_E5a/init.m index fa8ab69..2636077 100644 --- a/GAL/GAL_E5a/init.m +++ b/GAL/GAL_E5a/init.m @@ -1,10 +1,11 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -53,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -67,13 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); - %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end diff --git a/GAL/GAL_E5a/initSettings.m b/GAL/GAL_E5a/initSettings.m index e540ece..28d5cee 100644 --- a/GAL/GAL_E5a/initSettings.m +++ b/GAL/GAL_E5a/initSettings.m @@ -18,7 +18,7 @@ % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = 'dataL5.bin'; +settings.fileName = '/media/gnss/Ext2TB/data/L5/B200Ex_1.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -29,8 +29,8 @@ settings.fileType = 2; % Intermediate, sampling and code frequencies -settings.IF = -0.05e6; % [Hz]4.1304e6 -settings.samplingFreq = 20e6; % [Hz]16.3676e6 +settings.IF = -20e3; % [Hz]4.1304e6 +settings.samplingFreq = 18e6; settings.codeFreqBasis = 10.23e6; % [Hz] % Define number of chips in a E5aI/E5aQ tiered code period settings.codeLength = 10230; @@ -101,4 +101,4 @@ % Accumulation interval for computing VSM C/No (in ms) settings.CNo.VSMinterval = 100;%20; %% E5b carrier frequency =================================================== -settings.carrFreqBasis = 1176.45e6; %[Hz] +settings.carrFreqBasis = 1176.45e6; %[Hz] \ No newline at end of file diff --git a/GAL/GAL_E5b/include/NAVdecoding.m b/GAL/GAL_E5b/include/NAVdecoding.m index 7a38c88..8097f3d 100644 --- a/GAL/GAL_E5b/include/NAVdecoding.m +++ b/GAL/GAL_E5b/include/NAVdecoding.m @@ -106,7 +106,6 @@ % Find at what index/ms the preambles start index = find(abs(tlmXcorrResult(xcorrLength : xcorrLength * 2 - 1)) > 39.99)'; if isempty(index) - eph = []; return else % Analyze detected preamble like patterns ================================ diff --git a/GAL/GAL_E5b/include/NAVdecodingv2.m b/GAL/GAL_E5b/include/NAVdecodingv2.m index 8a61f65..b6f5729 100644 --- a/GAL/GAL_E5b/include/NAVdecodingv2.m +++ b/GAL/GAL_E5b/include/NAVdecodingv2.m @@ -219,7 +219,7 @@ end %--- Update TOW when read - if ~isempty(eph.TOW) + if isempty(eph.TOW) % Correct TOW to time for first page part TOW = eph.TOW - (pagecnt*2); % Reset parameter in case it has to be re-read @@ -234,7 +234,7 @@ end % if CRC is OK ... else % Not enough symbols to decode this page. - print('Not enough symbols to decode current page.') + disp('Not enough symbols to decode current page.') end % If enough symbols to decode page end % For loop iterating sync patterns diff --git a/GAL/GAL_E5b/include/acquisition.m b/GAL/GAL_E5b/include/acquisition.m index ecb2b9c..cd1b778 100644 --- a/GAL/GAL_E5b/include/acquisition.m +++ b/GAL/GAL_E5b/include/acquisition.m @@ -248,4 +248,4 @@ end % for PRN = satelliteList %=== Acquisition is over ================================================== -fprintf(')\n'); +fprintf(')\n'); \ No newline at end of file diff --git a/GAL/GAL_E5b/include/postNavigation.m b/GAL/GAL_E5b/include/postNavigation.m index 04d75a1..1b5f7f8 100644 --- a/GAL/GAL_E5b/include/postNavigation.m +++ b/GAL/GAL_E5b/include/postNavigation.m @@ -73,16 +73,13 @@ % [eph(PRN),subFrameStart(channelNr),TOW(channelNr)] = ... % NAVdecoding(trackResults(channelNr).I_P); % - [teph,subFrameStart(channelNr),TOW(channelNr)] = ... + [eph(PRN),subFrameStart(channelNr),TOW(channelNr)] = ... NAVdecoding(trackResults(channelNr).I_P); - cond1 = ~isempty(teph); + cond1 = ~isempty(eph(PRN)); cond2 = 0; cond3 = 0; if cond1 - cond2 = teph.flag == 1; - if cond2 - eph(PRN) = teph; - end + cond2 = eph(PRN).flag == 1; cond3 = eph(PRN).E5b_HS ==0; end %--- Exclude satellite if it does not have the necessary nav data ----- diff --git a/GAL/GAL_E5b/include/postProcessing.m b/GAL/GAL_E5b/include/postProcessing.m index 8d0358d..f29d9a7 100644 --- a/GAL/GAL_E5b/include/postProcessing.m +++ b/GAL/GAL_E5b/include/postProcessing.m @@ -53,6 +53,7 @@ % % 5) Plot the results. + %% Initialization ========================================================= disp ('Starting processing...'); @@ -60,80 +61,95 @@ %Initialize the multiplier to adjust for the data type if (settings.fileType==1) -dataAdaptCoeff=1; + dataAdaptCoeff=1; else -dataAdaptCoeff=2; + dataAdaptCoeff=2; end %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, dataAdaptCoeff*settings.skipNumberOfBytes, 'bof'); + + % 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, dataAdaptCoeff*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) || ~exist('acqResults', 'var')) + % Do acquisition if it is not disabled in settings or if the variable + % acqResults does not exist. + if ((settings.skipAcquisition == 0) || ~exist('acqResults', 'var')) + + % Find number of samples per spreading code + samplesPerCode = round(settings.samplingFreq / ... + (settings.codeFreqBasis / settings.codeLength)); + + % Read data for acquisition. At least 102ms of signal are needed + % for the fine frequency estimation + codeLen = max(102,settings.acqNonCohTime+2); + % are needed for the fine frequency estimation + data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; - % Find number of samples per spreading code - samplesPerCode = round(settings.samplingFreq / ... - (settings.codeFreqBasis / settings.codeLength)); - - % Read data for acquisition. At least 102ms of signal are needed - % for the fine frequency estimation - codeLen = max(102,settings.acqNonCohTime+2); - % are needed for the fine frequency estimation - data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; - - if (dataAdaptCoeff == 2) - data1 = data(1:2:end); - data2 = data(2:2:end); - data = data1 + 1i .* data2; - end + if (dataAdaptCoeff == 2) + data1 = data(1:2:end); + data2 = data(2:2:end); + data = data1 + 1i .* data2; + end - %--- Do the acquisition ------------------------------------------- - disp (' Acquiring satellites...'); - acqResults = acquisition(data, settings); -end + %--- Do the acquisition ------------------------------------------- + disp (' Acquiring satellites...'); + acqResults = acquisition(data, settings); + + save('acqResults') + end %% 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.'); - trkResults = []; - navResults = []; - return; -end + % 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.'); + trkResults = []; + navResults = []; + return; + end %% Track the signal ======================================================= -startTime = now; -disp ([' Tracking started at ', datestr(startTime)]); + startTime = now; + disp ([' Tracking started at ', datestr(startTime)]); + + % Process all channels for given data block + [trkResults, ~] = tracking(fid, channel, settings); + %plotTracking(1:settings.numberOfChannels, trkResults, settings, trkPath); + save( 'trkResults') + % Close the data file + fclose(fid); + + disp([' Tracking is over (elapsed time ', ... + datestr(now - startTime, 13), ')']) -% Process all channels for given data block -[trkResults, ~] = tracking(fid, channel, settings); -fclose(fid); + -disp([' Tracking is over (elapsed time ', ... - datestr(now - startTime, 13), ')']) +%% Calculate navigation solutions ========================================= + disp(' Calculating navigation solutions...'); + [navResults, ~] = postNavigation(trkResults, settings); + %plotNavigation(navResults, settings, navPath) + save('navResults') - -%% Calculate navigation solutions ========================================= -disp(' Calculating navigation solutions...'); -[navResults, ~] = postNavigation(trkResults, settings); -disp('Post processing of the signal is over.'); -%% Plot all results =================================================== + disp('Post processing of the signal is over.'); + + %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking plotTracking(1:settings.numberOfChannels, trkResults, settings); end @@ -142,8 +158,8 @@ plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); + else -% Error while opening the data file. -error('Unable to read file %s: %s.', settings.fileName, message); + % Error while opening the data file. + error('Unable to read file %s: %s.', settings.fileName, message); end % if (fid > 0) - diff --git a/GAL/GAL_E5b/init.m b/GAL/GAL_E5b/init.m index fa8ab69..f2fe6cf 100644 --- a/GAL/GAL_E5b/init.m +++ b/GAL/GAL_E5b/init.m @@ -1,10 +1,11 @@ -function [navSolutions] = init(filename) +1 %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -53,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -67,13 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); - %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end diff --git a/GAL/GAL_E5b/initSettings.m b/GAL/GAL_E5b/initSettings.m index dc7169b..6e33632 100644 --- a/GAL/GAL_E5b/initSettings.m +++ b/GAL/GAL_E5b/initSettings.m @@ -56,7 +56,7 @@ % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = 'data.bin'; +settings.fileName = '/media/gnss/Ext2TB/data/L5b/B200Ex_2.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -67,8 +67,8 @@ settings.fileType = 2; % Intermediate, sampling and code frequencies -settings.IF = -0.14e6; % [Hz]4.1304e6 -settings.samplingFreq = 20e6; % [Hz]16.3676e6 +settings.IF = -20e3; % [Hz]4.1304e6 +settings.samplingFreq = 18e6; % [Hz]16.3676e6 settings.codeFreqBasis = 10.23e6; % [Hz] % Define number of chips in a code period @@ -80,7 +80,7 @@ % List of satellites to look for. Some satellites can be excluded to speed % up acquisition. As of June 2020, in-orbit Galileo SVs includes PRNs: % [1 2 3 4 5 7 8 9 11 12 13 14 15 18 19 21 22 24 25 26 27 30 31 33 36] -settings.acqSatelliteList = [1 2 3 4 5 7 8 9 11 12 13 14 15 18 19 21 22 24 25 26 27 30 31 33 36]; +settings.acqSatelliteList = [1:36]; % Band around IF to search for satellite signal. Depends on max Doppler. % It is single sideband, so the whole search band is tiwce of it. settings.acqSearchBand = 5000; % [Hz] @@ -127,8 +127,9 @@ settings.truePosition.U = nan; %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 1; % 0 - Off - % 1 - On +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] diff --git a/GLO/GLO_GL1/acquisition.m b/GLO/GLO_GL1/include/acquisition.m similarity index 100% rename from GLO/GLO_GL1/acquisition.m rename to GLO/GLO_GL1/include/acquisition.m diff --git a/GLO/GLO_GL1/include/plotNavigation.m b/GLO/GLO_GL1/include/plotNavigation.m index f42aece..b0c6ae3 100644 --- a/GLO/GLO_GL1/include/plotNavigation.m +++ b/GLO/GLO_GL1/include/plotNavigation.m @@ -129,7 +129,7 @@ function plotNavigation(navSolutions, settings) skyPlot(handles(3, 2), ... navSolutions.az, ... navSolutions.el, ... - navSolutions.K(:, 1)); + navSolutions.PRN(:, 1)); title (handles(3, 2), ['Sky plot (mean PDOP: ', ... num2str(mean(navSolutions.DOP(2,:))), ')']); diff --git a/GLO/GLO_GL1/include/plotTracking.m b/GLO/GLO_GL1/include/plotTracking.m index 5e99c73..1e7abd8 100644 --- a/GLO/GLO_GL1/include/plotTracking.m +++ b/GLO/GLO_GL1/include/plotTracking.m @@ -50,7 +50,7 @@ function plotTracking(channelList, trackResults, settings) clf(channelNr +200); set(channelNr +200, 'Name', ['Channel ', num2str(channelNr), ... ' (K ', ... - num2str(trackResults(channelNr).K), ... + num2str(trackResults(channelNr).PRN), ... ') results']); %% Draw axes============================================================ @@ -154,7 +154,7 @@ function plotTracking(channelList, trackResults, settings) clf(channelNr +300); set(channelNr +300, 'Name', ['Channel ', num2str(channelNr), ... ' (K ', ... - num2str(trackResults(channelNr).K), ... + num2str(trackResults(channelNr).PRN), ... ') CNo']); plot(trackResults(channelNr).CNo.VSMValue);hold on plot(trackResults(channelNr).CNo.VSMValue,'mo');hold off diff --git a/GLO/GLO_GL1/postNavigation.m b/GLO/GLO_GL1/include/postNavigation.m similarity index 100% rename from GLO/GLO_GL1/postNavigation.m rename to GLO/GLO_GL1/include/postNavigation.m diff --git a/GLO/GLO_GL1/postProcessing.m b/GLO/GLO_GL1/include/postProcessing.m similarity index 93% rename from GLO/GLO_GL1/postProcessing.m rename to GLO/GLO_GL1/include/postProcessing.m index 36937e4..9cd0bb8 100644 --- a/GLO/GLO_GL1/postProcessing.m +++ b/GLO/GLO_GL1/include/postProcessing.m @@ -131,15 +131,22 @@ disp(' Processing is complete for this data block'); %% Plot all results =================================================== + disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking -plotTracking(1:settings.numberOfChannels, trkResults, settings); + plotTracking(1:settings.numberOfChannels, trkResults, settings); end if settings.plotNavigation -plotNavigation(navResults, settings); + plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); +disp('Post processing of the signal is over.'); else % Error while opening the data file. diff --git a/GLO/GLO_GL1/tracking.m b/GLO/GLO_GL1/include/tracking.m similarity index 100% rename from GLO/GLO_GL1/tracking.m rename to GLO/GLO_GL1/include/tracking.m diff --git a/GLO/GLO_GL1/init.m b/GLO/GLO_GL1/init.m new file mode 100644 index 0000000..2636077 --- /dev/null +++ b/GLO/GLO_GL1/init.m @@ -0,0 +1,79 @@ + +%-------------------------------------------------------------------------- +% CU Multi-GNSS SDR +% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos +%-------------------------------------------------------------------------- + +%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. +%-------------------------------------------------------------------------- +% +%Script initializes settings and environment of the software receiver. +%Then the processing is started. + +%-------------------------------------------------------------------------- +% CVS record: +% $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ + +%% Clean up the environment first ========================================= +% clear all; close all; clc; + +format ('compact'); +format ('long', 'g'); + +%--- Include folders with functions --------------------------------------- +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers + +%% Print startup ========================================================== +% fprintf(['\n',... +% 'Welcome to: softGNSS\n\n', ... +% 'An open source GNSS SDR software project initiated by:\n\n', ... +% ' Danish GPS Center/Aalborg University\n\n', ... +% 'The code was improved by GNSS Laboratory/University of Colorado.\n\n',... +% 'The software receiver softGNSS comes with ABSOLUTELY NO WARRANTY;\n',... +% 'for details please read license details in the file license.txt. This\n',... +% 'is free software, and you are welcome to redistribute it under\n',... +% 'the terms described in the license.\n\n']); +% fprintf(' -------------------------------\n\n'); + +%% Initialize constants, settings ========================================= +settings = initSettings(); + +%% Generate plot of raw data and ask if ready to start processing ========= +try + fprintf('Probing data (%s)...\n', settings.fileName) + probeData(settings); +catch + % There was an error, print it and exit + errStruct = lasterror; + disp(errStruct.message); + disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') + return; +end + +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... + postProcessing +end diff --git a/GLO/GLO_GL1/initSettings.m b/GLO/GLO_GL1/initSettings.m index 7a05d4d..6a39151 100644 --- a/GLO/GLO_GL1/initSettings.m +++ b/GLO/GLO_GL1/initSettings.m @@ -58,7 +58,7 @@ % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = 'data.bin'; +settings.fileName = '/home/gnss/Downloads/GL1_IF20KHz_FS12MHz/B200Ex_1.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -75,7 +75,7 @@ % Intermediate, sampling and code frequencies settings.IF = 0; %[Hz]1602e6-1590e6 -settings.samplingFreq = 10e6; %[Hz] +settings.samplingFreq = 12e6; %[Hz] settings.codeFreqBasis = 0.511e6; %[Hz] % Define number of chips in a code period @@ -131,8 +131,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 1; % 0 - Off - % 1 - On +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] settings.startOffset = 65; %[ms] Initial sign. travel time diff --git a/GLO/GLO_GL2/tracking.m b/GLO/GLO_GL2/Common/tracking.m similarity index 100% rename from GLO/GLO_GL2/tracking.m rename to GLO/GLO_GL2/Common/tracking.m diff --git a/GLO/GLO_GL2/acquisition.m b/GLO/GLO_GL2/include/acquisition.m similarity index 100% rename from GLO/GLO_GL2/acquisition.m rename to GLO/GLO_GL2/include/acquisition.m diff --git a/GLO/GLO_GL2/include/plotNavigation.m b/GLO/GLO_GL2/include/plotNavigation.m index f42aece..b0c6ae3 100644 --- a/GLO/GLO_GL2/include/plotNavigation.m +++ b/GLO/GLO_GL2/include/plotNavigation.m @@ -129,7 +129,7 @@ function plotNavigation(navSolutions, settings) skyPlot(handles(3, 2), ... navSolutions.az, ... navSolutions.el, ... - navSolutions.K(:, 1)); + navSolutions.PRN(:, 1)); title (handles(3, 2), ['Sky plot (mean PDOP: ', ... num2str(mean(navSolutions.DOP(2,:))), ')']); diff --git a/GLO/GLO_GL2/include/plotTracking.m b/GLO/GLO_GL2/include/plotTracking.m index 5e99c73..1e7abd8 100644 --- a/GLO/GLO_GL2/include/plotTracking.m +++ b/GLO/GLO_GL2/include/plotTracking.m @@ -50,7 +50,7 @@ function plotTracking(channelList, trackResults, settings) clf(channelNr +200); set(channelNr +200, 'Name', ['Channel ', num2str(channelNr), ... ' (K ', ... - num2str(trackResults(channelNr).K), ... + num2str(trackResults(channelNr).PRN), ... ') results']); %% Draw axes============================================================ @@ -154,7 +154,7 @@ function plotTracking(channelList, trackResults, settings) clf(channelNr +300); set(channelNr +300, 'Name', ['Channel ', num2str(channelNr), ... ' (K ', ... - num2str(trackResults(channelNr).K), ... + num2str(trackResults(channelNr).PRN), ... ') CNo']); plot(trackResults(channelNr).CNo.VSMValue);hold on plot(trackResults(channelNr).CNo.VSMValue,'mo');hold off diff --git a/GLO/GLO_GL2/postNavigation.m b/GLO/GLO_GL2/include/postNavigation.m similarity index 100% rename from GLO/GLO_GL2/postNavigation.m rename to GLO/GLO_GL2/include/postNavigation.m diff --git a/GLO/GLO_GL2/postProcessing.m b/GLO/GLO_GL2/include/postProcessing.m similarity index 96% rename from GLO/GLO_GL2/postProcessing.m rename to GLO/GLO_GL2/include/postProcessing.m index 45e6e56..aaf39ca 100644 --- a/GLO/GLO_GL2/postProcessing.m +++ b/GLO/GLO_GL2/include/postProcessing.m @@ -134,6 +134,10 @@ disp('Post processing of the signal is over.'); %% Plot all results =================================================== disp (' Ploting results...'); +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking plotTracking(1:settings.numberOfChannels, trkResults, settings); end diff --git a/GLO/GLO_GL2/include/tracking.m b/GLO/GLO_GL2/include/tracking.m new file mode 100644 index 0000000..202b71b --- /dev/null +++ b/GLO/GLO_GL2/include/tracking.m @@ -0,0 +1,351 @@ +function [trackResults, channel]= tracking(fid, channel, settings) +% Performs code and carrier tracking for all channels. +% +%[trackResults, channel] = tracking(fid, channel, settings) +% +% Inputs: +% fid - file identifier of the signal record for I +% channel - frequency channel, 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 spreading +% code's starting positions, together with other +% observation data from the tracking loops. All are +% saved every millisecond. + +%-------------------------------------------------------------------------- +% CU Multi-GNSS SDR +% (C) Updated by Jakob Almqvist, Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos +%-------------------------------------------------------------------------- +%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.31 2006/08/14 11:38:22 dpl Exp $ + +%% Initialize result structure ============================================ + +% 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 PRN 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); + +% Remain code and carrier phase +trackResults.remCodePhase = inf(1, settings.msToProcess); +trackResults.remCarrPhase = inf(1, settings.msToProcess); + +%C/No +trackResults.CNo.VSMValue = ... + zeros(1,floor(settings.msToProcess/settings.CNo.VSMinterval)); +trackResults.CNo.VSMIndex = ... + zeros(1,floor(settings.msToProcess/settings.CNo.VSMinterval)); + +%--- Copy initial settings for all channels ------------------------------- +trackResults = repmat(trackResults, 1, settings.numberOfChannels); + +% Get a vector with the C/A code sampled 1x/chip +caCode = generateCAcode(0,settings.codeFreqBasis,511); +% Then make it possible to do early and late versions +caCode = [caCode(511) caCode caCode(1)]; + +%% Initialize tracking variables ========================================== + +codePeriods = settings.msToProcess; % For GLONASS one C/A code is one ms + +%--- DLL variables -------------------------------------------------------- +% Define early-late offset (in chips) +earlyLateSpc = settings.dllCorrelatorSpacing; + +% Summation interval +PDIcode = settings.intTime; + +% Calculate filter coefficient values +[tau1code, tau2code] = calcLoopCoef(settings.dllNoiseBandwidth, ... + settings.dllDampingRatio, ... + 1.0); + +%--- PLL variables -------------------------------------------------------- +% Calculate filter coefficient values +[pf3,pf2,pf1] = calcLoopCoefCarr(settings); + +% Start waitbar +hwb = waitbar(0,'Tracking...','Visible','off'); + +%Adjust the size of the waitbar to insert text +CNoPos=get(hwb,'Position'); +set(hwb,'Position',[CNoPos(1),CNoPos(2),CNoPos(3),90],'Visible','on'); + + + +if (settings.fileType==1) + dataAdaptCoeff=1; +else + dataAdaptCoeff=2; +end + +%% Start processing channels ============================================== +for channelNr = 1:settings.numberOfChannels + + % Only process if the channel's acquisition was successful + if (channel(channelNr).status ~= '-') + % Save additional information - each channel's tracked frequency + % channel + trackResults(channelNr).PRN = channel(channelNr).K; + + % 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, ... + dataAdaptCoeff*(settings.skipNumberOfSamples + channel(channelNr).codePhase-1), ... + 'bof'); + + + %--- 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 + d2CarrError = 0.0; + dCarrError = 0.0; + + %C/No computation + vsmCnt = 0;CNo = 0; + + %=== 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) + + Ln=newline; + trackingStatus=['Tracking: Ch ', int2str(channelNr), ... + ' of ', int2str(settings.numberOfChannels),Ln ... + 'K: ', int2str(channel(channelNr).K),Ln ... + 'Completed ',int2str(loopCnt), ... + ' of ', int2str(codePeriods), ' msec',Ln... + 'C/No: ',CNo,' (dB-Hz)']; + + try + waitbar(loopCnt/codePeriods, ... + hwb, ... + trackingStatus); + 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 ------------------------------------------------ + % Record sample number (based on 8bit samples) + trackResults(channelNr).absoluteSample(loopCnt) =(ftell(fid))/dataAdaptCoeff; + + % 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; + + % Find the size of a "block" or code period in whole samples + blksize = ceil((settings.codeLength-remCodePhase) / codePhaseStep); + + % Read in the appropriate number of samples to process this + % interation + [rawSignal, samplesRead] = fread(fid, ... + dataAdaptCoeff*blksize, settings.dataType); + + rawSignal = rawSignal'; + + if (dataAdaptCoeff==2) + rawSignal1=rawSignal(1:2:end); + rawSignal2=rawSignal(2:2:end); + rawSignal = rawSignal2 + 1i .* rawSignal1; %transpose vector + end + + % If did not read in enough samples, then could be out of + % data - better exit + if (samplesRead ~= dataAdaptCoeff*blksize) + disp('Not able to read the specified number of samples for tracking, exiting!') + delete(hwb); + return + end + + %% Set up all the code phase tracking information ------------------------- + % Save remCodePhase for current correlation + trackResults(channelNr).remCodePhase(loopCnt) = remCodePhase; + + % Define index into early code vector + tcode = (remCodePhase-earlyLateSpc) : ... + codePhaseStep : ... + ((blksize-1)*codePhaseStep+remCodePhase-earlyLateSpc); + tcode2 = ceil(tcode) + 1; + earlyCode = caCode(tcode2); + + % 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); + + remCodePhase = (tcode(blksize) + codePhaseStep) - settings.codeLength; + + %% Generate the carrier frequency to mix the signal to baseband ----------- + % Save remCarrPhase for current correlation + trackResults(channelNr).remCarrPhase(loopCnt) = remCarrPhase; + + % Get the argument to sin/cos functions + time = (0:blksize) ./ settings.samplingFreq; + 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 + carrsig = exp(-1i .* trigarg(1:blksize)); + + %% Generate the six standard accumulated values --------------------------- + % First mix to baseband + iBasebandSignal = real(carrsig .* rawSignal); + qBasebandSignal = imag(carrsig .* 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 + d2CarrError = d2CarrError + carrError * pf3; + dCarrError = d2CarrError + carrError * pf2 + dCarrError; + carrNco = dCarrError + carrError * pf1; + + % Save carrier frequency for current correlation + trackResults(channelNr).carrFreq(loopCnt) = carrFreq; + + % Modify carrier freq based on NCO command + carrFreq = carrFreqBasis + carrNco; + + + + %% 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; + % Save code frequency for current correlation + trackResults(channelNr).codeFreq(loopCnt) = codeFreq; + + % Modify code freq based on NCO command + codeFreq = settings.codeFreqBasis - codeNco; + + %% Record various measures to show in postprocessing ---------------------- + 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; + + if (rem(loopCnt,settings.CNo.VSMinterval)==0) + vsmCnt=vsmCnt+1; + CNoValue=CNoVSM(trackResults(channelNr).I_P(loopCnt-settings.CNo.VSMinterval+1:loopCnt),... + trackResults(channelNr).Q_P(loopCnt-settings.CNo.VSMinterval+1:loopCnt),settings.CNo.accTime); + trackResults(channelNr).CNo.VSMValue(vsmCnt)=CNoValue; + trackResults(channelNr).CNo.VSMIndex(vsmCnt)=loopCnt; + CNo=int2str(CNoValue); + end + + end % for loopCnt + + % 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 K is assigned +end % for channelNr + +% Close the waitbar +close(hwb) + diff --git a/GLO/GLO_GL2/init.m b/GLO/GLO_GL2/init.m new file mode 100644 index 0000000..2636077 --- /dev/null +++ b/GLO/GLO_GL2/init.m @@ -0,0 +1,79 @@ + +%-------------------------------------------------------------------------- +% CU Multi-GNSS SDR +% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% Based on the original work by Darius Plausinaitis,Peter Rinder, +% Nicolaj Bertelsen and Dennis M. Akos +%-------------------------------------------------------------------------- + +%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. +%-------------------------------------------------------------------------- +% +%Script initializes settings and environment of the software receiver. +%Then the processing is started. + +%-------------------------------------------------------------------------- +% CVS record: +% $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $ + +%% Clean up the environment first ========================================= +% clear all; close all; clc; + +format ('compact'); +format ('long', 'g'); + +%--- Include folders with functions --------------------------------------- +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers + +%% Print startup ========================================================== +% fprintf(['\n',... +% 'Welcome to: softGNSS\n\n', ... +% 'An open source GNSS SDR software project initiated by:\n\n', ... +% ' Danish GPS Center/Aalborg University\n\n', ... +% 'The code was improved by GNSS Laboratory/University of Colorado.\n\n',... +% 'The software receiver softGNSS comes with ABSOLUTELY NO WARRANTY;\n',... +% 'for details please read license details in the file license.txt. This\n',... +% 'is free software, and you are welcome to redistribute it under\n',... +% 'the terms described in the license.\n\n']); +% fprintf(' -------------------------------\n\n'); + +%% Initialize constants, settings ========================================= +settings = initSettings(); + +%% Generate plot of raw data and ask if ready to start processing ========= +try + fprintf('Probing data (%s)...\n', settings.fileName) + probeData(settings); +catch + % There was an error, print it and exit + errStruct = lasterror; + disp(errStruct.message); + disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') + return; +end + +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... + postProcessing +end diff --git a/GLO/GLO_GL2/initSettings.m b/GLO/GLO_GL2/initSettings.m index 22a1b3f..e59ace3 100644 --- a/GLO/GLO_GL2/initSettings.m +++ b/GLO/GLO_GL2/initSettings.m @@ -58,7 +58,7 @@ % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = 'data.bin'; +settings.fileName = '/home/gnss/Downloads/GL2_IF20KHz_FS12MHz/B200Ex_1.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -75,7 +75,7 @@ % Intermediate, sampling and code frequencies settings.IF = 0; %[Hz]1602e6-1590e6 -settings.samplingFreq = 10e6; %[Hz] +settings.samplingFreq = 12e6; %[Hz] settings.codeFreqBasis = 0.511e6; %[Hz] % Define number of chips in a code period @@ -131,8 +131,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 1; % 0 - Off - % 1 - On +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] settings.startOffset = 65; %[ms] Initial sign. travel time diff --git a/GPS/GPS_L1CA/include/postProcessing.m b/GPS/GPS_L1CA/include/postProcessing.m index 08ed3e9..01779f7 100644 --- a/GPS/GPS_L1CA/include/postProcessing.m +++ b/GPS/GPS_L1CA/include/postProcessing.m @@ -98,6 +98,7 @@ %--- Do the acquisition ------------------------------------------- disp (' Acquiring satellites...'); acqResults = acquisition(data, settings); + save("acqResults") end %% Initialize channels and prepare for the run ============================ @@ -121,6 +122,7 @@ % Process all channels for given data block [trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") % Close the data file fclose(fid); disp([' Tracking is over (elapsed time ', ... @@ -130,10 +132,17 @@ disp(' Calculating navigation solutions...'); [navResults, ~] = postNavigation(trkResults, settings); +save("navResults") + disp(' Processing is complete for this data block'); disp('Post processing of the signal is over.'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking plotTracking(1:settings.numberOfChannels, trkResults, settings); end diff --git a/GPS/GPS_L1CA/init.m b/GPS/GPS_L1CA/init.m index 3acacff..2636077 100644 --- a/GPS/GPS_L1CA/init.m +++ b/GPS/GPS_L1CA/init.m @@ -1,4 +1,4 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos @@ -54,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -68,13 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); - %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end diff --git a/GPS/GPS_L1CA/initSettings.m b/GPS/GPS_L1CA/initSettings.m index 489fa4d..3a5d317 100644 --- a/GPS/GPS_L1CA/initSettings.m +++ b/GPS/GPS_L1CA/initSettings.m @@ -55,7 +55,7 @@ %% 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 -settings.fileName = '/media/gnss/Ext2TB/JMBF/L1/Lab16Nov_L1schar_Fs5e6_IF20e3.bin'; +settings.fileName = '/home/jmbf/Documents/datasets/L1_IF20Khz_FS18MHz.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -66,7 +66,7 @@ % Intermediate, sampling and code frequencies settings.IF = 20e3; % [Hz] -settings.samplingFreq = 5e6; % [Hz] +settings.samplingFreq = 18e6; % [Hz] settings.codeFreqBasis = 1.023e6; % [Hz] % Define number of chips in a code period @@ -122,9 +122,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 0; % 0 - Off; 1 - On -settings.plotAcquisition = 0; -settings.plotNavigation = 0; +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] settings.startOffset = 68.802; %[ms] Initial sign. travel time diff --git a/GPS/GPS_L2C/include/postProcessing.m b/GPS/GPS_L2C/include/postProcessing.m index c6b94e0..01779f7 100644 --- a/GPS/GPS_L2C/include/postProcessing.m +++ b/GPS/GPS_L2C/include/postProcessing.m @@ -1,4 +1,3 @@ -function [acqResults, trkResults,navResults] = postProcessing(settings, acqPath, trkPath, navPath) % Script postProcessing.m processes the raw signal from the specified data % file (in settings) operating on blocks of 37 seconds of data. % @@ -11,11 +10,10 @@ %-------------------------------------------------------------------------- % CU Multi-GNSS SDR -% (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos +% (C) Updated by Jakob Almqvist, Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- - %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 @@ -44,9 +42,9 @@ % 3.1) Initialize channels (preRun.m). % 3.2) Pass the channel structure and the file identifier to the tracking % function. It will read and process the data. The tracking results are -% stored in the trackResults structure. The results can be accessed this +% stored in the trkResults structure. The results can be accessed this % way (the results are stored each millisecond): -% trackResults(channelNumber).XXX(fromMillisecond : toMillisecond), where +% trkResults(channelNumber).XXX(fromMillisecond : toMillisecond), where % XXX is a field name of the result (e.g. I_P, codePhase etc.) % % 4) Pass tracking results to the navigation solution function. It will @@ -55,86 +53,107 @@ % % 5) Plot the results. - %% Initialization ========================================================= - disp ('Starting processing...'); - - [fid, message] = fopen(settings.fileName, 'rb'); - - %Initialize the multiplier to adjust for the data type - if (settings.fileType==1) - dataAdaptCoeff=1; - else - dataAdaptCoeff=2; - end - - %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, dataAdaptCoeff*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) || ~exist('acqResults', 'var')) - - % Find number of samples per L2 CM spreading code(20 ms) - samplesPerCode = round(settings.samplingFreq / ... - (settings.codeFreqBasis / settings.codeLength)); - - - % Read data for acquisition (3 CM codes should be enough) - data = fread(fid, dataAdaptCoeff*60*samplesPerCode/20, settings.dataType)'; - - if (dataAdaptCoeff == 2) - data1 = data(1:2:end); - data2 = data(2:2:end); - data = data1 + 1i .* data2; - end - - %--- Do the acquisition ------------------------------------------- - disp (' Acquiring satellites...'); - acqResults = acquisition(data, settings); - save(acqPath,"acqResults") - end - - %% Initialize channels and prepare for the run ============================ +%% Initialization ========================================================= +disp ('Starting processing...'); + +[fid, message] = fopen(settings.fileName, 'rb'); + +%Initialize the multiplier to adjust for the data type +if (settings.fileType==1) +dataAdaptCoeff=1; +else +dataAdaptCoeff=2; +end + +%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, dataAdaptCoeff*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) || ~exist('acqResults', 'var')) - % 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.'); - trkResults = []; - navResults = []; - return; - end + % Find number of samples per spreading code + samplesPerCode = round(settings.samplingFreq / ... + (settings.codeFreqBasis / settings.codeLength)); - %% Track the signal ======================================================= - startTime = now; - disp ([' Tracking started at ', datestr(startTime)]); - % Process all channels for given data block - [trkResults, ~] = tracking(fid, channel, settings); - save(trkPath,"trkResults") - % Close the data file - fclose(fid); - - %% Calculate navigation solutions ========================================= - disp(' Calculating navigation solutions...'); - [navResults,~] = postNavigation(trkResults, settings); - save(navPath,"navResults") - disp(' Processing is complete for this data block'); + % At least 42ms of signal are needed for fine frequency estimation + codeLen = max(42,settings.acqNonCohTime+2); + % Read data for acquisition. + data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; + + if (dataAdaptCoeff==2) + data1=data(1:2:end); + data2=data(2:2:end); + data=data1 + 1i .* data2; + end + + %--- Do the acquisition ------------------------------------------- + disp (' Acquiring satellites...'); + acqResults = acquisition(data, settings); + save("acqResults") +end + +%% 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.'); + trkResults = []; + navResults = []; + return; +end + +%% Track the signal ======================================================= +startTime = now; +disp ([' Tracking started at ', datestr(startTime)]); + +% Process all channels for given data block +[trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") +% Close the data file +fclose(fid); +disp([' Tracking is over (elapsed time ', ... + datestr(now - startTime, 13), ')']) + +%% Calculate navigation solutions ========================================= +disp(' Calculating navigation solutions...'); + +[navResults, ~] = postNavigation(trkResults, settings); +save("navResults") + +disp(' Processing is complete for this data block'); +disp('Post processing of the signal is over.'); +%% Plot all results =================================================== +disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + +if settings.plotTracking + plotTracking(1:settings.numberOfChannels, trkResults, settings); +end + +if settings.plotNavigation + plotNavigation(navResults, settings); +end +disp('Post processing of the signal is over.'); + +else +% Error while opening the data file. +error('Unable to read file %s: %s.', settings.fileName, message); +end % if (fid > 0) - - else - % Error while opening the data file. - error('Unable to read file %s: %s.', settings.fileName, message); - end % if (fid > 0) -end \ No newline at end of file diff --git a/GPS/GPS_L2C/init.m b/GPS/GPS_L2C/init.m index 3acacff..453eb6f 100644 --- a/GPS/GPS_L2C/init.m +++ b/GPS/GPS_L2C/init.m @@ -1,4 +1,4 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos @@ -54,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -68,13 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); - %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end \ No newline at end of file diff --git a/GPS/GPS_L2C/initSettings.m b/GPS/GPS_L2C/initSettings.m index b225d57..fb6b200 100644 --- a/GPS/GPS_L2C/initSettings.m +++ b/GPS/GPS_L2C/initSettings.m @@ -1,4 +1,4 @@ - function settings = initSettings() +function settings = initSettings() %Functions initializes and saves settings. Settings can be edited inside of %the function, updated from the command line or updated using a dedicated %GUI - "setSettings". @@ -56,7 +56,7 @@ % This is a "default" name of the data file (signal record) to be used in % the post-processing mode -settings.fileName = 'B200_GPSL2C.bin'; +settings.fileName = '/home/jmbf/Documents/datasets/L2_20KHz_8MHz.bin'; % Data type used to store one sample settings.dataType = 'schar'; @@ -67,7 +67,7 @@ % Intermediate, sampling settings.IF = -20e3; % [Hz] -settings.samplingFreq = 5e6; % [Hz] +settings.samplingFreq = 8e6; % [Hz] %% Code parameter setting % Define number of chips in a code period and code frequencies @@ -129,10 +129,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 0; % 0 - Off - % 1 - On -settings.plotAcquisition = 0; -settings.plotNavigation = 0; +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] settings.startOffset = 68.802; %[ms] Initial sign. travel time diff --git a/GPS/GPS_L5C/include/postProcessing.m b/GPS/GPS_L5C/include/postProcessing.m index 2551e26..01779f7 100644 --- a/GPS/GPS_L5C/include/postProcessing.m +++ b/GPS/GPS_L5C/include/postProcessing.m @@ -60,88 +60,100 @@ %Initialize the multiplier to adjust for the data type if (settings.fileType==1) - dataAdaptCoeff=1; +dataAdaptCoeff=1; else - dataAdaptCoeff=2; +dataAdaptCoeff=2; end %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, dataAdaptCoeff*settings.skipNumberOfBytes, 'bof'); + +% 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, dataAdaptCoeff*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) || ~exist('acqResults', 'var')) - - % Find number of samples per L5 spreading code(1 ms) - samplesPerCode = round(settings.samplingFreq / ... - (settings.codeFreqBasis / settings.codeLength)); - - % Read data for acquisition. At least 22ms of signal are needed - % for the fine frequency estimation - codeLen = max(22,settings.acqNonCohTime+2); - data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; +% Do acquisition if it is not disabled in settings or if the variable +% acqResults does not exist. +if ((settings.skipAcquisition == 0) || ~exist('acqResults', 'var')) + + % Find number of samples per spreading code + samplesPerCode = round(settings.samplingFreq / ... + (settings.codeFreqBasis / settings.codeLength)); - if (dataAdaptCoeff == 2) - data1 = data(1:2:end); - data2 = data(2:2:end); - data = data1 + 1i .* data2; - end - - %--- Do the acquisition ------------------------------------------- - disp (' Acquiring satellites...'); - acqResults = acquisition(data, settings); + + % At least 42ms of signal are needed for fine frequency estimation + codeLen = max(42,settings.acqNonCohTime+2); + % Read data for acquisition. + data = fread(fid, dataAdaptCoeff*codeLen*samplesPerCode, settings.dataType)'; + + if (dataAdaptCoeff==2) + data1=data(1:2:end); + data2=data(2:2:end); + data=data1 + 1i .* data2; end + %--- Do the acquisition ------------------------------------------- + disp (' Acquiring satellites...'); + acqResults = acquisition(data, settings); + save("acqResults") +end + %% 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.'); - trkResults = []; - navResults = []; - return; - end +% 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.'); + trkResults = []; + navResults = []; + return; +end %% Track the signal ======================================================= - startTime = now; - disp ([' Tracking started at ', datestr(startTime)]); - +startTime = now; +disp ([' Tracking started at ', datestr(startTime)]); + +% Process all channels for given data block +[trkResults, ~] = tracking(fid, channel, settings); +save("trkResults") +% Close the data file +fclose(fid); +disp([' Tracking is over (elapsed time ', ... + datestr(now - startTime, 13), ')']) - % Process all channels for given data block - [trkResults, ~] = tracking(fid, channel, settings); - % Close the data file - fclose(fid); - - %% Calculate navigation solutions ========================================= - disp(' Calculating navigation solutions...'); - [navResults,~] = postNavigation(trkResults, settings); - disp(' Processing is complete for this data block'); +disp(' Calculating navigation solutions...'); +[navResults, ~] = postNavigation(trkResults, settings); +save("navResults") + +disp(' Processing is complete for this data block'); +disp('Post processing of the signal is over.'); %% Plot all results =================================================== disp (' Ploting results...'); + +if settings.plotAcquisition + plotAcquisition(acqResults); +end + if settings.plotTracking -plotTracking(1:settings.numberOfChannels, trkResults, settings); + plotTracking(1:settings.numberOfChannels, trkResults, settings); end if settings.plotNavigation -plotNavigation(navResults, settings); + plotNavigation(navResults, settings); end disp('Post processing of the signal is over.'); else - % Error while opening the data file. - error('Unable to read file %s: %s.', settings.fileName, message); +% Error while opening the data file. +error('Unable to read file %s: %s.', settings.fileName, message); end % if (fid > 0) + diff --git a/GPS/GPS_L5C/init.m b/GPS/GPS_L5C/init.m index 88fb344..453eb6f 100644 --- a/GPS/GPS_L5C/init.m +++ b/GPS/GPS_L5C/init.m @@ -1,10 +1,11 @@ -function [navSolutions] = init(filename) + %-------------------------------------------------------------------------- % CU Multi-GNSS SDR % (C) Updated by Yafeng Li, Nagaraj C. Shivaramaiah and Dennis M. Akos % Based on the original work by Darius Plausinaitis,Peter Rinder, % Nicolaj Bertelsen and Dennis M. Akos %-------------------------------------------------------------------------- + %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 @@ -35,9 +36,9 @@ format ('long', 'g'); %--- Include folders with functions --------------------------------------- -addpath include % The software receiver functions -addpath 'Common' % Common functions between differnt SDR receivers -% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers +addpath include % The software receiver functions +addpath Common % Common functions between differnt SDR receivers +% addpath ('../IF_Data_Set') % IF data sets for each SDR receivers %% Print startup ========================================================== % fprintf(['\n',... @@ -53,7 +54,6 @@ %% Initialize constants, settings ========================================= settings = initSettings(); -settings.fileName = filename; %% Generate plot of raw data and ask if ready to start processing ========= try @@ -67,13 +67,13 @@ return; end -% disp(' Raw IF data plotted ') -% disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') -% disp(' '); -% gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); -% -% if (gnssStart == 1) -% disp(' '); -% %start things rolling... +disp(' Raw IF data plotted ') +disp(' (run setSettings or change settings in "initSettings.m" to reconfigure)') +disp(' '); +gnssStart = input('Enter "1" to initiate GNSS processing or "0" to exit : '); + +if (gnssStart == 1) + disp(' '); + % start things rolling... postProcessing -% end +end \ No newline at end of file diff --git a/GPS/GPS_L5C/initSettings.m b/GPS/GPS_L5C/initSettings.m index 179283e..240848e 100644 --- a/GPS/GPS_L5C/initSettings.m +++ b/GPS/GPS_L5C/initSettings.m @@ -52,7 +52,7 @@ %% 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 -settings.fileName = '/media/gnss/Ext2TB/JMBF/L5_Band/B200Fs18MHzIF20KHzExL5_1.bin'; +settings.fileName = '/home/jmbf/Documents/datasets/L5_IF20KHz_FS18MHz/B200Ex_1.bin'; % Data type used to store one sample settings.dataType = 'schar'; % File Types @@ -60,8 +60,8 @@ %2 - 8 bit I/Q samples I0,Q0,I1,Q1,I2,Q2,... settings.fileType = 2; % Intermediate, sampling -settings.IF = -0.05e6; % [Hz] -settings.samplingFreq = 20e6; % [Hz] +settings.IF = -20e3; % [Hz] +settings.samplingFreq = 18e6; % [Hz] %% Code parameter setting % Define number of chips in a code period and code frequencies settings.codeLength = 10230; % L5 code @@ -117,9 +117,9 @@ %% Plot settings ========================================================== % Enable/disable plotting of the tracking results for each channel -settings.plotTracking = 0; -settings.plotNavigation = 0; -settings.plotAcquisition = 0; +settings.plotTracking = 1; % 0 - Off; 1 - On +settings.plotAcquisition = 1; +settings.plotNavigation = 1; %% Constants ============================================================== settings.c = 299792458; % The speed of light, [m/s] settings.startOffset = 68.802; %[ms] Initial sign. travel time