GPS软件接收机（1）——GPS信号捕获

 samplesPerCode = round(settings.samplingFreq / ...

    
                   (settings.codeFreqBasis / settings.codeLength));
    
 % 计算每个C/A码包含多少采样点，“...”表示续行
    
  
    
 signal1 = longSignal(1 : samplesPerCode);
    
 signal2 = longSignal(samplesPerCode+1 : 2*samplesPerCode);
    
 %取相邻的两个信号，信号长度为1ms，采样点数为samplesPerCode
    
  
    
 signal0DC = longSignal - mean(longSignal); 
    
 %去直流
    
  
    
 ts = 1 / settings.samplingFreq;
    
 %计算采样时间间隔
    
  
    
 phasePoints = (0 : (samplesPerCode-1)) * 2 * pi * ts;
    
 %计算1ms信号内的载波相位点
    
  
    
 numberOfFrqBins = round(settings.acqSearchBand * 2) + 1;
    
 %计算频率捕获的个数，500Hz步进 总共2kHz，共41个
    
  
    
 caCodesTable = makeCaTable(settings);
    
 % makeCaTable函数根据采样频率产生32颗指定卫星的CA码 ，此时为32*samplesPerCode矩阵
    
  
    
  
    
 %--- 初始化捕获结果矩阵 -------------------------------
    
  
    
 results     = zeros(numberOfFrqBins, samplesPerCode);
    
 %一颗卫星的包含搜索频率个数与一个C/A码采样个数的零向量
    
  
    
 % Carrier frequencies of the frequency bins
    
 % 
    
 frqBins     = zeros(1, numberOfFrqBins);
    
 %产生一颗卫星的载波偏移初始向量
    
  
    
 %--- Initialize acqResults ------------------------------------------------
    
 acqResults.carrFreq     = zeros(1, 32);
    
 %初始化32颗卫星的载波频率
    
  
    
 acqResults.codePhase    = zeros(1, 32);
    
 % C/A code phases of detected signals
    
  
    
 acqResults.peakMetric   = zeros(1, 32);
    
 % Correlation peak ratios of the detected signals
    
  
    
 fprintf('(');

 for PRN = settings.acqSatelliteList

    
 %settings.acqSatelliteList=1:32
    
  
    
     caCodeFreqDom = conj(fft(caCodesTable(PRN, :)));%某颗PRN卫星的CA码采样值的FFT区共轭 
    
  
    
     for frqBinIndex = 1:numberOfFrqBins
    
     %带宽为500Hz的频率遍历
    
    
    
     frqBins(frqBinIndex) = settings.IF - ...
    
                            (settings.acqSearchBand/2) * 1000 + ...
    
                            0.5e3 * (frqBinIndex - 1);
    
     %载波偏移--10kHz—+10kHz 产生numberOfFrqBins个500Hz间隔的中频载波频率 
    
  
    
     
    
     sinCarr = sin(frqBins(frqBinIndex) * phasePoints);
    
     cosCarr = cos(frqBins(frqBinIndex) * phasePoints);
    
  %----产生一个信号周期完整（1ms）的中频载波正弦和余弦信号-------------------------
    
  
    
  
    
     %---------- "Remove carrier" from the signal-----------------------
    
     %-----------------连续的两个1ms中频数据载波剥离----------------------
    
     I1      = sinCarr .* signal1;
    
     Q1      = cosCarr .* signal1;
    
     I2      = sinCarr .* signal2;
    
     Q2      = cosCarr .* signal2;
    
  
    
     %--- Convert the baseband signal to frequency domain ------
    
     %--------载波剥离后的1ms数据做FFT---------------------------
    
     IQfreqDom1 = fft(I1 + j*Q1);
    
     IQfreqDom2 = fft(I2 + j*Q2);
    
  
    
     convCodeIQ1 = IQfreqDom1 .* caCodeFreqDom;%改为转置
    
     convCodeIQ2 = IQfreqDom2 .* caCodeFreqDom;%改为转置
    
    %频域共轭相乘，时域卷积判断最大值
    
  
    
     acqRes1 = abs(ifft(convCodeIQ1)) .^ 2;
    
     acqRes2 = abs(ifft(convCodeIQ2)) .^ 2;
    
    %傅里叶反变化获得C/A码的相关值
    
  
    
     %--- Check which msec had the greater power and save that, will
    
     %"blend" 1st and 2nd msec but will correct data bit issues
    
     %寻找卷积相关最大值，复现码的相位值 完成每个频率域的最大值搜索 
    
     %连续的两个1ms数据为了防止数据位的翻转对于相关值的影响
    
     %形成frqBinIndex*samplesPerCode的矩阵
    
     if (max(acqRes1) > max(acqRes2))
    
         results(frqBinIndex, :) = acqRes1; %将获得的相关值赋给results矩阵的第
    
     else
    
         results(frqBinIndex, :) = acqRes2;
    
     end
    
     
    
     end

 %% 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
    
     %将第一主峰与第二主峰进行比较，当第二主峰远离最高主峰超过一个码片的话，则作为第二主峰    
    
    
    
     %--- 选择最大相关值与载波频率 --------------
    
     [peakSize frequencyBinIndex] = max(max(results, [], 2));
    
     % 获取最大相关值与最大相关值所对应的频率搜索序号
    
  
    
     %--- Find code phase of the same correlation peak ---------------------
    
     [peakSize codePhase] = max(max(results));
    
     %获取最大相关值与最大相关值所对应的码相位序号
    
  
    
     %--- 获取与最大峰超过一个码片的次最大峰 ----
    
     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 < 2
    
     codePhaseRange = excludeRangeIndex2 : ...
    
                      (samplesPerCode + excludeRangeIndex1);
    
                      
    
     elseif excludeRangeIndex2 >= samplesPerCode
    
     codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ...
    
                      excludeRangeIndex1;
    
     else
    
     codePhaseRange = [1:excludeRangeIndex1, ...
    
                       excludeRangeIndex2 : samplesPerCode];
    
     end
    
 %确定寻找次最大峰的搜索区域
    
    
    
     secondPeakSize = max(results(frequencyBinIndex, codePhaseRange));
    
    %在当前已经搜索到的多普勒频率下的codePhaseRange范围内寻找次大值
    
  
    
     %--- Store result -----------------------------------------------------
    
     acqResults.peakMetric(PRN) = peakSize/secondPeakSize;
    
     %存储相关值最大值与次大值的比值 每一颗卫星均做记录
    
     
    
   
    
     if (peakSize/secondPeakSize) > settings.acqThreshold
    
   % 如果该结果大于设定的阈值，则确认捕获到信号
    
     
    
     fprintf('%02d ', PRN);
    
     %输出捕获到的卫星编号
    
     
    
     caCode = generateCAcode(PRN);%确定已经捕获到的CA码
    
     
    
 %生成10ms长的捕获C/A码序列，并用10ms的C/A码对载波信号进行剥离
    
     codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
    
                            (1/settings.codeFreqBasis));
    
     %10ms内可采样的码片数量
    
                        
    
     longCaCode = caCode((rem(codeValueIndex, 1023) + 1));
    
     %模拟信号采样，使得本地C/A数据格式与捕获信号保持一致
    
  
    
     xCarrier = ...
    
         signal0DC(codePhase:(codePhase + 10*samplesPerCode-1)) ...
    
         .* longCaCode;
    
     %用10ms的码相位对捕获信号进行码剥离
    
  
    
     %--- Find the next highest power of two and increase by 8x --------
    
     fftNumPts = 8*(2^(nextpow2(length(xCarrier))));
    
      %将点数转化为2的幂次整数以便使用fft（否则为DFT，运算时间增加），且将FFT的频率分辨率提高8倍
    
     
    
     fftxc = abs(fft(xCarrier, fftNumPts)); 
    
 %fftNumPts超过xCarrier长度，补0会提高频率分辨率，见解释：
    
     
    
     uniqFftPts = ceil((fftNumPts + 1) / 2);
    
 %由于频率谱的对称性，选择一半频谱即可估计载波频率值
    
     
    
 [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5)); %返回最大频率值与相对应的序列值，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');
    
 save( 'acqResults');

    caCode = generateCAcode(PRN);%确定已经捕获到的CA码

    
     
    
 %生成10ms长的捕获C/A码序列，并用10ms的C/A码对载波信号进行剥离
    
     codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
    
                            (1/settings.codeFreqBasis));
    
     %10ms内可采样的码片数量
    
                        
    
     longCaCode = caCode((rem(codeValueIndex, 1023) + 1));
    
     %模拟信号采样，使得本地C/A数据格式与捕获信号保持一致
    
  
    
     xCarrier = ...
    
         signal0DC(codePhase:(codePhase + 10*samplesPerCode-1)) ...
    
         .* longCaCode;
    
     %用10ms的码相位对捕获信号进行码剥离
    
  
    
     %--- Find the next highest power of two and increase by 8x --------
    
     fftNumPts = 8*(2^(nextpow2(length(xCarrier))));
    
      %将点数转化为2的幂次整数以便使用fft（否则为DFT，运算时间增加），且将FFT的频率分辨率提高8倍
    
     
    
     fftxc = abs(fft(xCarrier, fftNumPts)); 
    
 %fftNumPts超过xCarrier长度，补0会提高频率分辨率，见解释：
    
     
    
     uniqFftPts = ceil((fftNumPts + 1) / 2);
    
 %由于频率谱的对称性，选择一半频谱即可估计载波频率值
    
     
    
 [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5)); %返回最大频率值与相对应的序列值，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');
    
 save( 'acqResults');