LoRa PHY CSS解调

目录

1. LoRa CSS调制回顾

2. LoRa CSS解调原理

3. LoRa CSS解调的一种简单实现

4. 示例与总结

4.1 示例代码

4.1 CSS调制函数示例代码

4.3 标准down-chirp生成示例代码

---

1. LoRa CSS调制回顾

        简要回顾，详细内容可以参考：

LoRa PHY CSS调制分析及Matlab实现_weixin_43270276的博客-CSDN博客_css调制

        LoRa物理层进行信号传输时，通过定义调制带宽BW和扩频因子，得到了码片总数为：

        而不同符号是通过控制不同Chirp信号的初始扫频频率进行区分。扫描起始频率与符号的关系可以定义为：

其中， 为传输symbol值。

        在整个调制符号周期内，基带频率可以表示为：

        或者可以用分段函数的形式表示为：

其中：

        转换到基带信号的相位，可以表示为：

        上式中，假设调制时的初始相位为零，而 为第一段调制结束的相位，作为第二段调制开始时的起始相位，这样可以保证信号相位的连续性： 

        因此LoRa的CSS调制，本质上是将传输带宽BW划分成 个离散频率，然后由每个扫频的初始频率代表需要传递的信息。

        根据《统计信号处理基础》例题3.14的推导分析，频率估计的CRLB下限随采样的点数的三次方进行下降，对数据记录长度最为敏感。因此相对于用幅度和相位传递信息的调制来说，FSK类的调制技术，把需要传递的信息调制到离散频率上，解调时所需要的信噪比更低，因此从本质上能够实现远距离传输。

        而LoRa的CSS调制，除了使用不同的频率表示符号外，还结合Chirp调制技术，将需要传递的信息以不同的初始频率作为扫描起始频率，在整个符号调制时间内，完成带宽范围内的频率扫描，使得LoRa通信具有更好的抗干扰等特性。

---

2. LoRa CSS解调原理

        LoRa CSS解调的核心思路是去除线性扫频分量部分信号后，获取信号的初始频率，最终得到相应的编码值。所以LoRa CSS解调的常规方法之一，就是首先将捕获的采样信号，与标准的down-chirp信号进行相乘运算，消除线性扫频信号，然后对处理后的信号进行FFT运算，以分辨编码后的频率，最后通过频率判断当前的符号值。

        标准的down-chirp信号，基带频率线性变化示意如图 1所示，该信号起始频率由 开始线性下降，最终终止于 ，因此该信号的频率表达式可以表示为：

对应的相位表达式为：

上式中，假设初始相位为0。

 图 1 标准down-chirp信号基带频率线性变化示意图

        理想状态下，符号S的CSS调制信号和down-chirp信号可以分别表示为：

        解调的第一步，就是需要将这两个信号进行相乘操作，即得到：

        由于 采用分段函数表示，因此对分段函数进行分别运算描述。

        假设此时需要传递的符号为S，那么在基带信号中，第一段频率由 扫描至 ，第二段由 扫描至 ，进而完成整个符号的CSS调制频率线性扫描。因此，第一段调制信号与标准down-chirp信号频率变化，如图 2所示。

 图 2 第一段频率线性递增相乘示意

        在 时间内，待解调的信号频率为：

        因此，信号为：

        而 时间内，down-chirp信号频率为

        对应的信号相位为：

        因此，信号为：

        上述两个信号相乘后，在 得到：

        由上式可知，在第一段时间内，调制信号与标准down-chirp信号相乘后，消除了频率线性扫频的效果，最终理想情况下得到了频率为 正弦信号。

 第二段频率线性递增相乘示意

        第二段调制信号与标准down-chirp信号频率变化，如图 3所示。

        在 时间内，带解调的信号频率为：

        对应的信号相位为：

        其中 是上一段调制的最终相位

        因此，信号为：

        而 时间内， down-chirp信号频率为

        对应的信号相位为：

        其中 是上一段调制的最终相位

        因此，信号为：

        上述两个信号相乘后，在 得到：

        由上式可知，在第二段时间内，调制信号与标准down-chirp信号相乘后，消除了频率线性扫频的效果，最终理想情况下得到了频率为 正弦信号。

