www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/acquisition_lj1.m
function acqResults = acquisition_lj1(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.
%% Initialization ================================================
% Find number of samples per spreading code
samplesPerCode = round(settings.samplingFreq / ...
(settings.codeFreqBasis / settings.codeLength));
%--------------------------------------------------------------------------
% 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('(');
%CVS record:
%$Id: acquisition.m,v 1.1.2.12 2006/08/14 12:08:03 dpl Exp $
N=settings.samplingFreq*300e-3;%0.05
N1=settings.samplingFreq*1e-3;
N2=settings.samplingFreq*300e-3;
% cnt=N2/N1;
timebegin=0;
L=2^(nextpow2(N));
L1=2^(nextpow2(N1));
L2=2^(nextpow2(N2));
timeend=N/settings.samplingFreq;
timeend1=N1/settings.samplingFreq;
timeend2=N2/settings.samplingFreq;
snr=-18; % dB
As=10^(snr/20);
zz=fft(longSignal-mean(longSignal));
longSignal=longSignal(1:N);
x2=longSignal.^2-mean(longSignal.^2);
% x2=x1.*conj(x)-mean(x1.*conj(x));
pxx=abs(fft(x2,L));
clear x2
for n=1:settings.numberOfChannels
[fftMax1, fftMaxIndex1] = max(pxx);
% if fftMaxIndex1>L/2
% fftMaxIndex1=L+4-fftMaxIndex1;
% end
% fdm1(n)=(fftMaxIndex1/2-1)/L*settings.samplingFreq;
% fdq1=settings.samplingFreq/2-fdq1;
omiga1=(1-(fftMaxIndex1-1)/L)*2*pi;
fdm1(1)=(omiga1/2)*settings.samplingFreq/(2*pi);
fdm1(2)=(omiga1/2+pi)*settings.samplingFreq/(2*pi);
fftMaxIndex2=L+4-fftMaxIndex1;
omiga1=(1-(fftMaxIndex2-1)/L)*2*pi;
fdm1(3)=(omiga1/2)*settings.samplingFreq/(2*pi);
fdm1(4)=(omiga1/2+pi)*settings.samplingFreq/(2*pi);
index= find(abs(fdm1-settings.IF)<10e3);%|abs(fdm1-(settings.samplingFreq-settings.IF))<10e3
fdm(n)=fdm1(index);
for dot=1:9
pxx(fftMaxIndex1-5+dot)=0;
pxx(L+4-fftMaxIndex1-5+dot)=0;
end
end
% for n=1:4
% [fftMax1, fftMaxIndex1] = max(pxx);
% if fftMaxIndex1>L/2
% pxx(fftMaxIndex1)=0;
% n=n+1;
% else
% fdm(n)=(fftMaxIndex1/2-1)/L*settings.samplingFreq;
% pxx(fftMaxIndex1)=0;
% end
% end
% % fd=sort(fdm,'descend');
% p= fdm~=0;
% fdm=fdm(p);
% clear pxx
t_sig1=timebegin:1/settings.samplingFreq:timebegin+timeend2-1/settings.samplingFreq;
t_sig=timebegin:1/settings.samplingFreq:timebegin+timeend1-1/settings.samplingFreq;
for m=1:32
sg=As*signalgen1(m,timebegin,settings.samplingFreq,timeend2,0,settings.codeFreqBasis);
for n=1:settings.numberOfChannels
s=sg.*exp(-1i*2*pi*fdm(n)*t_sig1);
% s=As*signalgen(m,timebegin,settings.samplingFreq,timeend1,0,fdm(n),settings.codeFreqBasis);
s=fft(s).';
% for h=1:cnt
% sf(:,h)=s(N1*(h-1)+1:N1*h);
% % sf(:,h)=s;
% yf(:,h)=zz(N1*(h-1)+1:N1*h);
% temp(:,h)=abs(fft(conj(sf(:,h)).*yf(:,h),L1));
% end
% temps=sum(temp,2);
% [B,k]=max(temps);
% omiga2=(1-(k-1)/L1)*2*pi;
%
% omiga2=mod(omiga2,settings.samplingFreq*2*pi*1e-3/N1);
%
% omiga(1)=omiga2;
%
% delay(m,n)=omiga(1)*N1/(2*pi*settings.samplingFreq)*1e3;
y1=conj(s).*zz(1:N2);
y2=abs(fft(y1,L2));
[B,k]=max(y2);%(1:uniqFftPt-1));
omiga2=(1-(k-1)/L2)*2*pi;
omiga2=mod(omiga2,settings.samplingFreq*2*pi*1e-3/N2);
omiga(1)=omiga2;
delay(m,n)=omiga(1)*N2/(2*pi*settings.samplingFreq)*1e3;
clear Z yt y1 y2
I1=longSignal(1:samplesPerCode).'.*sin(2*pi*fdm(n)*(t_sig-delay(m,n)*1e-3));
Q1=longSignal(1:samplesPerCode).'.*cos(2*pi*fdm(n)*(t_sig-delay(m,n)*1e-3));
IQ=I1 + 1i*Q1;
% IQfreqDom1 = fft(IQ);
sj2=signalgen1(m,timebegin,settings.samplingFreq,timeend1,delay(m,n)*1e-3,settings.codeFreqBasis);
peakSize=abs(sum(IQ.*sj2));
codePhase=round(delay(m,n)*samplesPerCode);
% SC=IQfreqDom1.*conj(fft(sj2));
% C=abs(ifft(SC)).^2;
% [peakSize,codePhase]=max(C);
%--- 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(C(codePhaseRange));
%--- Store result -----------------------------------------------------
acqResults.peakMetric(m) = peakSize/secondPeakSize;
if (peakSize/secondPeakSize) > settings.acqThreshold
acqResults.carrFreq(m) = fdm(n);
acqResults.codePhase(m) = codePhase;
end
end
end
% [peak,num]=max(C);
% acqResults.peakMetric(num) = peak;
% acqResults.carrFreq(num) = fdm;
% for ii=1:length(num)
% delay1(ii)=delay(num(ii),ii);
% end
% acqResults.codePhase(num) = round(delay1*samplesPerCode);
% %% Initialization ================================================
%
% % Find number of samples per spreading code
%
%
% % 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, :)));
%
% %--- 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+acqRes2+acqRes3+acqRes4+acqRes5;%+acqRes6;
% 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');