Chirp扩频信号原理验证及抗干扰性能仿真

1 Chirp扩频收发仿真

1.1 Chirp扩频信号产生

1.1.1 Chirp信号生成

LoRa 使用 CSS （Chirp Spread Spectrum）线性扩频调制。基本通信单元是线性调频信号（Linear Frequency Modulation, LFM），即频率随时间线性增加（或减小）的信号，也被称为Chirp信号。频率随着时间线性增加的chirp符号称为upchirp，将频率随着时间线性减小的chirp符号称为downchirp。最高频率和最低频率之间的差值为LoRa的带宽B。
设basic upchirp的最低频率为$f_0=-B/2$，最高频率为 $f_1=B/2$，chirp时间长度为T。其频率可以表示为
$$f(t)=f_0+kt$$
其中$k=B/T$为扫频速度。线性变化的频率对时间做积分可得相位
$$\varphi \left(t\right)=2\pi \left(f_{0}t+\frac{1}{2}k{t}^2\right)$$
basic upchirp表示为
$$C(t)=e^{j2\pi (f_0+\frac{1}{2}kt)t}$$
basic downchirp表示为
$$C^{\ast }(t)=e^{j2\pi (f_1-\frac{1}{2}kt)t}$$
设SF = 7，带宽B = 10 Hz，采样频率Fs = 100 Hz，根据表达式画出basic upchirp和basic downchirp的时域波形及时频图如图 3 1所示。从时频图可以直观看出频率随时间线性变化的关系，对于upchirp，频率从$f_0=-B/2$到$f_1=B/2$，对于downchirp，频率从$f_1=B/2$到$f_0=-B/2$。从时域波形直观来看，upchirp和downchirp在T/2时刻频率取0，另外两者的实部时域波形一样，虚部时域波形互为相反数，与其定义表达式一致。

1.1.2 CSS调制

CSS调制基于频率移位。当嵌入数据时，首先令basic upchirp乘上一个固定频率的偏移分量，偏移后的信号可以表示为
$$C(t)e^{j2\pi \Delta ft}$$
将所有频率高于$f_1$的信号段循环频移至$f_0$频率处。如果定义了$M=2^{SF}$种不同的偏移频率，最多可以表示SF比特的数据。SF即扩频因子（Spreading Factor，），定义为
$$2^{SF}=BT$$
符号速率为$R_s=1/T$，比特速率为$R_b=SF/T$。给定带宽B，SF越大，每一个Chirp的长度T越长，相应符号速率越小。SF用于调节传输速率和接收灵敏度，越大的SF传输速率越小，但支持更远的通讯距离。
考虑将复包络中心频率置为0，将频率加上-B/2的偏移，得到1个符号周期内调制符号 m ∈ { 0 , 1 , … , M − 1 } m \in \{ 0,1, \ldots ,M - 1\} m∈{0,1,…,M−1}对应LoRa信号瞬时频率为
$$f(t;m)=-\frac{B}{2}+m\frac{B}{M}+kt-Bu(t-{\tau }_m),0\le t\le T_s$$
其中$u(t)$为阶跃函数，${\tau }_m$为瞬时频率从最大值降为最小值对应的时刻，为
$${\tau }_m=\left(1-\frac{m}{M}\right)T$$
则有符号$m$对应LoRa信号时域表达式
$$c(t;m)=e^{j2\pi \left[\frac{1}{2}k{t}^2+m\frac{B}{M}t-\frac{B}{2}t-Btu(t{-}_{{\tau }_m})\right]}$$

1.2. Chirp扩频信号解调

