www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/acquisition_1ms.m
function acqResults = acquisition_1ms(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 - 11 ms of raw 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.
%--------------------------------------------------------------------------
% SoftGNSS v3.0
%
% Copyright (C) Darius Plausinaitis and Dennis M. Akos
% Written by Darius Plausinaitis and Dennis M. Akos
% Based on Peter Rinder and Nicolaj Bertelsen
%--------------------------------------------------------------------------
%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 $
%% Initialization ================================================
% Find number of samples per spreading code
samplesPerCode = round(settings.samplingFreq / ...
(settings.codeFreqBasis / settings.codeLength));
% Create two 1msec vectors of data to correlate with and one with zero DC
signal1 = longSignal(1 : samplesPerCode);
% signal2 = longSignal(samplesPerCode+1 : 2*samplesPerCode);
% signal3 = longSignal(2*samplesPerCode+1 : 3*samplesPerCode);
% signal4 = longSignal(3*samplesPerCode+1 : 4*samplesPerCode);
% signal5 = longSignal(4*samplesPerCode+1 : 5*samplesPerCode);
% signal6 = longSignal(5*samplesPerCode+1 : 6*samplesPerCode);
signal0DC = longSignal - mean(longSignal);
%--------------------------------------------------------------------------
% Find sampling period
ts = 1 / settings.samplingFreq;
%--------------------------------------------------------------------------
% Find phase points of the local carrier wave
phasePoints = (0 : (samplesPerCode-1)) * 2 * pi * ts;
%--------------------------------------------------------------------------
% Number of the frequency bins for the given acquisition band (500Hz steps)
numberOfFrqBins = round(settings.acqSearchBand * 2) + 1;
%--------------------------------------------------------------------------
% Generate all C/A codes and sample them according to the sampling freq.
caCodesTable = makeCaTable(settings);
%--------------------------------------------------------------------------
% Initialize arrays to speed up the code
% Search results of all frequency bins and code shifts (for one satellite)
results = zeros(numberOfFrqBins, samplesPerCode);
% Carrier frequencies of the frequency bins
frqBins = zeros(1, numberOfFrqBins);
%--------------------------------------------------------------------------
% Initialize acqResults
% Carrier frequencies of detected signals
acqResults.carrFreq = zeros(1, 32);
%--------------------------------------------------------------------------
% C/A code phases of detected signals
acqResults.codePhase = zeros(1, 32);
%--------------------------------------------------------------------------
% Correlation peak ratios of the detected signals
acqResults.peakMetric = zeros(1, 32);
%--------------------------------------------------------------------------
fprintf('(');
%--------------------------------------------------------------------------
% Perform search for all listed PRN numbers ...
for PRN = settings.acqSatelliteList
%% Correlate signals =============================================
%--- Perform DFT of C/A code ------------------------------------------
caCodeFreqDom = conj(fft(caCodesTable(PRN, :)));
% caCodeFreqDom=repmat(caCodeFreqDom,1,280);
%--- Make the correlation for whole frequency band (for all freq. bins)
for frqBinIndex = 1:numberOfFrqBins
%--- Generate carrier wave frequency grid (0.5kHz step) -----------
frqBins(frqBinIndex) = settings.IF - ...
(settings.acqSearchBand/2) * 1000 + ...
0.5e3 * (frqBinIndex - 1);
%--- Generate local sine and cosine -------------------------------
sinCarr = sin(frqBins(frqBinIndex) * phasePoints);
cosCarr = cos(frqBins(frqBinIndex) * phasePoints);
%--- "Remove carrier" from the signal -----------------------------
I1 = sinCarr .* signal1;
Q1 = cosCarr .* signal1;
% I2 = sinCarr .* signal2;
% Q2 = cosCarr .* signal2;
% I3 = sinCarr .* signal3;
% Q3 = cosCarr .* signal3;
% I4 = sinCarr .* signal4;
% Q4 = cosCarr .* signal4;
% I5 = sinCarr .* signal5;
% Q5 = cosCarr .* signal5;
% I6 = sinCarr .* signal6;
% Q6 = cosCarr .* signal6;
%--- Convert the baseband signal to frequency domain --------------
IQfreqDom1 = fft(I1 + j*Q1);
% IQfreqDom2 = fft(I2 + j*Q2);
% IQfreqDom3 = fft(I3 + j*Q3);
% IQfreqDom4 = fft(I4 + j*Q4);
% IQfreqDom5 = fft(I5 + j*Q5);
% IQfreqDom6 = fft(I6 + j*Q6);
%--- Multiplication in the frequency domain (correlation in time domain)
convCodeIQ1 = IQfreqDom1 .* caCodeFreqDom;
% convCodeIQ2 = IQfreqDom2 .* caCodeFreqDom;
% convCodeIQ3 = IQfreqDom3 .* caCodeFreqDom;
% convCodeIQ4 = IQfreqDom4 .* caCodeFreqDom;
% convCodeIQ5 = IQfreqDom5 .* caCodeFreqDom;
% convCodeIQ6 = IQfreqDom6 .* caCodeFreqDom;
%--- Perform inverse DFT and store correlation results ------------
acqRes1 = abs(ifft(convCodeIQ1)) .^ 2;
% acqRes2 = abs(ifft(convCodeIQ2)) .^ 2;
% acqRes3 = abs(ifft(convCodeIQ3)) .^ 2;
% acqRes4 = abs(ifft(convCodeIQ4)) .^ 2;
% acqRes5 = abs(ifft(convCodeIQ5)) .^ 2;
% acqRes6 = abs(ifft(convCodeIQ6)) .^ 2;
%--- Check which msec had the greater power and save that, will
%"blend" 1st and 2nd msec but will correct data bit issues
% if (max(acqRes1) > max(acqRes2))
% results(frqBinIndex, :) = acqRes1;
% else
% results(frqBinIndex, :) = acqRes2;
% end
results(frqBinIndex, :) = acqRes1;%+acqRes6;+acqRes2+acqRes3+acqRes4+acqRes5;
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 --------------
[peakSize frequencyBinIndex] = max(max(results, [], 2));
%--- Find code phase of the same correlation peak ---------------------
[peakSize codePhase] = max(max(results));
%--- Find 1 chip wide C/A code phase exclude range around the peak ----
samplesPerCodeChip = round(settings.samplingFreq / settings.codeFreqBasis);
excludeRangeIndex1 = codePhase - samplesPerCodeChip;
excludeRangeIndex2 = codePhase + samplesPerCodeChip;
%--- Correct C/A code phase exclude range if the range includes array
%boundaries
if excludeRangeIndex1 < 1
codePhaseRange = excludeRangeIndex2 : ...
(samplesPerCode + excludeRangeIndex1);
elseif excludeRangeIndex2 > samplesPerCode
codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ...
excludeRangeIndex1;
else
codePhaseRange = [1:excludeRangeIndex1, ...
excludeRangeIndex2 : samplesPerCode];
end
%--- Find the second highest correlation peak in the same freq. bin ---
secondPeakSize = max(results(frequencyBinIndex, codePhaseRange));
%--- Store result -----------------------------------------------------
acqResults.peakMetric(PRN) = peakSize/secondPeakSize;
% If the result is above threshold, then there is a signal ...
if (peakSize/secondPeakSize) > settings.acqThreshold
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fine resolution frequency search ==============================
%--- Indicate PRN number of the detected signal -------------------
fprintf('%02d ', PRN);
%--- Generate 8 msec long C/A codes sequence for given PRN --------
caCode = generateCAcode(PRN);
codeValueIndex = floor((ts * (1:8*samplesPerCode)) / ...
(1/settings.codeFreqBasis));
longCaCode = caCode((rem(codeValueIndex, 1023) + 1));
%--- Remove C/A code modulation from the original signal ----------
% (Using detected C/A code phase)
xCarrier = ...
signal0DC(codePhase:(codePhase + 8*samplesPerCode-1)) ...
.* longCaCode;
%--- Find the next highest power of two and increase by 8x --------
fftNumPts = 8*(2^(nextpow2(length(xCarrier))));
%--- Compute the magnitude of the FFT, find maximum and the
%associated carrier frequency
fftxc = abs(fft(xCarrier, fftNumPts));
uniqFftPts = ceil((fftNumPts + 1) / 2);
[fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5));
fftFreqBins = (0 : uniqFftPts-1) * settings.samplingFreq/fftNumPts;
%--- Save properties of the detected satellite signal -------------
acqResults.carrFreq(PRN) = fftFreqBins(fftMaxIndex);
acqResults.codePhase(PRN) = codePhase;
else
%--- No signal with this PRN --------------------------------------
fprintf('. ');
end % if (peakSize/secondPeakSize) > settings.acqThreshold
end % for PRN = satelliteList
%=== Acquisition is over ==================================================
fprintf(')\n');