图 4 与标准down-chirp信号相乘后的LoRa CSS调制波形 

        因此，将LoRa CSS调制后的波形，与标准down-chirp信号相乘，得到的输出波形和信号频率，总结于图 4。由于LoRa传输时限定了传输的BW，因此不同符号由不同起始频率的扫频信号，并在整个符号传输时间内完成整个BW带宽的线性扫频。因此存在需要将超过BW的部分直接取模后减去BW。在利用标准down-chirp信号与调制信号相乘后，直接去除CSS调制时的线性扫频效果，因此在频率取模减去BW的时刻，前后分别为两段正弦信号，且彼此频率上相差BW。显然，如果直接对相乘后的信号做FFT运算，会在-BW至BW这个信号带宽内，出现频谱峰值，该峰值对应的频率，就是LoRa CSS调制对应的起始扫频频率。最终可以通过该频率，计算得到调制的符号，恢复需要传递的信息。

        用示例代码进行demo，将调制符号值设置为20时，其他参数用示例代码中的设置，得到的FFT频谱示意图如图 5所示。此时如果按照示例条件，两个信号的频率应该分别为：156.25kHz和-343.75kHz，与图 5中的信号频率Peak值接近。

图 5 符号值为20时，±BW带宽内信号频谱

        而如果将调制符号值设置为64-20=44时，其他参数用示例代码中的设置，得到的FFT频谱示意图如图 6所示。此时如果按照示例条件，两个信号的频率应该分别为：343.75kHz和-156.25kHz，与图 6中的信号频率Peak值接近。

图 6 符号值为44时，±BW带宽内信号频谱

        比较两个示例，可以发现结果与上述理论计算一致。但从上述图形中也可以发现一些结论：

• ±BW带宽内，存在两个信号峰值，两个峰值相差BW
• 不同符号的扫频起始频率不同，导致两段信号所占据的时间长度不同，根据DFT计算的处理方式，尽管两个信号实际的幅度是一样的，但由于在整个计算周期的时间长度不同，倒是幅度峰值大小不同
• 理论上来说，符号为0至 时，0~BW频带内的信号FFT后的幅度较大，而符号为 至 ，-BW~0频带内的信号FFT后的幅度较大
• 由于采集到的两个频率信号，不是完整周期的，因此DFT运算时存在频谱泄露

 3. LoRa CSS解调的一种简单实现

        下面介绍的一种简单实现方法，基于FFT原理，相关FFT介绍，可以参考：

深入理解离散傅里叶变换(DFT) - 知乎

快速傅里叶变换(FFT)中为什么要“补零”？ - 知乎

快速理解FFT算法（完整无废话） - 知乎

        假设对LoRa CSS调制解调数字信号处理时，每个Code可以采集 个点，由于整个符号有 个，因此整个符号可以采集到的信号长度为：

        由于整个符号传输的时间为：

        因此，可以得到离散化该LoRa CSS信号时的采样率为：

        根据DFT的性质， 个点的时域数据，DFT之后仍然得到 个点的频域数据，此时每个点的频谱分辨率为：

        通过比较发现，该频谱分辨率，正好就是LoRa CSS调制时，不同码片对应扫频起始频率的递增值。也就是说，当确定码片的符号为S时，对应的初始频率可以表示为：

        正是由于这个巧合，可以得到一种基于FFT的简单LoRa CSS解调方法。如果我们将图 4中与标准down-chirp信号相乘后的信号，直接做 个点的FFT，Matlab中输出数组元素下标与对应元素幅度值的关系，如图 7所示。

        根据 个点（ 是偶数时）FFT的性质，可以得到，数组下标为 元素分别存放频率为 的信号信息，也就是 的频率信号。而数组下标为 的元素分别存放频率为 的信号信息,也就是 的信号信息。