LoRa解调过程，实质就是求出chirp符号的起始频率，共分两步。
（1）解扩：将接收信号与downchirp信号相乘，得单频信号，频率等于编码chirp的频率偏移量：
$$C^{\ast }(t)\cdot C(t)e^{j2\pi \Delta ft}=e^{j2\pi \Delta ft}$$
将发送信号表达式带入，有
$$\Delta f=m\frac{B}{M}-Bu(t{-}_{{\tau }_m})=\{\begin{array}{rl}& m\frac{B}{M},t<{\tau }_m\\ & \left(\frac{m}{M}-1\right)B,t\ge {\tau }_m\end{array}$$
也即解扩后信号有2个频点。
（2）解调：对相乘信号做FFT将时域信号转化为频域，得到单边频谱绝对值最大值对应下标，即对应信号编码的数据，即
$$\hat{m}=\mathrm{mod}\left(\underset{k=0,1,...,M-1}{\mathrm{argmax}}(|Y_k|),M\right)$$
设带宽为10Hz，采样频率为100Hz，扩频因子为7。如图 3 3，设传输符号为30，则频率起始为$-B/2+B\times 30/2^7=-2.6562\mathrm{Hz}$ ，在$t=\left(1-30/128\right)T=9.7997$s处频率切换回 -5 Hz。解调时首先进行解扩，也即用downchirp进行相乘，得到由两个频率，频率1为$mB/M=2.3438Hz$，频率2为$-B/2+mB/M=-2.6562Hz$，映射到符号，频率1对应为30，频率2对应为 -98，取FFT后最大值对应的频率1对应符号，对128取模，得到恢复的传输符号为30。
需要注意的是，该方法没有利用到负频率信息，因此进行改进，由于正负频率波峰之间相差$-B/2$的带宽，因此将相差$-B/2$ 的频谱相加，之后再求最大值，相当于利用了全部的波峰能量。

2 Chirp物理层信号结构

2.1 LoRa帧结构

前导码：长度2~65535，全部为basic upchirp符号。
帧同步：长度2，调制符号分别为 { x , 2 S F − x } \{ x,{2^{SF}} - x\} {x,2SF−x}，定义所在网络符号。接收机回忽略不属于自己网络的帧。
频率同步：长度2.25，为downchirp符号，最后0.25个符号用于时间和频率同步。
头部：非必须，包含有效负载长度、前向纠错码率、CRC。
有效负载：传输的信息数据。
CRC：非必须，循环校验码。

实验：带宽500 KHz，采样频率5 MHz，扩频因子为7。前导码长度5，帧同步信息为{44,84}，不加载头部信息、CRC，有效负载长度为5。负载信息为[37,92,48,113,85]。经过解扩，恢复出频点映射到符号为[37,92,48,113,85]。

2.2 LoRa物理层模型

信道编码 ：编码类型为Hamming码，提高抵抗噪声的能力。
交织：采用对角交织，将码字中比特分散到不同符号中，减少突发错误带来的影响。
编码交织将传输的二进制比特流划分为信息比特块，每个块的大小为$SF\times 4$，编码后对每一行的4个比特增加CR个校验位，得到$SF\times (4+CR)$大小比特块。
进制转换：对二进制比特块进行进制转化（格雷映射），映射到十进制符号。

3 Chirp扩频抗干扰性能仿真

这里仅对CSS调制进行抗干扰性能分析，不考虑信道编码和交织。

3.1 单音干扰下抗干扰原理验证

产生单音干扰（实信号），接收机接收信号表示为
$$r(t)=s(t)+J(t)$$
仿真参数如表所示。
表 3 1 Chirp扩频调制单音干扰下抗干扰原理验证仿真参数表

带宽	500 KHz
采样率	5 MHz
扩频因子	6
单音干扰频率	125 KHz
信干比	20 dB
发送符号	10

发送信号信息是10。未加干扰时的接收信号频谱和解扩信号频谱如图 3 8所示。在正频率处，出现单音频率0.078125MHz，映射到符号上即为10。

加入当加入20 dB干信比、频率为125 KHz的单音干扰后，经过解扩，单音干扰的能量扩散到频谱上，频率峰值仍为0.078125MHz，能够恢复出传输符号。

3.2 单音干扰下抗干扰性能仿真

固定单音干扰频率为125 KHz，在干信比范围为-10到21 dB进行仿真，扩频因子分别设置为6，7，8。不考虑噪声，得到不同干信比下的误码率如图所示。

固定干信比，得到AWGN信道下误码率随比特信噪比变化的曲线。显然，随着扩频因子增大，抗单音干扰性能提升。

4 代码

4.1 Chirp信号产生

clc;clear;close all;
%% 参数设置
SF = 7;% 10
BW = 10; % Hz 100
Fs = 100; % Hz 50
%% basic chirp 产生
T = 2^SF/BW;        % Chirp长度
t = 0:1/Fs:T-1/Fs;       % 采样时刻
k = BW / T;         % 频率变化率（线性）
fUp = -BW/2 + k*t'; % upchirp 频率变化 -Bw/2~BW/2
phaseUp = 2*pi*(-BW/2 + k/2*t) .* t;     % upchirp 时变相位
fDown = BW/2 - k*t'; % downchirp 频率变化 Bw/2~-Bw/2
phaseDown = 2*pi*(BW/2 - k/2*t) .* t;  % downchirp 时变相位
BasicUpChirp = exp(1j*phaseUp);
BasicDownChirp = exp(1j*phaseDown);

figure;set(gcf,'Color','w','Position',[0 0 500 400]);
subplot(221);
plot(t,real(BasicUpChirp),'Color',[25,79,151]/255,'LineWidth',0.6);hold on;
plot(t,imag(BasicUpChirp),'Color',[0,150,107]/255,'LineWidth',0.6);hold off;
xlim([0,T]);ylim([-1.2 1.2]);lgd=legend('real','imag');lgd.NumColumns = 2;
lgd.Location = 'southoutside';lgd.Box = "off";
title('basic upchirp time domain');xlabel('t(s)');
text(-1.5,1,'(a)','FontName','Euclid','FontSize',8);
subplot(223);
plot(t,fUp,'Color',[200,45,49]/255,'LineWidth',0.8);
title('downchirp freq-time gram');xlim([0,T]);xlabel('t(s)');
text(-1.5,5,'(b)','FontName','Euclid','FontSize',8);

subplot(222);
plot(t,real(BasicDownChirp),'Color',[25,79,151]/255,'LineWidth',0.6);hold on;
plot(t,imag(BasicDownChirp),'Color',[0,150,107]/255,'LineWidth',0.6);hold off;
xlim([0,T]);ylim([-1.2 1.2]);lgd=legend('real','imag');lgd.NumColumns =2;
lgd.Location = 'southoutside';lgd.Box = "off";title('downchirp time domain');
xlabel('t(s)');text(-1.5,1,'(c)','FontName','Euclid','FontSize',8);
subplot(224);
plot(t,fDown,'Color',[200,45,49]/255,'LineWidth',0.8);xlim([0,T]);xlabel('t(s)');
text(-1.5,5,'(d)','FontName','Euclid','FontSize',8);title('downchirp freq-time gram');

%% 调制符号
% test on symbols
m_pool = [1,10,20,40,60];
figure(2);set(gcf,'Color','w','Position',[0 0 700 350]);
f1 = figure(2);
ax1 = axes(f1,"Position",[0.1 0.1 0.4 0.8]);
lgdwords = {};
for i = 1:length(m_pool)
    m = m_pool(i); % symbol in (0,2^SF -1)
    M = 2^SF;
    taum = (1-m/M)*T;
    u_taum = zeros(1,length(t));
    index = find(t>=taum);
    u_taum(index(1):end) = 1; % step function
    Dataf = -BW/2 + m/M*BW+ k*t-BW*u_taum; % -BW/2 ~ BW/2 当t = tau时取0
    DataPhase = 2*pi*(-BW/2 + m/M*BW+ k/2*t-BW*u_taum).*t;
    DataChirp = exp(1j*DataPhase); % Chirp时变相位
    plot(ax1,t,Dataf,'LineWidth',0.8);hold on;
    lgdwords = [lgdwords;['m=',num2str(m)]];
    axi = axes(f1,"Position",[0.6 0.1+0.175*(i-1) 0.3 0.1]);
    plot(axi,t,real(DataChirp),'LineWidth',0.5,'Color',[25,79,151]/255);hold on;
    plot(axi,t,imag(DataChirp),'LineWidth',0.5,'Color',[0,150,107]/255);hold off;
    title(axi,['m=',num2str(m)]);    axi.XLim = [0,T];
end
leg1 = legend(ax1,lgdwords);
% leg1.Location = 'southoutside';
leg1.Location = 'northwest';
leg1.NumColumns = 1;leg1.Box = "off";
ax1.YLim = [-1.2*BW/2,1.2*BW/2];
ax1.XLabel.String = '\it{t} (s)';
ax1.YLabel.String = 'Frequency(Hz)';
ax1.FontName = 'Euclid';
set(ax1,'FontSize',8,'Fontname', 'Euclid');
set(ax1,'linewidth',0.8); 
ax1.Title.String = 'Freq-Time Gram';
%% dechirp
m = 30; % symbol in (0,2^SF -1)
M = 2^SF;
taum = (1-m/M)*T;
u_taum = zeros(1,length(t));
index = find(t>=taum);
u_taum(index(1):end) = 1; % step function
Dataf = -BW/2 + m/M*BW+ k*t-BW*u_taum; % -BW/2 ~ BW/2 当t = tau时取0
DataPhase = 2*pi*(-BW/2 + m/M*BW+ k/2*t-BW*u_taum).*t;
DataChirp = exp(1j*DataPhase); % Chirp时变相位
ChirpMult2 = DataChirp.*BasicDownChirp;
ChirpMult2(find(isnan(ChirpMult2))) = 0;
L = length(t);
fshift = (-L/2:L/2-1)*(Fs/L);
CodeShift = (-L/2:L/2-1);
Y = fft(ChirpMult2,L);
Z = fftshift(Y);
P2 = abs(Z/L);
index = find(P2 ==max(P2));
b = CodeShift(index(1));%搜索最大值
sym_dec_DeMod=mod(b,M);% 恢复符号

figure;set(gcf,'Color','w','Position',[0 0 500 500]); 
subplot(511);
plot(t,Dataf,'Color',Mycolor(5),'LineWidth',0.6);
xlim([0,T]);text(-1.5,0,'(a)','FontName','Euclid','FontSize',8);
box on;xlabel('t(s)');
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.8)
subplot(512);
plot(t,real(DataChirp),'Color',[25,79,151]/255,'LineWidth',0.6);hold on;
plot(t,imag(DataChirp),'Color',[98,91,161]/255,'LineWidth',0.6);hold off;
xlim([0,T]);xlabel('t(s)');ylim([-1.5,1.5]);
text(-1.5,0,'(b)','FontName','Euclid','FontSize',8);
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.8)
subplot(513);
plot(t,real(ChirpMult2),'Color',[189,107,8]/255,'LineWidth',0.6);
ylim([-2,2]);xlim([0,T]);box on;xlabel('t(s)');
text(-1.5,0,'(c)','FontName','Euclid','FontSize',8);
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.8)
subplot(514);
semilogy(fshift,P2,'Color',[0,127,84]/255,'LineWidth',0.6);
xlabel('Freq(Hz)');xlim([fshift(1),fshift(end)]);ylim([10^-3,1]);box on;
text(1.25*fshift(1),0.5,'(d)','FontName','Euclid','FontSize',8);
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.8)
subplot(515);
plot(CodeShift,P2,'Color',[0,127,84]/255,'LineWidth',0.6);
xlabel('code');xlim([-M+1,M+1]);box on;
text(1.25*CodeShift(1),0.5,'(e)','FontName','Euclid','FontSize',8);
set(gca,'FontName','Times New Roman','FontSize',10,'LineWidth',0.8)

