GPS卫星信号接收机并行捕获程序--针对真实数据（附真实数据集）

【实例简介】

【实例截图】

【核心代码】

%% Correlate signals 关联信号=====   该捕获方法是并行码相位搜索=================================================   
    %--- Perform DFT of C/A code ------------------------------------------
    caCodeFreqDom = conj(fft(caCodesTable(PRN, :)));                          %conj为求复共轭的函数，fft对信号进行快速Fourier变换

    %--- Make the correlation for whole frequency band (for all freq. bins)   对整个频段进行相关（适用于所有频段）
    for frqBinIndex = 1:numberOfFrqBins

        %--- Generate carrier wave frequency grid (0.5kHz step) -----------生成载波频率网格（0.5kHz步长）
        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;

        %--- Convert the baseband signal to frequency domain -------将基带信号转换为频域
        IQfreqDom1 = fft(I1 j*Q1);
        IQfreqDom2 = fft(I2 j*Q2);

        %--- Multiplication in the frequency domain (correlation in time domain)    频域乘法（时域相关性）
        convCodeIQ1 = IQfreqDom1 .* caCodeFreqDom;
        convCodeIQ2 = IQfreqDom2 .* caCodeFreqDom;

        %--- Perform inverse DFT and store correlation results -----执行逆DFT并存储相关结果----   傅里叶反变换---
        acqRes1 = abs(ifft(convCodeIQ1)) .^ 2;          %一旦接收机通过傅里叶反变换计算得到相关值，那么接下来的信号检测与线性搜索捕获法一样，找出所有自相关渎职的峰值
        acqRes2 = abs(ifft(convCodeIQ2)) .^ 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  
%    end  
  frqBinIndex = 1:numberOfFrqBins
 
%  t=toc;
    
%% 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 --在峰值周围找到1个芯片宽的C / A码相位排除范围--
    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   如果范围包括数组边界，则纠正C / A代码阶段排除范围
    if excludeRangeIndex1 < 2
        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 ----peakMetric:峰度测量-------------------------------------------------
    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 ------指示检测到的信号的prn号-------------
        fprintf('%02d ', PRN);                             %%d就是输出整型，用%02d，打印时显示至少两位，前面补零
        
        %--- Generate 10msec long C/A codes sequence for given PRN ----为给定的PRN生成10毫秒长的C / A码序列----
        caCode = generateCAcode(PRN);
        
        codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
                               (1/settings.codeFreqBasis));
                           
        longCaCode = caCode((rem(codeValueIndex, 1023) 1));
    
%！！      %--- Remove C/A code modulation from the original signal --------从原始信号中删除C / A码调制--
        % (Using detected C/A code phase)
        xCarrier = ...
            signal0DC(codePhase:(codePhase 10*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                                     计算FFT的幅度，找到最大值和相关的载波频率
        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 ==================================================