GPSbuhuo.zip_GPS 捕获_GPS捕获_GPS捕获仿真_gps定位_捕获

clear; close all; clc; format ('compact'); format ('long', 'g'); addpath include addpath geoFunctions settings = initSettings(); CN0 = 35; Sig_BW = 2e6; SNR = CN0 - 10*log10(Sig_BW); SNR = 10.0^(SNR/10.0); NOISE_DENSITY = -203; Noise_power = Sig_BW*10.0^(NOISE_DENSITY/10); Noise_amp = sqrt(Noise_power); Sig_power = Noise_power*SNR; Sig_amp = sqrt(2*Sig_power); samplesPerCode = round(settings.samplingFreq * (settings.BOCcodeLength/ ... settings.BOCcodeFreq)); signalI = zeros(1,settings.msToProcess*1e-3*settings.samplingFreq); signalQ = zeros(1,settings.msToProcess*1e-3*settings.samplingFreq); ts = 1 / settings.samplingFreq; phasePoints = (1 : (settings.msToProcess*1e-3*settings.samplingFreq)) * 2 * pi * ts; % Visible_sv = [1 3 6 15 20 26]; %可见星 Visible_sv = 1; frqshift = zeros(size(Visible_sv)); frqBins = zeros(size(Visible_sv)); codePhase = zeros(size(Visible_sv)); navdatLength = settings.msToProcess/20; for ii = 1:length(Visible_sv) PRN = Visible_sv(ii); fprintf('%02d ', PRN); frqshift(ii) = settings.acqSearchBand * 1e3 * rand(1,1) - ... (settings.acqSearchBand/2) * 1e3; frqBins(ii) = settings.IF + frqshift(ii); NavData = randi([0,1],1,ceil(navdatLength+1))*2 - 1; dataValueIndex = floor((ts * (1:settings.msToProcess*1e-3*settings.samplingFreq)) / ...(1/50)); NavDataS = NavData(dataValueIndex + 1); clear dataValueIndex caCode1023 = generateCAcode(PRN); codeValueIndex = floor((ts * (1:settings.msToProcess*1e-3*settings.samplingFreq)) / ...(1/settings.codeFreq)); codePhase(ii) = floor(1023*rand(1,1)); codePhase_o=codePhase(ii); longCaCode = caCode1023((rem(codeValueIndex + settings.CAcodeLength ... - codePhase(ii), settings.CAcodeLength) + 1)); clear caCode1023 subCarr = square(settings.subcarrFreq * phasePoints); longBocCode = longCaCode.*subCarr; clear longCaCode subCarr sinCarr = sin(frqBins(ii) * phasePoints); cosCarr = cos(frqBins(ii) * phasePoints); I1 = cosCarr .* longBocCode .* NavDataS; Q1 = sinCarr .* longBocCode .* NavDataS; clear NavDataS longBocCode sinCarr cosCarr signalI = signalI + I1; signalQ = signalQ + Q1; clear I1 Q1 end %% % randn('state',100); Noise = randn(size(signalI))+1i*randn(size(signalI)); Noise = Noise_amp * Noise; signal = Sig_amp*(signalI + 1i*signalQ) + Noise; % signal = signalI + 1i*signalQ; clear signalI signalQ Noise signal_r = real(signal.'); % signal_r = int8(signal_r/max(abs(signal_r))*2^3); signal_r = signal_r/max(abs(signal_r))*2^3; filename = strcat(settings.fileName, '.dat'); fid = fopen(filename,'wb'); count = fwrite(fid,signal_r, settings.dataType); fclose(fid); filename1 = strcat(settings.fileName, '_obsvalue.mat'); save(filename1,'Visible_sv','frqshift','codePhase'); timeScale = 0 : 1/settings.samplingFreq : 5e-3; data = real(signal); longSignal = data; signal1 = longSignal(1 : samplesPerCode); signal2 = longSignal(samplesPerCode+1 : 2*samplesPerCode); signal0DC = longSignal - mean(longSignal); phasePoints = (0 : (samplesPerCode-1)) * 2 * pi * ts; numberOfFrqBins = round(settings.acqSearchBand * 2) + 1; [caCodesTable,bocCodesTable] = makeCaTable(settings); results = zeros(numberOfFrqBins, samplesPerCode); frqBins = zeros(1, numberOfFrqBins); acqResults.carrFreq = zeros(1, 32); acqResults.codePhase = zeros(1, 32); acqResults.peakMetric = zeros(1, 32); fprintf('('); PRNdet=[2 1]; for PRNdetindex=1:2; PRN=PRNdet(PRNdetindex); % PRN=2; bocCodeFreqDom = conj(fft(bocCodesTable(PRN, :))); tic; for frqBinIndex = 1:numberOfFrqBins frqBins(frqBinIndex) = settings.IF - ... (settings.acqSearchBand/2) * 1000 + ... 0.5e3 * (frqBinIndex - 1); sinCarr = sin(frqBins(frqBinIndex) * phasePoints); cosCarr = cos(frqBins(frqBinIndex) * phasePoints); I1 = sinCarr .* signal1; Q1 = cosCarr .* signal1; I2 = sinCarr .* signal2; Q2 = cosCarr .* signal2; %fft IQfreqDom1 = fft(I1 + j*Q1); IQfreqDom2 = fft(I2 + j*Q2); convCodeIQ1 = IQfreqDom1 .* bocCodeFreqDom;%%改caCodeFreqDom; convCodeIQ2 = IQfreqDom2 .* bocCodeFreqDom;%%caCodeFreqDom; acqRes1 = abs(ifft(convCodeIQ1)) .^ 2; acqRes2 = abs(ifft(convCodeIQ2)) .^ 2; if (max(acqRes1) > max(acqRes2)) results(frqBinIndex, :) = acqRes1; else results(frqBinIndex, :) = acqRes2; end end % frqBinIndex = 1:numberOfFrqBins toc; figure(16); % clf(16) subplot(2,2,[3-PRNdetindex 5-PRNdetindex]); phase=0 : (samplesPerCode-1); [X,Y]=meshgrid(phase*settings.CAcodeLength/samplesPerCode,frqBins/1e6); surf(X,Y, results); % title (['第(' int2str(PRN) ')颗可见星的直接捕获结果']); ylabel('载波频率 (MHz)'); xlabel('相位（码片）');zlabel('幅度'); % axis tight; axis([min(min(X)) max(max(X)) min(min(Y)) max(max(Y)) min(min(results)) 1.5e-8]); % hold on end [peakSize, frequencyBinIndex] = max(max(results, [], 2)); [peakSize codePhase] = max(max(results)); figure(19); clf(19) plot((codePhase-37*2:codePhase+37*2)*settings.CAcodeLength/samplesPerCode,... results(frequencyBinIndex, codePhase-37*2:codePhase+37*2)) xlabel('相位（码片）');ylabel('幅度'); samplesPerCACodeChip = round(settings.samplingFreq / settings.codeFreq); excludeRangeIndex1 = codePhase - samplesPerCACodeChip; excludeRangeIndex2 = codePhase + samplesPerCACodeChip; if excludeRangeIndex1 < 1 codePhaseRange = excludeRangeIndex2 : ... (samplesPerCode + excludeRangeIndex1); elseif excludeRangeIndex2 >= samplesPerCode codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ... excludeRangeIndex1; else codePhaseRange = [1:excludeRangeIndex1, ... excludeRangeIndex2 : samplesPerCode]; end secondPeakSize = max(results(frequencyBinIndex, codePhaseRange)); acqResults.peakMetric(PRN) = peakSize/secondPeakSize; if (peakSize/secondPeakSize) > settings.acqThreshold fprintf('%02d ', PRN); caCode1023 = generateCAcode(PRN); codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ... (1/settings.codeFreq)); longCaCode = caCode1023((rem(codeValueIndex, settings.CAcodeLength) + 1)); clear caCode1023 subCarr = square(settings.subcarrFreq * ts * (1:10*samplesPerCode) * 2 * pi); longBocCode = longCaCode.*subCarr; clear longCaCode subCarr xCarrier = ... signal0DC(codePhase:(codePhase + 10*samplesPerCode-1)) ... .* longBocCode;%longCaCode; fftNumPts = 8*(2^(nextpow2(length(xCarrier)))); fftxc = abs(fft(xCarrier, fftNumPts)); uniqFftPts = ceil((fftNumPts + 1) / 2); [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5)); fftFreqBins = (0 : uniqFftPts-1) * settin