4.2 Chirp信号组帧

clear;close all;
%% 参数设定
% CSS调制参数
% BW = 10; % Hz 1004
% Fs = 100; % Hz 50
SF = 7;% 扩频因子
BW = 500e3; % 带宽
Fs = 10*BW; % 采样频率
T = 2^SF/BW; % 符号周期
M = 2^SF;
Nsamples = M*Fs/BW; % 采样点数
%组帧参数
Npre = 5; % 前导码个数
Ndata = 5; % 负载符号个数
NetIndex = 44; % 帧同步符号值
NFreqSync = 2.25;
%信道参数
EbN0dB = 100;
SNRdB = EbN0dB - 10*log10(0.5*T*Fs/SF);% 信噪比 /dB
%% 发射部分
% basic chirp 产生
SegUpChirp = UpChirpGen(SF,BW,Fs);
SegDownChirp = DownChirpGen(SF,BW,Fs);
% Preamble chirp前导码产生
SigPreChirp = repmat(SegUpChirp',Npre,1)';
% 帧同步 chirp 产生
SegFrameChirp =  [DataChirpGen(SF,BW,Fs,NetIndex),...
                 DataChirpGen(SF,BW,Fs,M-NetIndex)];
% 频率同步 chirp产生
SegFreSyncChirp = repmat(SegDownChirp',3,1)';
SegFreSyncChirp = [SegFreSyncChirp(1:floor(NFreqSync*Nsamples)),zeros(1,3*Nsamples - floor(NFreqSync*Nsamples))];
% data chirp 产生
symTx = randi(M-1,1,Ndata); % symbol in (0,2^SF -1)
SigDataChirp = zeros(1,Nsamples*Ndata);% 负载信号
for i = 1:Ndata
    SegDataChirp =  DataChirpGen(SF,BW,Fs,symTx(i));
    SigDataChirp((i-1)*Nsamples+1:i*Nsamples) = SegDataChirp;
end
% 发射信号产生
SigTx = [SigPreChirp,SegFrameChirp,SegFreSyncChirp,SigDataChirp];
tFrame = (1:length(SigTx))/Fs;
figure;set(gcf,'Color','w','Position',[100 100 500 125]);
stft(SigTx,Fs,Window=kaiser(Nsamples/8,5),OverlapLength=128,FFTLength=Nsamples);
ylim([-20*BW/2/10e6,20*BW/2/10e6]);xlim([0,tFrame(end)*1e3]);
title('Transmitted Frame Signal');
%% 接收部分
Ps = 1;SNR = 10^(SNRdB/10);Pn = Ps/SNR;
SigRx = SigTx + sqrt(Pn/2) * (randn(size(SigTx))*(1+1j));% add Noise
SigDownChirp = repmat(SegDownChirp',floor(length(SigRx)/Nsamples),1)';
SigDeChirp = SigRx.*SigDownChirp; % DeChirp

figure;set(gcf,'Color','w','Position',[100 100 500 125]);
stft(SigDeChirp,Fs,Window=kaiser(Nsamples/8,5),OverlapLength=128,FFTLength=Nsamples);
ylim([0,20*BW/2/10e6]);xlim([0,tFrame(end)*1e3]);
title('DeChirp Frame Signal');

symRec = zeros(1,Ndata);
for i = 1 : Ndata
    range = (i-1)*Nsamples+1:i*Nsamples;
    symRec(i) = ChirpDeMod(SigDataChirp(range),SegDownChirp,M);
end
%% 解调
function sym_dec_DeMod = ChirpDeMod(DataChirp,DownChirp,M)
    ChirpMult2 = DataChirp.*DownChirp;
    ChirpMult2(find(isnan(ChirpMult2))) = 0;
    L = length(DataChirp);
    CodeShift = (-L/2:L/2-1);
    Y = fft(ChirpMult2,L);
    Z = fftshift(Y);
    P2 = abs(Z/L);
    index=find(P2 ==max(P2));
    b = CodeShift(index(1));%搜索最大值
    sym_dec_DeMod=mod(b,M);
end
%% CSS 调制
function DataChirp = DataChirpGen(SF,BW,Fs,m)
    T = 2^SF/BW;        % Chirp长度
    t = 0:1/Fs:T-1/Fs;   % 采样时刻
    k = BW / T;         % 频率变化率（线性）
    M = 2^SF;
    taum = (1-m/M)*T;
    u_taum = zeros(1,length(t));
    index = find(t>=taum);
    u_taum(index(1):end) = 1; % step function
    DataPhase = 2*pi*(-BW/2 + m/M*BW+ k/2*t-BW*u_taum).*t;
    DataChirp = exp(1j*DataPhase); % Chirp时变相位
end
%% Basic UpChirp
function UpChirp = UpChirpGen(SF,BW,Fs)
    T = 2^SF/BW;        % Chirp长度
    t = 0:1/Fs:T-1/Fs;   % 采样时刻
    k = BW / T;         % 频率变化率（线性）
    phaseDown = 2*pi*(-BW/2 + k/2*t) .* t;  % downchirp 时变相位
    UpChirp = exp(1j*phaseDown);
end
%% Basic DownChirp
function DownChirp = DownChirpGen(SF,BW,Fs)
    T = 2^SF/BW;        % Chirp长度
    t = 0:1/Fs:T-1/Fs;   % 采样时刻
    k = BW / T;         % 频率变化率（线性）
    phaseDown = 2*pi*(BW/2 - k/2*t) .* t;  % downchirp 时变相位
    DownChirp = exp(1j*phaseDown);
end