图 7 N点（偶数）数据DFT之后数组下标与频率对应关系 

         由于我们仅需要关注±BW内的信息，而FFT之后每个频点对应的频谱分辨率为 ，因此符号为0，1，…, 的扫频初始频率，分别为 ， ，…, ,而这些信号根据此信号FFT的设定，分别在数值下标为1，2,…, 的 个数组中。而相应的负频率， ,  ,…, ，也就是 ，…, 这 符号对应的初始频率，则分别在数值下标为 , ，…, 的 个数组中，如图 8所示。因此，理想的LoRa CSS波形解调，仅需要关注上述数组序号中的信号信息即可，其他数组元素中的信息，可以忽略。

图 8 ±BW带宽内每个数组下标对应的频率 

         同时，由于DFT性质的原因，使得前 个数组中保留正频率信号信息，而后 个数组中保留的是负频率信号信息，因此可以通过Matlab中的类似fftshift函数，将负频率信息搬移至左边，同时忽略超过±BW频率的信号分量，可以得到类似图 9的结果。

图 9 fftshift之后的频谱图 

 因此，当 时，也就是一个码片至少有2个采样点时，我们仅需要在FFT后的数组中，分别提取出下标为1，2,…, 和 , ，…, 个元素对应的幅度值即可，其他频率信息可以直接丢弃。而上述数组元素的下标值，直接和符号S映射关系。

        例如，符号S为0时，此时初始扫频频率为0，对应上述FFT之后的数组，理论上的频率峰值应该出现在数组下标为1的元素中，其他数组元素应该都是噪声，因此数组元素1存在峰值信号。而根据LoRa CSS解调原理，不为零的符号S，初始频率分布为 和 。此时理论上会在序号为 和 对应数组元素中，因此上述两个数组元素理论上存在峰值，而其他数组元素包含的是噪声信号。

        因此，如果我们对FFT之后的数组，搜索峰值对应的下标值，就可以得到当前符号S的可能值。但由于S不为零时，存在两个相差BW的频率信号，因此会有两个峰值。根据上述分析，这两个峰值还会因为S值不同，有时出现在正频率数组下标中，有时出现在负频率数组下标中。而两个频率都相差了BW频率，因此为了简化操作，我们直接将负频率信号向右平移BW个单位后，两组数据相加，或者说对应数组下标的信号幅度直接相加，如图 10所示。最后在叠加后的信号频谱中进行幅度峰值搜索，对应的数组元素下标，就是待传递信号的符号S+1，进而得到解调后的符号。

 图 10 负频率信号向右平移BW后，直接与正频率信号幅度叠加

4. 示例与总结

最后给出了一种LoRa CSS解调的示例代码，以供参考。仿真显示，该方法可以在SNR=-20dB的时候，仍然能够恢复调制的符号信息。

4.1 示例代码

close all; clear all; clc;
BW=500e3;                                                                   %CSS调制带宽
SF=6;                                                                       %CSS调制SF因子
samp_per_code=20;                                                           %CSS调制每个扩频code时间段内的采样点数
symbol_value=60                                                             %当前待调试的符号值
N=2^SF;                                                                     %code 总数
T_symbol=N/BW                                                               %symbol占用时间
SampleRate=samp_per_code*BW                                                 %实际采样率
Npts=samp_per_code*N                                                        %总采样点数
t=(0:Npts)/SampleRate;                                                      %CSS调制抽样时刻，需要的采样时刻序号应该为0：Npts-1，最后一个用于验证相位连续性
init_freq=symbol_value/N*BW;                                                %初始频率
k=BW/T_symbol;                                                              %频率增加斜率

[s0,SampleRate]=LoRa_CSSModDemo(symbol_value, SF, BW, samp_per_code);       %CSS调制
s0=awgn(s0,20,'measured');                                                  %添加高斯白噪声   

