268 lines
12 KiB
Mathematica
268 lines
12 KiB
Mathematica
|
function acqResults = acquisition(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 (số mẫu trong mỗi mã trải phổ)
|
||
|
|
samplesPerCode = round(settings.samplingFreq / ...
|
||
|
|
(settings.codeFreqBasis / settings.codeLength));
|
||
|
|
|
||
|
|
% Create two 1msec vectors of data to correlate with and one with zero DC
|
||
|
|
% (tạo ra 2 vector 1ms để tương quan và 1 vector có DC = 0)
|
||
|
|
signal1 = longSignal(1 : samplesPerCode);
|
||
|
|
signal2 = longSignal(samplesPerCode+1 : 2*samplesPerCode);
|
||
|
|
% Tạo ra 2 signal để
|
||
|
|
|
||
|
|
signal0DC = longSignal - mean(longSignal);
|
||
|
|
|
||
|
|
% Find sampling period: khoảng thời gian lấy mẫu
|
||
|
|
ts = 1 / settings.samplingFreq;
|
||
|
|
|
||
|
|
% Find phase points of the local carrier wave: điểm pha của sóng mang cục
|
||
|
|
% bộ
|
||
|
|
phasePoints = (0 : (samplesPerCode-1)) * 2 * pi * ts;
|
||
|
|
|
||
|
|
% Number of the frequency bins for the given acquisition band (500Hz steps)
|
||
|
|
% số lượng ngăn tần số cho băng tần acquisition
|
||
|
|
numberOfFrqBins = round(settings.acqSearchBand * 2) + 1;
|
||
|
|
|
||
|
|
% Generate all C/A codes and sample them according to the sampling freq.
|
||
|
|
% Tạo các mã C/A và lấy mẫu theo tần số lấy mẫu
|
||
|
|
caCodesTable = makeCaTable(settings);
|
||
|
|
|
||
|
|
|
||
|
|
%--- Initialize arrays to speed up the code -------------------------------
|
||
|
|
% Khởi tạo mảng tăng tốc mã
|
||
|
|
% Search results of all frequency bins and code shifts (for one satellite)
|
||
|
|
% Kết quả tìm kiếm của tất cả ngăn tần số và dịch chuyển mã
|
||
|
|
results = zeros(numberOfFrqBins, samplesPerCode);
|
||
|
|
|
||
|
|
% Carrier frequencies of the frequency bins
|
||
|
|
% Tần số sóng mang của ngăn tần số
|
||
|
|
frqBins = zeros(1, numberOfFrqBins);
|
||
|
|
|
||
|
|
|
||
|
|
%--- Initialize acqResults ------------------------------------------------
|
||
|
|
% Khởi tạo acqResults
|
||
|
|
% Carrier frequencies of detected signals: tần số sóng mang của tín hiệu
|
||
|
|
% được phát hiện
|
||
|
|
acqResults.carrFreq = zeros(1, 32);
|
||
|
|
% C/A code phases of detected signals: pha mã C/A của tín hiệu được phát
|
||
|
|
% hiện
|
||
|
|
acqResults.codePhase = zeros(1, 32);
|
||
|
|
% Correlation peak ratios of the detected signals: tỉ lệ cực đại tương quan
|
||
|
|
% của các tín hiệu được phát hiện
|
||
|
|
acqResults.peakMetric = zeros(1, 32);
|
||
|
|
|
||
|
|
fprintf('(');
|
||
|
|
|
||
|
|
% Perform search for all listed PRN numbers ... : thực hiện tìm kiếm các
|
||
|
|
% PRN được liệt kê
|
||
|
|
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)
|
||
|
|
% Tạo mối tương quan cho toàn bộ dải tần
|
||
|
|
for frqBinIndex = 1:numberOfFrqBins
|
||
|
|
|
||
|
|
%--- Generate carrier wave frequency grid (0.5kHz step) -----------
|
||
|
|
% Tạo lưới tần số sóng mang (0.5 kHz)
|
||
|
|
frqBins(frqBinIndex) = settings.IF - ...
|
||
|
|
(settings.acqSearchBand/2) * 1000 + ...
|
||
|
|
0.5e3 * (frqBinIndex - 1);
|
||
|
|
|
||
|
|
%--- c -------------------------------
|
||
|
|
% Tạo sin và cos cục bộ
|
||
|
|
sinCarr = sin(frqBins(frqBinIndex) * phasePoints);
|
||
|
|
cosCarr = cos(frqBins(frqBinIndex) * phasePoints);
|
||
|
|
|
||
|
|
%--- "Remove carrier" from the signal -----------------------------
|
||
|
|
% Loại bỏ sóng mang
|
||
|
|
IQ1 = (sinCarr+1i*cosCarr).* signal1;
|
||
|
|
I1 = real(IQ1);
|
||
|
|
Q1 = imag(IQ1);
|
||
|
|
|
||
|
|
IQ2 = (sinCarr+1i*cosCarr).* signal2;
|
||
|
|
I2 = real(IQ2);
|
||
|
|
Q2 = imag(IQ2);
|
||
|
|
|
||
|
|
%--- Convert the baseband signal to frequency domain --------------
|
||
|
|
% Chuyển đổi tín hiệu băng cơ sở sang miền tần số
|
||
|
|
IQfreqDom1 = fft(I1 + 1i*Q1);
|
||
|
|
IQfreqDom2 = fft(I2 + 1i*Q2);
|
||
|
|
|
||
|
|
%--- Multiplication in the frequency domain (correlation in time
|
||
|
|
%domain)
|
||
|
|
% Phép nhân trong miền tần số (tương quan trong miền thời gian)
|
||
|
|
convCodeIQ1 = IQfreqDom1 .* caCodeFreqDom;
|
||
|
|
convCodeIQ2 = IQfreqDom2 .* caCodeFreqDom;
|
||
|
|
|
||
|
|
%--- Perform inverse DFT and store correlation results ------------
|
||
|
|
% Thực hiện DFT nghịch đảo và lưu kết quả tương quan
|
||
|
|
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
|
||
|
|
% Kiểm tra ms nào mạnh hơn và lưu lại, trộn ms thứ 1 và 2 nhưng sẽ
|
||
|
|
% khắc phục vấn đề về bit dữ liệu
|
||
|
|
if (max(acqRes1) > max(acqRes2))
|
||
|
|
results(frqBinIndex, :) = acqRes1; %Lưu acqRes1 vào hàng thứ frqBinIndex
|
||
|
|
else
|
||
|
|
results(frqBinIndex, :) = acqRes2;
|
||
|
|
end
|
||
|
|
|
||
|
|
end % frqBinIndex = 1:numberOfFrqBins
|
||
|
|
|
||
|
|
%% Look for correlation peaks in the results ==============================
|
||
|
|
% Find the highest peak and compare it to the second highest peak
|
||
|
|
% Tìm đỉnh cao nhất và so sánh với đỉnh cao thứ 2
|
||
|
|
% The second peak is chosen not closer than 1 chip to the highest peak
|
||
|
|
% Đỉnh thứ 2 được chọn không gần hơn 1 chip so với đỉnh cao nhất
|
||
|
|
|
||
|
|
%--- Find the correlation peak and the carrier frequency --------------
|
||
|
|
% Tìm đỉnh tương quan và tần số sóng mang
|
||
|
|
[peakSize frequencyBinIndex] = max(max(results, [], 2)); %Giá trị lớn nhất trên mỗi hàng
|
||
|
|
|
||
|
|
%--- Find code phase of the same correlation peak ---------------------
|
||
|
|
% Tìm pha mã của cùng 1 đỉnh tương quan
|
||
|
|
[peakSize codePhase] = max(max(results));
|
||
|
|
|
||
|
|
%--- Find 1 chip wide C/A code phase exclude range around the peak ----
|
||
|
|
% Tìm 1 chip mã C/A loại trừ khoảng xung quanh đỉnh
|
||
|
|
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
|
||
|
|
% Đúng vị trí mã C/A và loại bỏ khoảng nếu khoảng đó bao gồm các biên
|
||
|
|
% của mảng
|
||
|
|
if excludeRangeIndex1 < 2
|
||
|
|
codePhaseRange = excludeRangeIndex2 : ...
|
||
|
|
samplesPerCode;%(samplesPerCode + excludeRangeIndex1);
|
||
|
|
|
||
|
|
elseif excludeRangeIndex2 >= samplesPerCode
|
||
|
|
codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ...
|
||
|
|
excludeRangeIndex1;
|
||
|
|
else
|
||
|
|
codePhaseRange = [1:excludeRangeIndex1, ...
|
||
|
|
excludeRangeIndex2 : samplesPerCode];
|
||
|
|
end
|
||
|
|
try
|
||
|
|
%--- Find the second highest correlation peak in the same freq. bin ---
|
||
|
|
% Tìm đỉnh tương quan cao thứ 2 trong cùng freq. bin
|
||
|
|
secondPeakSize = max(results(frequencyBinIndex, codePhaseRange));
|
||
|
|
catch exception
|
||
|
|
msgbox(exception.message);
|
||
|
|
end;
|
||
|
|
%--- Store result -----------------------------------------------------
|
||
|
|
acqResults.peakMetric(PRN) = peakSize/secondPeakSize;
|
||
|
|
|
||
|
|
% If the result is above threshold, then there is a signal ...
|
||
|
|
% Nếu kết quả vượt ngưỡng thì có tín hiệu
|
||
|
|
if (peakSize/secondPeakSize) > settings.acqThreshold
|
||
|
|
|
||
|
|
freq_Range=1:numberOfFrqBins;
|
||
|
|
xR=settings.IF - (settings.acqSearchBand/2) * 1000 + 0.5e3 * (freq_Range - 1);
|
||
|
|
yR=[0:samplesPerCode-1];
|
||
|
|
figure;
|
||
|
|
surf(results);shading interp, title("PRN "+PRN);
|
||
|
|
|
||
|
|
%% Fine resolution frequency search =======================================
|
||
|
|
|
||
|
|
%--- Indicate PRN number of the detected signal -------------------
|
||
|
|
% Cho biết số PRN của tín hiệu được phát hiện
|
||
|
|
fprintf('%02d ', PRN);
|
||
|
|
|
||
|
|
%--- Generate 10msec long C/A codes sequence for given PRN --------
|
||
|
|
% Tạo chuỗi mã C/A dài 10 ms cho PRN đã cho
|
||
|
|
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 ----------
|
||
|
|
% Loại bỏ điều chế mã C/A khỏi tín hiệu gốc
|
||
|
|
% (Using detected C/A code phase) (Sử dụng pha mã C/A được phát
|
||
|
|
% hiện)
|
||
|
|
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
|
||
|
|
% Tính toán độ lớn của FFT, tìm tần số sóng mang tối đa và liên
|
||
|
|
% quan
|
||
|
|
fftxc = abs(fft(xCarrier, fftNumPts));
|
||
|
|
|
||
|
|
uniqFftPts = ceil((fftNumPts + 1) / 2);
|
||
|
|
[fftMax, fftMaxIndex] = max(fftxc);
|
||
|
|
fftFreqBins = (0 : fftNumPts) * settings.samplingFreq/fftNumPts;
|
||
|
|
|
||
|
|
|
||
|
|
%--- Save properties of the detected satellite signal -------------
|
||
|
|
% Lưu thuộc tính tín hiệu vệ tinh được phát hiện
|
||
|
|
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');
|