Gold码扩频序列串行捕获MATLAB代码

clear all;
clc;
%发射端
%设置是否有信号
TxCodeOnOff = 1;%1 - On, 0 - Off
%设置仿真时间
TxTotalTime = 1e-3;%1ms
%设置伪码频率
TxCodeFrequency = 1.023e6;%1.023MHz
%设置载波频率
TxCarrierFrequency = 2 * 1.023e6;%2.046MHz
%设置采样频率
TxSampleFrequency = 16 * 1.023e6;%MHz
%发射端产生伪码并计算码长
TxGoldSequence0 = dyc_gold_gen_v3();
TxGoldSequence0Len = length(TxGoldSequence0);
%用载波调制
TxTemp1 = TxSampleFrequency / TxCodeFrequency;
%将伪码序列扩展
TxGoldSequence1 = repmat(TxGoldSequence0, TxTemp1, 1);
TxGoldSequence1 = reshape(TxGoldSequence1, 1, TxTemp1 * TxGoldSequence0Len);
TxGoldSequence2 = TxCodeOnOff * TxGoldSequence1;
TxTiming = 0 : (1 / TxSampleFrequency) : (TxTotalTime - 1 / TxSampleFrequency);
TxSignalSequence = sin(2 * pi * TxCarrierFrequency * TxTiming + TxGoldSequence2 * pi);

% figure(11);%伪码序列
% plot(TxGoldSequence1(3450:3620));
% axis tight;
% figure(12);%调制后发射的信号
% plot(TxSignalSequence(3450:3620));
% axis tight;

%信道
%设置多普勒频移
ChDiffFrequency = -4500;%Hz
%设置信噪比
ChSN = -20;%dB
%设置码相移
ChDelayValue = 800;%码片相移量
ChDelay = ChDelayValue * TxTemp1;%此值为采样结果的码相移，是伪码相移的 TxTemp1 倍

%将序列进行移位（需要先移位在加多普勒频移，否则的话移位操作会影响多普勒频移）
ChGoldSequence = circshift(TxGoldSequence2, [0, ChDelay]);
%生成带有多普勒频移的序列
ChTiming = TxTiming;
ChSignalSequence1 = sin(2 * pi * (TxCarrierFrequency + ChDiffFrequency) * ChTiming + ChGoldSequence * pi);
%加入噪声
ChSingalSequence3 = awgn(ChSignalSequence1, ChSN, 'measured');%(此函数还可添加更多参数)

% figure(21);%产生码相移的信号
% plot(ChGoldSequence(3450:3620));
% figure(22);%产生频移的信号
% plot(ChSignalSequence1(2500:4500));
% axis tight;
% axis tight;
% figure(23);%加入噪声的信号
% plot(ChSingalSequence3(3450:3620));
% axis tight;

%接收端
%设置频率搜索步长
RxCarrierFrequencyDelta = 500;%Hz
%设置频率搜索范围
RxDiffFrequencyMax = 20 * RxCarrierFrequencyDelta;%Hz
%设置码相位搜索步长
RxDelayDelta = 10;
%设采样率和码速率之比与发送端一致
RxDelayRatio = TxTemp1;
%设接收端采样率与发送端一致
RxTiming = TxTiming;

%
%计算接收端频率搜索次数
RxFrequencySearchTimes = fix(2 * RxDiffFrequencyMax / RxCarrierFrequencyDelta) + 1;%(可能有误)
%计算码相位搜索次数
RxGoldLen = TxGoldSequence0Len;
RxGoldSearchTimes = ceil(RxGoldLen / RxDelayDelta);
%
RxSignalAmplitude = zeros(RxGo

ldSearchTimes, RxFrequencySearchTimes);

%
% %测试1：初始化不经过调制解调时接收端码相干结果
% RxTestRes1 = zeros(1, RxGoldSearchTimes);


%
for RxTemp11 = 1 : RxGoldSearchTimes

%
display(['Processing... ' num2str(RxTemp11 / RxGoldSearchTimes * 100, '%.2f') '%']);

RxTemp12 = RxTemp11 - 1;
%产生接收伪码
RxGoldSequence1 = TxGoldSequence1;
RxGoldSequence2 = circshift(RxGoldSequence1, [0, RxTemp12 * RxDelayDelta * RxDelayRatio]);%移动码相位
RxGoldSequence3 = RxGoldSequence2 .* 2 - 1;%

% %测试1：计算不经过调制解调时接收端码相干
% RxTestRes1(RxTemp11) = TxGoldSequence2 * RxGoldSequence3';


%计算码相关
RxSignalSequence1 = ChSingalSequence3 .* RxGoldSequence3;


%
RxCarrierFrequency1 = TxCarrierFrequency - RxDiffFrequencyMax;
for RxTemp21 = 1 : RxFrequencySearchTimes
RxTemp22 = RxTemp21 - 1;
%产生接收载波
RxCarrierFrequency = RxCarrierFrequency1 + RxTemp22 * RxCarrierFrequencyDelta;
RxCarrierSequenceI = cos(2 * pi * RxCarrierFrequency * RxTiming);
RxCarrierSequenceQ = sin(2 * pi * RxCarrierFrequency * RxTiming);
%I/Q两路分别与接收载波做相关
RxSignalXcorr2I = RxSignalSequence1 * RxCarrierSequenceI';
RxSignalXcorr2Q = RxSignalSequence1 * RxCarrierSequenceQ';
%幅度检波(平方相加再开根号)
RxSignalAmplitude(RxTemp11, RxTemp21) = sqrt(RxSignalXcorr2I ^ 2 + RxSignalXcorr2Q ^ 2);
end
end


% figure(32); %测试1：显示不经过调制解调时接收端码相干结果
% plot(RxTestRes1);

%生成坐标轴
%纵坐标（行）
RxLabelCode = ((1 : RxGoldSearchTimes) - 1) * RxDelayDelta;
%横坐标（列）
RxLabelFrequency = ((1 : RxFrequencySearchTimes) - 1) * RxCarrierFrequencyDelta - RxDiffFrequencyMax;
%生成三维捕捉结果
figure(31);
mesh(RxLabelFrequency, RxLabelCode, RxSignalAmplitude);


RxSignalAmplitude1 = RxSignalAmplitude(:);%(将结果转化为一维向量便于后续计算)
%计算捕捉到的最大幅度
RxSignalResMax = max(RxSignalAmplitude1);
%计算捕捉结果的平均幅度
RxSignalResAverage = mean(RxSignalAmplitude1);
%
RxSignalResDeviation = std(RxSignalAmplitude1);



%判断是否捕获
%设置门限
RxCatchThreshold = RxSignalResAverage + 5 * RxSignalResDeviation;
%
if RxSignalResMax > RxCatchThreshold
%捕获成功
display('Success!!!');
%查找最大值对应的坐标
[RxResRow RxResColumn] = find(RxSignalAmplitude == RxSignalResMax);%(此方法仅

限二位数组)
display(['Doppler Frequency: ' num2str(RxLabelFrequency(RxResColumn), '%.0f') 'Hz']);
display(['Code Phase: ' num2str(RxLabelCode(RxResRow), '%.0f')]);
else
%捕获失败
display('Fail!!!');
end