freqin1=(symbol_value/N*BW)/1e3                                             %符号对应起始频率，单位kHz
freqin2=(symbol_value/N*BW-BW)/1e3                                          %符号对应起始频率减带宽，单位kHz
freq_res=SampleRate/Npts                                                    %fft频谱分辨率
BW_npts=BW/freq_res                                                         %频率为BW时对应的采样点数,数值上应该等于2^SF
[s1,SampleRate]=LoRa_DownChirpDemo(SF, BW, samp_per_code);                  %创建标准的down-chirp信号

ft=fft(s0.*s1)/Npts;                                                        %相乘后做FFT运算
fvalue=[-BW_npts+1:BW_npts-1]*freq_res/1e3;                                 %频率刻度，从-BW+fr至BW-fr
fftvalue=[ft(end-BW_npts+2:end) ft(1:BW_npts)];                             %根据FFT性质，抽取上述频率对应数组下标
figure('Name', 'CSS FFT Freq','NumberTitle', 'off')
plot(fvalue,abs(fftvalue))                                                  %作图

Codevalue=[-BW_npts+1:BW_npts-1];                                           %频率刻度，从-BW+fr至BW-fr
figure('Name', 'CSS FFT Code','NumberTitle', 'off')                         %频率幅度与码片对应关系
plot(Codevalue,abs(fftvalue))  
ft1=abs(ft(1:BW_npts));                                                     %获取正频率BW内信号分量
ft2=abs(ft(end-BW_npts+2:end));                                             %获取负频率BW内信号分量
ft2=[0 ft2];                                                                %补充直流分量，使得两个数组长度像等
ft_add=ft1+ft2;                                                             %直接幅度相加
[a b]=max(ft_add);                                                          %搜索最大值
symbol_value_DeMod=b-1

4.1 CSS调制函数示例代码

function [s,SampleRate]=LoRa_CSSModDemo(S_Value, SF, BW, samp_per_code)
%S_Value: 待调制的符号，取值0~2^SF-1
%SF;扩频因子
%BW：调制带宽
%samp_per_code：每个code的采样点数
NCode=2^SF;                                                                 %code 总数
T_symbol=NCode/BW;                                                          %symbol占用时间
SampleRate=samp_per_code*BW;                                                %实际采样率
Npts=samp_per_code*NCode;                                                   %总采样点数
init_freq=S_Value/NCode*BW;                                                 %初始频率
k=BW/T_symbol;                                                              %频率增加斜率
Npts1=samp_per_code*(NCode-S_Value);                                        %t1时间段对应总采样点数


t1=(0:Npts1)/SampleRate;                                                    %t1时间段实际相位采样时刻，需要的采样时刻序号应该为0：Npts1-1，最后一个点为第二段的起始相位
tmp=(init_freq-BW/2)*t1+1/2*k*t1.*t1;                                       %根据第一段相位公式计算
Theta1=tmp(1:end-1);                                                        %第一段相位取值
Theta_init=tmp(end);                                                        %第一段相位在t1时刻的取值
t2=(Npts1:Npts)/SampleRate;                                                 %第二段相位采样时刻
t2=t2-t1(end);                              
tmp=1/2*k*t2.*t2-BW/2*t2+Theta_init;                                        %根据第二段相位公式计算
Theta_All=[Theta1 tmp(1:end-1)];
s=exp(j*2*pi*Theta_All);

end

4.3 标准down-chirp生成示例代码

function [s,SampleRate]=LoRa_DownChirpDemo(SF, BW, samp_per_code)
%产生从BW/2至-BW/2的线性扫频
%SF;扩频因子
%BW：调制带宽
%samp_per_code：每个code的采样点数
N=2^SF;                                     %code 总数
T_symbol=N/BW;                              %symbol占用时间
SampleRate=samp_per_code*BW;                %实际采样率
k=BW/T_symbol;                              %频率增加斜率
Npts=samp_per_code*N;                       %总采样点数
t=(0:Npts-1)/SampleRate;
tmp=BW/2*t-1/2*k*t.*t;        %根据第二段相位公式计算
s=exp(j*2*pi*tmp);

end