diff --git a/2to3 b/2to3 new file mode 100755 index 0000000..5244470 --- /dev/null +++ b/2to3 @@ -0,0 +1,8 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- +import re +import sys +from cmd_2to3.__main__ import main +if __name__ == '__main__': + sys.argv[0] = re.sub(r'(-script\.pyw|\.exe)?$', '', sys.argv[0]) + sys.exit(main()) diff --git a/README.md b/README.md index 076a073..2e323ea 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ This is a GNSS Software Defined Radio (SDR) implemented in python, based on Soft # System requirements -* python 2.7 +* python 3 * matplotlib * scipy * numpy @@ -15,12 +15,12 @@ Coming soon! # Running the GNSS SDR -1. Examine "main.py" -2. Tweak parameters for the "settings" class if necessary -3. Specify the binary file to be processed and run "main.py" -4. Wait until it is finished +1. run ```python main.py -f=your\file\path -s=(your sampling frequency in Hz) -d=(numpy dtypes, float32 .etc.)``` +2. if you want to modify other settings like skipNumofSamples or set file source to be real signal not complex signal (default complex) check "initialize.py". -# Resources -* The official homepage of the textbook +# Updates +* I/Q file support added +* fileType added to settings, **all datetypes** in **numpy.dtype** are supported +* little bug solved diff --git a/__pycache__/acquisition.cpython-38.pyc b/__pycache__/acquisition.cpython-38.pyc new file mode 100644 index 0000000..1f35440 Binary files /dev/null and b/__pycache__/acquisition.cpython-38.pyc differ diff --git a/__pycache__/acquisition.cpython-39.pyc b/__pycache__/acquisition.cpython-39.pyc new file mode 100644 index 0000000..f4d4695 Binary files /dev/null and b/__pycache__/acquisition.cpython-39.pyc differ diff --git a/__pycache__/ephemeris.cpython-38.pyc b/__pycache__/ephemeris.cpython-38.pyc new file mode 100644 index 0000000..2b250fb Binary files /dev/null and b/__pycache__/ephemeris.cpython-38.pyc differ diff --git a/__pycache__/ephemeris.cpython-39.pyc b/__pycache__/ephemeris.cpython-39.pyc new file mode 100644 index 0000000..5eb427f Binary files /dev/null and b/__pycache__/ephemeris.cpython-39.pyc differ diff --git a/__pycache__/initialize.cpython-38.pyc b/__pycache__/initialize.cpython-38.pyc new file mode 100644 index 0000000..87f12cd Binary files /dev/null and b/__pycache__/initialize.cpython-38.pyc differ diff --git a/__pycache__/initialize.cpython-39.pyc b/__pycache__/initialize.cpython-39.pyc new file mode 100644 index 0000000..a9b614f Binary files /dev/null and b/__pycache__/initialize.cpython-39.pyc differ diff --git a/__pycache__/postNavigation.cpython-38.pyc b/__pycache__/postNavigation.cpython-38.pyc new file mode 100644 index 0000000..003d774 Binary files /dev/null and b/__pycache__/postNavigation.cpython-38.pyc differ diff --git a/__pycache__/postNavigation.cpython-39.pyc b/__pycache__/postNavigation.cpython-39.pyc new file mode 100644 index 0000000..1be1685 Binary files /dev/null and b/__pycache__/postNavigation.cpython-39.pyc differ diff --git a/__pycache__/tracking.cpython-38.pyc b/__pycache__/tracking.cpython-38.pyc new file mode 100644 index 0000000..04472f7 Binary files /dev/null and b/__pycache__/tracking.cpython-38.pyc differ diff --git a/__pycache__/tracking.cpython-39.pyc b/__pycache__/tracking.cpython-39.pyc new file mode 100644 index 0000000..be2792a Binary files /dev/null and b/__pycache__/tracking.cpython-39.pyc differ diff --git a/acquisition.py b/acquisition.py index bc201a5..209a385 100644 --- a/acquisition.py +++ b/acquisition.py @@ -87,13 +87,14 @@ def acquire(self, longSignal): # Correlation peak ratios of the detected signals peakMetric = np.zeros(32) - print '(' + print('(') # Perform search for all listed PRN numbers ... for PRN in range(len(settings.acqSatelliteList)): # Correlate signals ====================================================== # --- Perform DFT of C/A code ------------------------------------------ caCodeFreqDom = np.fft.fft(caCodesTable[PRN, :]).conj() + #Make the correlation for whole frequency band (for all freq. bins) for frqBinIndex in range(numberOfFrqBins): # --- Generate carrier wave frequency grid (0.5kHz step) ----------- frqBins[frqBinIndex] = settings.IF - \ @@ -104,13 +105,15 @@ def acquire(self, longSignal): cosCarr = np.cos(frqBins[frqBinIndex] * phasePoints) - I1 = sinCarr * signal1 + sigCarr=np.exp(1j*frqBins[frqBinIndex] * phasePoints) - Q1 = cosCarr * signal1 + I1 = np.real(sigCarr * signal1) - I2 = sinCarr * signal2 + Q1 = np.imag(sigCarr * signal1) - Q2 = cosCarr * signal2 + I2 = np.real(sigCarr * signal2) + + Q2 = np.imag(sigCarr * signal2) IQfreqDom1 = np.fft.fft(I1 + 1j * Q1) @@ -121,9 +124,10 @@ def acquire(self, longSignal): convCodeIQ2 = IQfreqDom2 * caCodeFreqDom - acqRes1 = abs(np.fft.ifft(convCodeIQ1)) ** 2 + #Perform inverse DFT and store correlation results + acqRes1 = abs(np.fft.ifft(convCodeIQ1)) - acqRes2 = abs(np.fft.ifft(convCodeIQ2)) ** 2 + acqRes2 = abs(np.fft.ifft(convCodeIQ2)) # "blend" 1st and 2nd msec but will correct data bit issues if acqRes1.max() > acqRes2.max(): @@ -142,14 +146,14 @@ def acquire(self, longSignal): peakSize = results.max(0).max() codePhase = results.max(0).argmax() - samplesPerCodeChip = long(round(settings.samplingFreq / settings.codeFreqBasis)) + samplesPerCodeChip = int(round(settings.samplingFreq / settings.codeFreqBasis)) excludeRangeIndex1 = codePhase - samplesPerCodeChip excludeRangeIndex2 = codePhase + samplesPerCodeChip # boundaries - if excludeRangeIndex1 <= 0: + if excludeRangeIndex1 <= 1: codePhaseRange = np.r_[excludeRangeIndex2:samplesPerCode + excludeRangeIndex1 + 1] elif excludeRangeIndex2 >= samplesPerCode - 1: @@ -166,7 +170,7 @@ def acquire(self, longSignal): if (peakSize / secondPeakSize) > settings.acqThreshold: # Fine resolution frequency search ======================================= # --- Indicate PRN number of the detected signal ------------------- - print '%02d ' % (PRN + 1) + print('%02d ' % (PRN + 1)) caCode = settings.generateCAcode(PRN) codeValueIndex = np.floor(ts * np.arange(1, 10 * samplesPerCode + 1) / (1.0 / settings.codeFreqBasis)) @@ -178,26 +182,44 @@ def acquire(self, longSignal): fftNumPts = 8 * 2 ** (np.ceil(np.log2(len(xCarrier)))) + #Compute the magnitude of the FFT, find maximum and the # associated carrier frequency fftxc = np.abs(np.fft.fft(xCarrier, np.long(fftNumPts))) uniqFftPts = np.long(np.ceil((fftNumPts + 1) / 2.0)) - fftMax = fftxc[4:uniqFftPts - 5].max() - fftMaxIndex = fftxc[4:uniqFftPts - 5].argmax() + #fftMax = fftxc[4:uniqFftPts - 5].max() + #fftMaxIndex = fftxc[4:uniqFftPts - 5].argmax() + fftMax = fftxc.max() + fftMaxIndex = fftxc.argmax() fftFreqBins = np.arange(uniqFftPts) * settings.samplingFreq / fftNumPts - carrFreq[PRN] = fftFreqBins[fftMaxIndex] + print('fftxcsize: %4d' %np.size(fftxc)) + if fftMaxIndex>uniqFftPts: + if fftNumPts%2==0: + fftFreqBinsRev=-fftFreqBins[(uniqFftPts-2):0:-1] + fftMax=fftxc[(uniqFftPts):np.size(fftxc)].max() + fftMaxIndex = fftxc[(uniqFftPts):np.size(fftxc)].argmax() + carrFreq[PRN]=-fftFreqBinsRev[fftMaxIndex-1] + else: + fftFreqBinsRev=-fftFreqBins[(uniqFftPts-1):0:-1] + fftMax=fftxc[(uniqFftPts):np.size(fftxc)].max() + fftMaxIndex = fftxc[(uniqFftPts):np.size(fftxc)].argmax() + carrFreq[PRN]=fftFreqBinsRev[fftMaxIndex-1] + else: + carrFreq[PRN]=(-1)**(settings.fileType-1)*fftFreqBins[fftMaxIndex-1] + + #carrFreq[PRN] = fftFreqBins[fftMaxIndex] codePhase_[PRN] = codePhase else: # --- No signal with this PRN -------------------------------------- - print '. ' + print('. ') # === Acquisition is over ================================================== - print ')\n' + print(')\n') acqResults = np.core.records.fromarrays([carrFreq, codePhase_, peakMetric], names='carrFreq,codePhase,peakMetric') self._results = acqResults @@ -209,17 +231,18 @@ def plot(self): import matplotlib.pyplot as plt # from scipy.io.matlab import loadmat + print("start to plot acqusition") # %% configure matplotlib - mpl.rcdefaults() - # mpl.rcParams['font.sans-serif'] - # mpl.rcParams['font.family'] = 'serif' - mpl.rc('savefig', bbox='tight', transparent=False, format='png') - mpl.rc('axes', grid=True, linewidth=1.5, axisbelow=True) - mpl.rc('lines', linewidth=1.5, solid_joinstyle='bevel') - mpl.rc('figure', figsize=[8, 6], autolayout=False, dpi=120) - mpl.rc('text', usetex=True) - mpl.rc('font', family='serif', serif='Computer Modern Roman', size=16) - mpl.rc('mathtext', fontset='cm') + # mpl.rcdefaults() + # # mpl.rcParams['font.sans-serif'] + # # mpl.rcParams['font.family'] = 'serif' + # mpl.rc('savefig', bbox='tight', transparent=False, format='png') + # mpl.rc('axes', grid=True, linewidth=1.5, axisbelow=True) + # mpl.rc('lines', linewidth=1.5, solid_joinstyle='bevel') + # mpl.rc('figure', figsize=[8, 6], autolayout=False, dpi=120) + # mpl.rc('text', usetex=True) + # mpl.rc('font', family='serif', serif='Computer Modern Roman', size=16) + # mpl.rc('mathtext', fontset='cm') # mpl.rc('font', size=16) # mpl.rc('text.latex', preamble=r'\usepackage{cmbright}') @@ -249,7 +272,7 @@ def plot(self): hAxes.xaxis.grid() # Mark acquired signals ================================================== - acquiredSignals = self.peakMetric * (self.carrFreq > 0) + acquiredSignals = self.peakMetric * (self.peakMetric > self._settings.acqThreshold) plt.bar(range(1, 33), acquiredSignals, FaceColor=(0, 0.8, 0)) plt.legend(['Not acquired signals', 'Acquired signals']) @@ -291,7 +314,7 @@ def preRun(self): # --- Load information about each satellite -------------------------------- # Maximum number of initialized channels is number of detected signals, but # not more as the number of channels specified in the settings. - for ii in range(min(settings.numberOfChannels, sum(self.carrFreq > 0))): + for ii in range(min(settings.numberOfChannels, sum(self.peakMetric > self._settings.acqThreshold))): PRN[ii] = PRNindexes[ii][0] + 1 acquiredFreq[ii] = self.carrFreq[PRNindexes[ii][0]] @@ -318,22 +341,22 @@ def showChannelStatus(self): channel = self._channels settings = self._settings assert isinstance(channel, np.recarray) - print ('\n*=========*=====*===============*===========*=============*========*') - print ('| Channel | PRN | Frequency | Doppler | Code Offset | Status |') - print ('*=========*=====*===============*===========*=============*========*') + print('\n*=========*=====*===============*===========*=============*========*') + print('| Channel | PRN | Frequency | Doppler | Code Offset | Status |') + print('*=========*=====*===============*===========*=============*========*') for channelNr in range(settings.numberOfChannels): if channel[channelNr].status != '-': - print '| %2d | %3d | %2.5e | %5.0f | %6d | %1s |' % ( + print('| %2d | %3d | %2.5e | %5.0f | %6d | %1s |' % ( channelNr, channel[channelNr].PRN, channel[channelNr].acquiredFreq, channel[channelNr].acquiredFreq - settings.IF, channel[channelNr].codePhase, - channel[channelNr].status) + channel[channelNr].status)) else: - print '| %2d | --- | ------------ | ----- | ------ | Off |' % channelNr + print('| %2d | --- | ------------ | ----- | ------ | Off |' % channelNr) - print '*=========*=====*===============*===========*=============*========*\n' + print('*=========*=====*===============*===========*=============*========*\n') if __name__ == '__main__': diff --git a/geoFunctions/__init__.py b/geoFunctions/__init__.py index 4b3023d..668cf24 100644 --- a/geoFunctions/__init__.py +++ b/geoFunctions/__init__.py @@ -54,7 +54,7 @@ def cart2geo(X, Y, Z, i, *args, **kwargs): iterations += 1 if iterations > 100: - print 'Failed to approximate h with desired precision. h-oldh: %e.' % (h - oldh) + print ('Failed to approximate h with desired precision. h-oldh: %e.' % (h - oldh)) break phi *= (180 / np.pi) @@ -264,7 +264,7 @@ def cart2utm(X, Y, Z, zone, *args, **kwargs): iterations += 1 if iterations > 100: - print 'Failed to approximate U with desired precision. U-oldU: %e.' % (U - oldU) + print ('Failed to approximate U with desired precision. U-oldU: %e.' % (U - oldU)) break # Normalized meridian quadrant, KW p. 50 (96), p. 19 (38b), p. 5 (21) @@ -991,7 +991,7 @@ def togeod(a, finv, X, Y, Z, *args, **kwargs): break # Not Converged--Warn user if i == maxit - 1: - print ' Problem in TOGEOD, did not converge in %2.0f iterations' % i + print (' Problem in TOGEOD, did not converge in %2.0f iterations' % i) dphi *= rtd return dphi, dlambda, h @@ -1192,4 +1192,4 @@ def tropo(sinel, hsta, p, tkel, hum, hp, htkel, hhum): # print "This program is being run by itself" pass else: - print 'Importing functions from ./geoFunctions/' + print ('Importing functions from ./geoFunctions/') diff --git a/geoFunctions/__pycache__/__init__.cpython-38.pyc b/geoFunctions/__pycache__/__init__.cpython-38.pyc new file mode 100644 index 0000000..deccde2 Binary files /dev/null and b/geoFunctions/__pycache__/__init__.cpython-38.pyc differ diff --git a/geoFunctions/__pycache__/__init__.cpython-39.pyc b/geoFunctions/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000..5135126 Binary files /dev/null and b/geoFunctions/__pycache__/__init__.cpython-39.pyc differ diff --git a/initialize.py b/initialize.py index 027e761..308d788 100644 --- a/initialize.py +++ b/initialize.py @@ -82,31 +82,36 @@ def __init__(self): # Processing settings ==================================================== # Number of milliseconds to be processed used 36000 + any transients (see # below - in Nav parameters) to ensure nav subframes are provided - self.msToProcess = 37000.0 + self.msToProcess = 10000.0 # Number of channels to be used for signal processing - self.numberOfChannels = 8 + self.numberOfChannels = 2 - # Move the starting point of processing. Can be used to start the signal - # processing at any point in the data record (e.g. for long records). fseek - # function is used to move the file read point, therefore advance is byte - # based only. - self.skipNumberOfBytes = 0 + # modifyed by mortarboard in 2022 09 20 + # self.skipNumberOfBytes are calculated according to the number of samples and datatype + self.skipNumberOfSamples=0 + self.skipNumberOfSamples=int(self.skipNumberOfSamples) + + self.fileType=2 + #1 for real signal, 2 for I/Q signal # 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 - self.fileName = '/Users/yangsu/Downloads/GNSS_signal_records/GPSdata-DiscreteComponents-fs38_192-if9_55.bin' + self.fileName = 'test.bin' # Data type used to store one sample - self.dataType = 'int8' + self.dataType = 'float32' + + # calculate skip bytes + self.skipNumberOfBytes = int(self.skipNumberOfSamples*self.fileType*np.dtype(self.dataType).itemsize) # Intermediate, sampling and code frequencies - self.IF = 9548000.0 + self.IF = 0 - self.samplingFreq = 38192000.0 + self.samplingFreq = 2.6e6 - self.codeFreqBasis = 1023000.0 + self.codeFreqBasis = 1.023e6 # Define number of chips in a code period self.codeLength = 1023 @@ -120,10 +125,10 @@ def __init__(self): self.acqSatelliteList = range(1, 33) # Band around IF to search for satellite signal. Depends on max Doppler - self.acqSearchBand = 14.0 + self.acqSearchBand = 20.0 #(KHz) # Threshold for the signal presence decision rule - self.acqThreshold = 2.5 + self.acqThreshold = 1.4 # Tracking loops settings ================================================ # Code tracking loop parameters @@ -131,7 +136,7 @@ def __init__(self): self.dllNoiseBandwidth = 2.0 - self.dllCorrelatorSpacing = 0.5 + self.dllCorrelatorSpacing = 0.5 #(chips) # Carrier tracking loop parameters self.pllDampingRatio = 0.7 @@ -331,7 +336,7 @@ def probeData(self, fileNameStr=None): import matplotlib.pyplot as plt from scipy.signal import welch - from scipy.signal.windows.windows import hamming + from scipy.signal.windows import hamming # Function plots raw data information: time domain plot, a frequency domain # plot and a histogram. @@ -362,59 +367,118 @@ def probeData(self, fileNameStr=None): # 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). - fid.seek(settings.skipNumberOfBytes, 0) + fid.seek(int(settings.skipNumberOfBytes), 0) samplesPerCode = settings.samplesPerCode + if self.fileType==1: + dataAdaptCoeff=1 + else: + dataAdaptCoeff=2 + + try: data = np.fromfile(fid, settings.dataType, - 10 * samplesPerCode) + dataAdaptCoeff*100 * samplesPerCode) + dataLength=np.size(data) + + if self.fileType==2:#if IQ file + data1=data[0:dataLength:2] + data2=data[1:dataLength:2] + data2=data2.astype(complex) + data2.imag=data2.real + data2.real=data1 + data=data2 + except IOError: # The file is too short - print 'Could not read enough data from the data file.' + print('Could not read enough data from the data file.') # --- Initialization --------------------------------------------------- plt.figure(100) plt.clf() timeScale = np.arange(0, 0.005, 1 / settings.samplingFreq) - plt.subplot(2, 2, 1) - plt.plot(1000 * timeScale[1:samplesPerCode / 50], - data[1:samplesPerCode / 50]) - plt.axis('tight') - plt.grid() - plt.title('Time domain plot') - plt.xlabel('Time (ms)') - plt.ylabel('Amplitude') - plt.subplot(2, 2, 2) - f, Pxx = welch(data - np.mean(data), - settings.samplingFreq / 1000000.0, - hamming(16384, False), - 16384, - 1024, - 16384) - plt.semilogy(f, Pxx) - plt.axis('tight') - plt.grid() - plt.title('Frequency domain plot') - plt.xlabel('Frequency (MHz)') - plt.ylabel('Magnitude') - plt.show() - plt.subplot(2, 2, 3.5) - plt.hist(data, np.arange(- 128, 128)) - dmax = np.max(np.abs(data)) + 1 - - plt.axis('tight') - adata = plt.axis() - - plt.axis([-dmax, dmax, adata[2], adata[3]]) - plt.grid('on') - plt.title('Histogram') - plt.xlabel('Bin') - plt.ylabel('Number in bin') + #for real signal + if self.fileType==1: + #Time Domain plot + plt.subplot(2, 2, 1) + plt.plot(1000 * timeScale[1:samplesPerCode], + data[1:samplesPerCode]) + plt.axis('tight') + plt.grid() + plt.title('Time domain plot') + plt.xlabel('Time (ms)') + plt.ylabel('Amplitude') + + #Frequency Domain plot + plt.subplot(2, 2, 2) + f, Pxx = welch(data - np.mean(data), + settings.samplingFreq / 1000000.0, + hamming(16384, False), + 16384, + 1024, + 16384) + plt.semilogy(f, Pxx) + plt.axis('tight') + plt.grid() + plt.title('Frequency domain plot') + plt.xlabel('Frequency (MHz)') + plt.ylabel('Magnitude') + plt.show() + + #Histogram plot + plt.subplot(2, 2, 3.5) + plt.hist(data, np.arange(- 128, 128)) + dmax = np.max(np.abs(data)) + 1 + + plt.axis('tight') + adata = plt.axis() + + plt.axis([-dmax, dmax, adata[2], adata[3]]) + plt.grid('on') + plt.title('Histogram') + plt.xlabel('Bin') + plt.ylabel('Number in bin') + # for IQ signal + else: + #Time Domain plot + plt.subplot(2, 2, 1) + plt.plot(1000 * timeScale[1:samplesPerCode], + data.real[1:samplesPerCode]) + plt.axis('tight') + plt.grid() + plt.title('Time domain plot (Real/I)') + plt.xlabel('Time (ms)') + plt.ylabel('Amplitude') + + plt.subplot(2, 2, 2) + plt.plot(1000 * timeScale[1:samplesPerCode], + data.imag[1:samplesPerCode]) + plt.axis('tight') + plt.grid() + plt.title('Time domain plot (Imag/Q)') + plt.xlabel('Time (ms)') + plt.ylabel('Amplitude') + + #Frequency Domain plot + plt.subplot(2, 2, 3.5) + f, Pxx = welch(data - np.mean(data), + settings.samplingFreq / 1000000.0, + hamming(16384, False), + 16384, + 1024, + 16384) + plt.semilogy(f, Pxx) + plt.axis('tight') + plt.grid() + plt.title('Frequency domain plot') + plt.xlabel('Frequency (MHz)') + plt.ylabel('Magnitude') + plt.show() # === Error while opening the data file ================================ - except IOError as e: - print 'Unable to read file "%s": %s' % (fileNameStr, e) + except: + print("unable to read file") # ./postProcessing.m @@ -456,7 +520,7 @@ def postProcessing(self, fileNameStr=None): import acquisition import postNavigation import tracking - print 'Starting processing...' + print('Starting processing...') settings = self if not fileNameStr: fileNameStr = settings.fileName @@ -469,7 +533,8 @@ def postProcessing(self, fileNameStr=None): # 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). - fid.seek(settings.skipNumberOfBytes, 0) + fid.seek(int(settings.skipNumberOfBytes), 0) + # Acquisition ============================================================ # Do acquisition if it is not disabled in settings or if the variable # acqResults does not exist. @@ -477,13 +542,29 @@ def postProcessing(self, fileNameStr=None): # Find number of samples per spreading code samplesPerCode = settings.samplesPerCode + if self.fileType==1: + dataAdaptCoeff=1 + else: + dataAdaptCoeff=2 + # frequency estimation - data = np.fromfile(fid, settings.dataType, 11 * samplesPerCode) + data = np.fromfile(fid, settings.dataType, dataAdaptCoeff*11 * samplesPerCode) + + dataLength=np.size(data) + + if self.fileType==2:#if IQ file + data1=data[0:dataLength:2] + data2=data[1:dataLength:2] + data2=data2.astype(complex) + data2.imag=data2.real + data2.real=data1 + data=data2 - print ' Acquiring satellites...' + + #print' Acquiring satellites...' acqResults = acquisition.AcquisitionResult(settings) acqResults.acquire(data) - acqResults.plot() + #acqResults.plot() # 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) @@ -492,36 +573,42 @@ def postProcessing(self, fileNameStr=None): acqResults.showChannelStatus() else: # No satellites to track, exit - print 'No GNSS signals detected, signal processing finished.' + #print'No GNSS signals detected, signal processing finished.' trackResults = None + acqResults.plot() + print("Press Enter to continue...") # Track the signal ======================================================= startTime = datetime.datetime.now() - print ' Tracking started at %s' % startTime.strftime('%X') + #print' Tracking started at %s' % startTime.strftime('%X') trackResults = tracking.TrackingResult(acqResults) + trackResults.track(fid) + trackResults.plot() try: trackResults.results = np.load('trackingResults_python.npy') - except IOError: + + except : trackResults.track(fid) np.save('trackingResults_python', trackResults.results) - print ' Tracking is over (elapsed time %s s)' % (datetime.datetime.now() - startTime).total_seconds() + #print' Tracking is over (elapsed time %s s)' % (datetime.datetime.now() - startTime).total_seconds() # Auto save the acquisition & tracking results to save time. - print ' Saving Acquisition & Tracking results to storage' + #print' Saving Acquisition & Tracking results to storage' + input("Press Enter to continue...") # Calculate navigation solutions ========================================= - print ' Calculating navigation solutions...' + #print' Calculating navigation solutions...' navResults = postNavigation.NavigationResult(trackResults) navResults.postNavigate() - print ' Processing is complete for this data block' + #print' Processing is complete for this data block' # Plot all results =================================================== - print ' Plotting results...' + #print' Plotting results...' # TODO turn off tracking plots for now if not settings.plotTracking: trackResults.plot() navResults.plot() - print 'Post processing of the signal is over.' + #print'Post processing of the signal is over.' except IOError as e: # Error while opening the data file. - print 'Unable to read file "%s": %s.' % (settings.fileName, e) + print('Unable to read file "%s": %s.' % (settings.fileName, e)) diff --git a/initialize.pyc b/initialize.pyc new file mode 100644 index 0000000..3bac503 Binary files /dev/null and b/initialize.pyc differ diff --git a/main.py b/main.py index a8f92da..137a902 100644 --- a/main.py +++ b/main.py @@ -1,5 +1,5 @@ import initialize - +import numpy as np # ./init.m # -------------------------------------------------------------------------- @@ -39,34 +39,50 @@ # --- Include folders with functions --------------------------------------- # addpath('include') # addpath('geoFunctions') -# Print startup ========================================================== -print '\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', \ - ' -------------------------------\n\n' +# #printstartup ========================================================== +#print'\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', \ + # ' -------------------------------\n\n' # Initialize settings class========================================= settings = initialize.Settings() +# add arguments parsing to override settings +import argparse +parser=argparse.ArgumentParser() +parser.add_argument("-f","--file",help="the source file") +parser.add_argument("-s","--samp_rate",help="sampling frequency",type=int) +parser.add_argument("-d","--dtype",help="datatype, numpy dtypes") +args=parser.parse_args() + +if(args.file): + settings.fileName=args.file +if(args.samp_rate): + settings.samplingFreq=args.samp_rate +if(args.dtype): + settings.dataType=args.dtype + settings.skipNumberOfBytes=int(settings.skipNumberOfSamples*settings.fileType*np.dtype(settings.dataType).itemsize) + # Generate plot of raw data and ask if ready to start processing ========= try: - print 'Probing data "%s"...' % settings.fileName + print('Probing data "%s"...' % settings.fileName) settings.probeData() - settings.probeData('/Users/yangsu/Downloads/GNSS_signal_records/GPS_and_GIOVE_A-NN-fs16_3676-if4_1304.bin') + #settings.probeData('/Users/yangsu/Downloads/GNSS_signal_records/GPS_and_GIOVE_A-NN-fs16_3676-if4_1304.bin') finally: pass -print ' Raw IF data plotted ' -print ' (run setSettings or change settings in "initialize.py" to reconfigure)' -print ' ' +print(' Raw IF data plotted ') +print(' (run setSettings or change settings in "initialize.py" to reconfigure)') +print(' ') gnssStart = True # gnssStart = int(raw_input('Enter "1" to initiate GNSS processing or "0" to exit : ').strip()) if gnssStart: - print ' ' + #print' ' settings.postProcessing() diff --git a/postNavigation.py b/postNavigation.py index f9a714f..91f16d1 100644 --- a/postNavigation.py +++ b/postNavigation.py @@ -10,7 +10,7 @@ def __init__(self, trackResult): self._results = trackResult.results self._channels = trackResult.channels self._settings = trackResult.settings - self._solutions = None + self._solutions = No self._eph = None @property @@ -103,7 +103,7 @@ def postNavigate(self): if settings.msToProcess < 36000 or sum(trackResults.status != '-') < 4: # Show the error message and exit - print 'Record is to short or too few satellites tracked. Exiting!' + #print'Record is to short or too few satellites tracked. Exiting!' navSolutions = None self._solutions = navSolutions eph = None @@ -148,7 +148,7 @@ def postNavigate(self): # Check if the number of satellites is still above 3 ===================== if activeChnList.size == 0 or activeChnList.size < 4: # Show error message and exit - print 'Too few satellites with ephemeris data for position calculations. Exiting!' + #print'Too few satellites with ephemeris data for position calculations. Exiting!' navSolutions = None self._solutions = navSolutions eph = None @@ -263,7 +263,7 @@ def postNavigate(self): else: # --- There are not enough satellites to find 3D position ---------- - print ' Measurement No. %d' % currMeasNr + ': Not enough information for position solution.' + #print' Measurement No. %d' % currMeasNr + ': Not enough information for position solution.' # excluded automatically in all plots. For DOP it is easier to use # zeros. NaN values might need to be excluded from results in some # of further processing to obtain correct results. @@ -436,7 +436,7 @@ def plot(self): h32.set_title('Sky plot (mean PDOP: %f )' % np.mean(navSolutions[0].DOP[1, :])) f.show() else: - print 'plotNavigation: No navigation data to plot.' + print('plotNavigation: No navigation data to plot.') @staticmethod # navPartyChk.m @@ -627,7 +627,7 @@ def findPreambles(self): # valid preamble and therefore nothing more can be done for it. activeChnList = np.setdiff1d(activeChnList, channelNr) - print 'Could not find valid preambles in channel %2d !' % channelNr + #print'Could not find valid preambles in channel %2d !' % channelNr return firstSubFrame, activeChnList diff --git a/tracking.py b/tracking.py index 649de40..c692b75 100644 --- a/tracking.py +++ b/tracking.py @@ -52,12 +52,20 @@ def track(self, fid): tau1carr, tau2carr = settings.calcLoopCoef(settings.pllNoiseBandwidth, settings.pllDampingRatio, 0.25) # hwb=waitbar(0,'Tracking...') + if settings.fileType==1: + dataAdaptCoeff=1 + else: + dataAdaptCoeff=2 + # Initialize a temporary list of records rec = [] # Start processing channels ============================================== for channelNr in range(settings.numberOfChannels): + + msToProcess = np.long(settings.msToProcess) + # Initialize fields for record(structured) array of tracked results status = '-' @@ -97,24 +105,44 @@ def track(self, fid): # Only process if PRN is non zero (acquisition was successful) if channel[channelNr].PRN != 0: + #print' starting to trackc PRN: %2d' %channel[channelNr].PRN # Save additional information - each channel's tracked PRN PRN = channel[channelNr].PRN + + # 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) - fid.seek(settings.skipNumberOfBytes + channel[channelNr].codePhase, 0) + + sbCoeff=np.dtype(settings.dataType).itemsize + # if settings.dataType=='float32': + # sbCoeff=4 + # if settings.dataType=='single': + # sbCoeff=4 + # if settings.dataType=='byte': + # sbCoeff=1 + #if settings.dataType=='float': + # sbCoeff=4 + #if settings.dataType=='float': + # sbCoeff=4 + #if settings.dataType=='float': + # sbCoeff=4 + + fid.seek(int(np.round(sbCoeff*dataAdaptCoeff*(settings.skipNumberOfSamples + channel[channelNr].codePhase))), 0) # Here PRN is the actual satellite ID instead of the 0-based index - caCode = settings.generateCAcode(channel[channelNr].PRN - 1) + caCode = settings.generateCAcode(channel[channelNr].PRN-1) - caCode = np.r_[caCode[-1], caCode, caCode[0]] + caCode = np.r_[caCode[-1], caCode, caCode[0], caCode[1]] # 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 @@ -136,9 +164,9 @@ def track(self, fid): # all the time with GUI task. if loopCnt % 50 == 0: try: - print 'Tracking: Ch %d' % (channelNr + 1) + ' of %d' % settings.numberOfChannels + \ + print('Tracking: Ch %d' % (channelNr + 1) + ' of %d' % settings.numberOfChannels + \ '; PRN#%02d' % channel[channelNr].PRN + \ - '; Completed %d' % loopCnt + ' of %d' % codePeriods + ' msec' + '; Completed %d' % loopCnt + ' of %d' % codePeriods + ' msec') finally: pass # Read next block of data ------------------------------------------------ @@ -151,37 +179,48 @@ def track(self, fid): blksize = np.long(blksize) # interaction - rawSignal = np.fromfile(fid, settings.dataType, blksize) + rawSignal = np.fromfile(fid, settings.dataType, dataAdaptCoeff*blksize) samplesRead = len(rawSignal) + + if dataAdaptCoeff==2: + rawSignal1=rawSignal[0:samplesRead:2] + rawSignal2=rawSignal[1:samplesRead:2] + rawSignal2=rawSignal2.astype(complex) + rawSignal2.imag=rawSignal2.real + rawSignal2.real=rawSignal1 + rawSignal=rawSignal2 + # If did not read in enough samples, then could be out of # data - better exit - if samplesRead != blksize: - print 'Not able to read the specified number of samples for tracking, exiting!' + if samplesRead != dataAdaptCoeff*blksize: + print('Not able to read the specified number of samples for tracking, exiting!') fid.close() trackResults = None return trackResults # Set up all the code phase tracking information ------------------------- # Define index into early code vector tcode = np.linspace(remCodePhase - earlyLateSpc, - blksize * codePhaseStep + remCodePhase - earlyLateSpc, - blksize, endpoint=False) + (blksize-1) * codePhaseStep + remCodePhase - earlyLateSpc, + blksize, endpoint=True) tcode2 = np.ceil(tcode) earlyCode = caCode[np.longlong(tcode2)] - + + # Define index into late code vector tcode = np.linspace(remCodePhase + earlyLateSpc, - blksize * codePhaseStep + remCodePhase + earlyLateSpc, - blksize, endpoint=False) + (blksize-1) * codePhaseStep + remCodePhase + earlyLateSpc, + blksize, endpoint=True) tcode2 = np.ceil(tcode) lateCode = caCode[np.longlong(tcode2)] + # Define index into prompt code vector tcode = np.linspace(remCodePhase, - blksize * codePhaseStep + remCodePhase, - blksize, endpoint=False) + (blksize-1) * codePhaseStep + remCodePhase, + blksize, endpoint=True) tcode2 = np.ceil(tcode) @@ -189,6 +228,8 @@ def track(self, fid): remCodePhase = tcode[blksize - 1] + codePhaseStep - 1023.0 + + # Generate the carrier frequency to mix the signal to baseband ----------- time = np.arange(0, blksize + 1) / settings.samplingFreq @@ -200,11 +241,12 @@ def track(self, fid): carrSin = np.sin(trigarg[0:blksize]) + carrsig=np.exp(1j*trigarg[0:blksize]) # Generate the six standard accumulated values --------------------------- # First mix to baseband - qBasebandSignal = carrCos * rawSignal + qBasebandSignal = np.real(carrsig * rawSignal) - iBasebandSignal = carrSin * rawSignal + iBasebandSignal = np.imag(carrsig * rawSignal) I_E = (earlyCode * iBasebandSignal).sum() @@ -230,6 +272,7 @@ def track(self, fid): oldCarrError = carrError + #Modify carrier freq based on NCO command carrFreq = carrFreqBasis + carrNco carrFreq_[loopCnt] = carrFreq @@ -238,6 +281,7 @@ def track(self, fid): codeError = (np.sqrt(I_E * I_E + Q_E * Q_E) - np.sqrt(I_L * I_L + Q_L * Q_L)) / ( np.sqrt(I_E * I_E + Q_E * Q_E) + np.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) @@ -246,6 +290,7 @@ def track(self, fid): oldCodeError = codeError + #Modify code freq based on NCO command codeFreq = settings.codeFreqBasis - codeNco codeFreq_[loopCnt] = codeFreq @@ -281,6 +326,7 @@ def track(self, fid): rec.append((status, absoluteSample, codeFreq_, carrFreq_, I_P_, I_E_, I_L_, Q_E_, Q_P_, Q_L_, dllDiscr, dllDiscrFilt, pllDiscr, pllDiscrFilt, PRN)) + ##print'tracked a signal %2d' %np.size(rec) trackResults = np.rec.fromrecords(rec, dtype=[('status', 'S1'), ('absoluteSample', 'object'), ('codeFreq', 'object'), @@ -298,16 +344,16 @@ def plot(self): import matplotlib as mpl # %% configure matplotlib - mpl.rcdefaults() - # mpl.rcParams['font.sans-serif'] - # mpl.rcParams['font.family'] = 'serif' - mpl.rc('savefig', bbox='tight', transparent=False, format='png') - mpl.rc('axes', grid=True, linewidth=1.5, axisbelow=True) - mpl.rc('lines', linewidth=1.5, solid_joinstyle='bevel') - mpl.rc('figure', figsize=[8, 6], dpi=120) - mpl.rc('text', usetex=True) - mpl.rc('font', family='serif', serif='Computer Modern Roman', size=8) - mpl.rc('mathtext', fontset='cm') + # mpl.rcdefaults() + # # mpl.rcParams['font.sans-serif'] + # # mpl.rcParams['font.family'] = 'serif' + # mpl.rc('savefig', bbox='tight', transparent=False, format='png') + # mpl.rc('axes', grid=True, linewidth=1.5, axisbelow=True) + # mpl.rc('lines', linewidth=1.5, solid_joinstyle='bevel') + # mpl.rc('figure', figsize=[8, 6], dpi=120) + # mpl.rc('text', usetex=True) + # mpl.rc('font', family='serif', serif='Computer Modern Roman', size=8) + # mpl.rc('mathtext', fontset='cm') # mpl.rc('font', size=16) # mpl.rc('text.latex', preamble=r'\usepackage{cmbright}') @@ -318,21 +364,22 @@ def plot(self): settings = self._settings channelList = range(settings.numberOfChannels) + import matplotlib as mpl import matplotlib.gridspec as gs import matplotlib.pyplot as plt # %% configure matplotlib - mpl.rcdefaults() - # mpl.rcParams['font.sans-serif'] - # mpl.rcParams['font.family'] = 'serif' - mpl.rc('savefig', bbox='tight', transparent=False, format='png') - mpl.rc('axes', grid=True, linewidth=1.5, axisbelow=True) - mpl.rc('lines', linewidth=1.5, solid_joinstyle='bevel') - mpl.rc('figure', figsize=[8, 6], dpi=120) - mpl.rc('text', usetex=True) - mpl.rc('font', family='serif', serif='Computer Modern Roman', size=8) - mpl.rc('mathtext', fontset='cm') + # mpl.rcdefaults() + # # mpl.rcParams['font.sans-serif'] + # # mpl.rcParams['font.family'] = 'serif' + # mpl.rc('savefig', bbox='tight', transparent=False, format='png') + # mpl.rc('axes', grid=True, linewidth=1.5, axisbelow=True) + # mpl.rc('lines', linewidth=1.5, solid_joinstyle='bevel') + # mpl.rc('figure', figsize=[8, 6], dpi=120) + # mpl.rc('text', usetex=True) + # mpl.rc('font', family='serif', serif='Computer Modern Roman', size=8) + # mpl.rc('mathtext', fontset='cm') # mpl.rc('font', size=16) # mpl.rc('text.latex', preamble=r'\usepackage{cmbright}') @@ -350,14 +397,16 @@ def plot(self): # Protection - if the list contains incorrect channel numbers channelList = np.intersect1d(channelList, range(settings.numberOfChannels)) - + channel = self._channels # === For all listed channels ============================================== - for channelNr in channelList: + for channelNr in range(np.size(trackResults)): # Select (or create) and clear the figure ================================ # The number 200 is added just for more convenient handling of the open # figure windows, when many figures are closed and reopened. # Figures drawn or opened by the user, will not be "overwritten" by # this function. + print(' %2d' %np.size(trackResults)) + f = plt.figure(channelNr + 200) f.set_label('Channel ' + str(channelNr) + ' (PRN ' + str(trackResults[channelNr].PRN) + ') results') @@ -424,3 +473,6 @@ def plot(self): ylabel='Amplitude', title='Filtered DLL discriminator') f.show() + + #something just get past the period + diff --git a/trackingResults_python.npy b/trackingResults_python.npy new file mode 100644 index 0000000..2b88593 Binary files /dev/null and b/trackingResults_python.npy differ