用于低成本接收机的LoRa SF11 500KHz波形检测解调算法

前一篇里，获取了LORAwan的物理层波形，并通过Octave查看了它的瞬时频率。LoRa是私有协议，网上已经有了很不错的开源的实现，如：

S2_LoRa通信实验

LoRaPhy

以及GNU-Radio的Lora模块、LimeSDR的Lora实现。当我试图修改上述代码实现上一篇文章里获取的波形的还原时，却无法得到正确的结果，也找不到原因。为了知其所以然，还是准备从头动手实现。由于是业余学习，不打算解决所有LORA协议的接收处理问题，只就事论事，只做默认出厂配置SF-11 500K一类。

接收机是国产B205-mini或者B-210，开发环境是开源的计算工具 GNU-Octave 。本文代码参考：
a4_lora/octave

0 .实验学习结论摘要

虽然是站在LoRa-PHY、GNU-Radio巨人的肩膀上，但是在近1个月的不断踩坑、尝试中，还是非常的不顺利。最终得到正确结果时，感觉就像攀过了一座陡峭的小山。

1. 我们发现对于Lora这种使用Chirp的起始频率携带信息的低成本硬件，其固有的时钟差、频率差对结果影响极大。尤其是使用另一个不靠谱的硬件来接收（比如我的山寨 B205mini），两者的钟差是不靠谱+不靠谱=离谱。能在Low SNR下解析出正确的数据，需要对波形的深刻理解和认识，每一个步骤都值得思考。
2. 对我学习的清华大学项目 LoRa-PHY里，用于同步时钟误差的方法，是假设主要的误差都来自LoRa设备，用频差去折算钟差（时差）。在本实验里，我们的接收设备也很不靠谱，使得这个举措变得不稳定。我们通过观察相关峰值的最大值的周期变化，发现了隐含在头部14组up-chirp里的低频周期，提出了基于相关峰起伏周期的钟差估计算法，并使用它对钟差进行补偿，效果很好。
3. 对收发时钟、频率稳定度都很低的情形，常见的二维时频相关算法耗时太大。提出基于 up-down chirp 边缘峰值交汇的时频差精确估计算法，只在窗口的边缘检测峰值，并两两连线，计算交汇坐标。

本人已经退休，不需要发Paper。谁如果论文里想引用这些idea，只要标明参考出处即可。

1. 波形样本来源

样本来源就是上一节使用山寨B205mini/B210采集的IQ路数据。由于靠的很近，几乎没有噪声。采样率是1MHz，双字节 int16。LoRa的设备的带宽是500KHz，SR=11。

波形样本已经跟随代码，一起放在Git服务器上。

2. 粗略显示LoRa物理层突发的结构

对于无噪声的波形，可以采用相位差分的方式，直接观察瞬时频率。读取一个突发数据，
并利用复数求取角度差分，就是相位差，也就是瞬时频率。

上图表明，我们购买的Lora模块默认开启了14个首部up-chirp，2.25个down-chirp

%%%%%%%%%%%%%%%读取一个2倍速率采集的无噪声文件。
%% USRP B210 430MHz, 16bit IQ BW500KHz
%fid = fopen('d:/lora/test.wav',"rb");
fid = fopen('data/iq_sf11_500KHz_spr1MHz.wav',"rb");
raw_iq_1d = fread(fid,'int16')';
fclose(fid);
%%组合复波形
raw_iq = raw_iq_1d(1:2:end)+1j*raw_iq_1d(2:2:end);
%幅度归一化
max_amp = sum(abs(raw_iq))/length(raw_iq);
raw_iq = raw_iq/max_amp;figure(1);
hold on
draw_freqs = zeros(1,length(raw_iq));
draw_freqs(2:end) = angle(raw_iq(2:end)) - angle(raw_iq(1:end-1));
draw_freqs(draw_freqs<=-pi ) = draw_freqs(draw_freqs<=-pi ) + 2*pi ;
draw_freqs(draw_freqs>=pi ) = draw_freqs(draw_freqs>=pi ) - 2*pi ;
plot(draw_freqs,'r');

2.1 chirp的数学描述

根据LoRa的定义，一个chirp 跨域的带宽是 B Hz，持续的时间是 T 秒，他们满足:

$$BT=2^{SF}$$

在T秒内，基带IQ向量旋转的频率从 -B/2 变化到 B/2。

由于频率

$$f(t)=2\pi (f_0+\frac{B}{T}t)$$

而向量的相位是频率的积分，计算：

$$p(t)=\int f(t)dt=2\pi (f_0+\frac{B}{2T}t)t,t\subset [0,T)$$

此时，样点的复平面向量可以表示为：

$$s(t)=e^{2\pi j(f_0+\frac{B}{2T}t)t}$$

这个t的取值是 0 ~ T秒。

2.2 chirp的归一化离散形式

使用上述模拟波形的样式来计算chirp, 不容易出错。但如果换算为离散形式，则更方便写程序。假设我们采用高于LoRa带宽B的K倍速率采样，即 Sr = BK，样点为 n，则参考本文开头的第一个式子里SF、B、T的关系，一个chirp的样点数为：
$$SYM=N(SF,K)=\frac{2^{SF}}{B}Sr=2^{SF}\cdot K$$
比如SF=11, K=2时，N就是4096点。根据离散、连续的关系，

$$t=\frac{n}{Sr}=\frac{n}{BK}$$

$$p(t)=2\pi (f_0+\frac{B}{2T}\frac{n}{BK})\frac{n}{BK}=2\pi (f_0+\frac{n}{2T\cdot K})\frac{n}{BK}$$

由于 $BT=2^{SF}$, 有：

$$p(t)=2\pi (\frac{f_0}{B}+\frac{n}{2\cdot 2^{SF}\cdot K})\frac{n}{K}$$

剩下一个 f0/B，也可以兑换为比例，其实f0取-B/2时，他就是 -1/2。

经过上面的推导，我们把 chirp 表示为和具体采样率、频率无关的值，只和SF,K有关。

$$p_{up-chirp}(t)=2\pi (-\frac{1}{2}+\frac{n}{2\cdot 2^{SF}\cdot K})\frac{n}{K}$$

$$S_{up-chirp}(n)=e^{2\pi j(-\frac{1}{2}+\frac{n}{2\cdot 2^{SF}\cdot K})\frac{n}{K}}$$

我们可以使用Octave生成这样的chirp:

chirp_n  =0:1:SYM-1;
chirp_up    = exp(2j*pi*(-1/2+chirp_n/2./(2^SF)./K).*chirp_n./K);
plot(real(chirp_up));

3. 在频差条件下检测头部

使用与上面的 chirp 对应的相关滤波器，即可从噪声中检测头部的存在。也可以利用down-chirp 相乘后，fft查看有无单音来检测，那样速度更快。

检测这个步骤，使用共轭，能够和原始波形直接合成为很强的峰。

$$\overline{S_{up-chirp}}(n)=e^{-2\pi j(-\frac{1}{2}+\frac{n}{2\cdot 2^{SF}\cdot K})\frac{n}{K}}$$

这个波形的瞬时频率为：

$$\overline{f_{up-chirp}}(n)=-2\pi (-\frac{1}{2}+\frac{n}{2^{SF}\cdot K})\frac{1}{K}$$

我们在第一阶段，主要的工作就是用上面这个相关波形去和原始波形做相关，两两相乘而后求和，检测峰值。

3.1 钟差频差对头部检测几乎没有影响

即使在有频偏的情况下，这种检测依然是有效的。只要频率偏移不是很离谱，通过样点的移动，总能用时间延迟的损失来抵消频偏：

设存在一个时间偏移 $n^{\prime }=\delta +n$, 以及一个固有的频率偏移 $\Delta$：

$$\overline{f_{up-chirp}}(n^{\prime })=\overline{f_{up-chirp}}(\delta +n)=-2\pi j(\Delta -\frac{1}{2}+\frac{n+\delta }{2\cdot 2^{SF}\cdot K})\frac{1}{K}$$

当这个$\delta$满足下式，使得频偏$\Delta$其恰好能被时间推移抵消：

$$\Delta +\frac{\delta }{2\cdot 2^{SF}\cdot K}=0$$

则结果变为

$$\overline{f_{up-chirp}}(n^{\prime })=-2\pi (-\frac{1}{2}+\frac{n^{\prime }-\delta }{2\cdot 2^{SF}\cdot K})\frac{1}{K}$$

相关峰与无频偏时相比，只是提前或者推迟了$\delta$个样点，幅度稍微有些降低而已（可用的重叠频段不再是B Hz了）。

3.2 检测结果

我们的头部检测结果如下：

绿色的是瞬时频率，蓝色的是相关峰，黑色的是峰值标记位置。请格外注意这些蓝色峰值的起伏，它很重要！

4. 基于相关峰起伏周期的钟差估计算法

在第3节，我们发现即使不去补偿时频差，也能顺利检测到头部。但是，LoRa靠频率的绝对起点f0来表示符号, 解扩、解调不做钟差、频差，就得不到正确的结果。

4.1 低成本接收机的困境

学习资源 LoRa-PHY里，认为钟差基本源自LoRa设备本身。这或许因为大学实力雄厚，采集设备的时钟稳定而标准。此时，可以沿用LoRa-PHY的算法，使用频差/射频频率来兑换钟差。

然而，我们的接收机成本很低，本身的频差、钟差就很大，和LoRa这种60多元的板子旗鼓相当。此时，上述算法就不行了。LoRa-PHY抵消后存在一个微弱的残差，使得在发送64字节数据时，符号依旧会漂移3-4次。

4.2 相关峰起伏周期的成因

注意看，3.2节的头部检测相关峰，每个峰值的幅度并不一致！ 讲起来，按照Lora的标准定义，采用K倍速率采样，不会有起伏。即使具有频偏，相关位置的相位也是定值，每次相关的幅度衰减是恒定的。为啥周期起伏？还是个COS函数的样子？那是因为时钟存在误差，每次对齐的相位会发生漂移。

这个小细节蕴含着非常重要的一个参数，正是LoRa设备的时钟和我们的接收机时钟之间的误差。这个误差，使得原本应该 N 点1个符号的chirp，接收到后可能变成 N+0.12点，如果不去补偿这个误差，以N为周期去卡符号，就会慢慢错位。错位后，每次相关时，相位初始值发生缓慢平移，造成峰值周期性的起伏。

提取这个误差，非常简单，对这十几个点，做FFT，如对 14个点做1048576点的FFT，并求取峰值，峰值与窗口的比例偏移就是时钟偏移。我们把上图的顶部放大：

可以看到这就是个余弦。用上面的数据做FFT，求取补偿因子 cv_pay :

   vpk = sqrt(var(prepeaks (:,2)));prepeaks (:,2)  = (prepeaks (:,2) - mpk)/vpk * 100;ck_fft = prepeaks (:,2);cm_fft = abs(fft(ck_fft,1048576));[ci_fft,cv_fft] = max(cm_fft(1:1048576/2));cv_pay = cv_fft / 1048576/SYM;

cv_pay 的意思就是每个样点要补偿多少点，是个很小的数。太大了就不对了。

4.3 计算结果

本例子里是 3.9700e-05，1个符号产生 3.9700e-05 * N = 0.1626点偏移，在K=2时，每6个chirp就会冒一个点，每12个chip就会发生1次FFT峰值的移动，造成错误。用LoRaPHY的方法，计算出的是0.127，二者还是存在相当的差异。

5 基于 up-down chirp 边缘峰值交汇的时频差精确估计算法

由于存在第三部分的原因，无论是单独采用up chirp，还是 down chirp，都不能准确测定频差。但是，Lora存在一个up-down chirp的符号反转，构成了所谓的双线性调频(dual-chirp)，就用这个波形进行匹配滤波。

5.1 笨办法观察规律

最笨的办法，就是两层 for 循环，把原始波形进行时间搬移、频率搬移，而后求取各个搬移参数下的相关峰。

下图显示的是搜索的窗口。根据特定的频差，首先把载波乘以 $e^{2\pi j\Delta n}$,而后，再平移$\delta$样点，计算共轭相乘的和。

在窗口内，遍历所有的$\Delta ,\delta$,
$$D(m,n)=e^{2\pi j(-\frac{1}{2}+\Delta +\frac{m+\delta }{2\cdot 2^{SF}\cdot K})\frac{m+\delta }{K}}\cdot e^{-2\pi j(-\frac{1}{2}+\frac{n}{2\cdot 2^{SF}\cdot K})\frac{n}{K}}$$

while(pos < end_updownpos)x_ti = pos - start_updownpos + 1;for (x_fi = 1:searchf)sumsig = raw_iq(pos + chirp_n2).*exp(1j.*2*pi*(x_fi-searchf/2)/2./(2^SF)./K./K.*chirp_n2) .*  chirp_updown(1+chirp_n2);sumv = abs(sum(sumsig));updown_tf(x_ti,x_fi) = sumv;endforpos = pos + 1;end %while

结果如下：

这个计算极其耗时，实用性不高。不过，这个图却很好的说明了up-down chirp做时频二维相关的特点。

• 中间的尖峰：时频都对准了，up,down同时相关上了，产生双倍的扩频增益。
• 两条相交的山脊：由于时频可以互补，发生了第三章阐述的平移，部分相关上1条（up或者down），另一条没有增益。
• 山脊是线性的、相交于尖峰。

5.2 基于边缘最大值直线交点快速交汇的时频差估计

上面的这个算法，计算实在太慢了。注意到其实upchip，downchip各自的最大值在上述搜索窗口的边缘能够准确捕获，只要在边缘4个边上进行计算，而后通过连线交点估计位置即可。

   if (pos < end_updownpos)x_ti = pos - start_updownpos + 1;for (x_fi = 1:searchf)if (pos>start_updownpos && pos+1<end_updownpos)if (x_fi>1 && x_fi<searchf)continue;endifendifsumsig = raw_iq(pos + chirp_n2).*exp(1j.*2*pi*(x_fi-searchf/2)/2./(2^SF)./K./K.*chirp_n2) .*  chirp_updown(1+chirp_n2);sumv = abs(sum(sumsig));updown_tf(x_ti,x_fi) = sumv;endforpos = pos + 1;elsex_f = zeros(1,4);x_t = zeros(1,4);max_ct = 0;raw_tf = updown_tf;%找到4个峰值while (max_ct < 4)[max_xtv,max_xti] = max(updown_tf);[max_xfv,max_xfi] = max(max(updown_tf));curr_f = max_xfi;curr_t = max_xti(curr_f);updown_tf(curr_t,curr_f) = 0;check_pk = 1;fake_peak = 0;while (check_pk <= max_ct && fake_peak==0)distance = sqrt((x_f(check_pk) - curr_f)^2 +(x_t(check_pk) - curr_t)^2);if (distance < 16 * K)fake_peak = 1;endifcheck_pk = check_pk + 1;endwhileif (fake_peak==1)continue;endifmax_ct = max_ct + 1;x_f(max_ct) = curr_f;x_t(max_ct) = curr_t;endwhile%求直线交点xrge = size(updown_tf);[x_it1,y_it1] = intersection_point(x_t(1),x_f(1),x_t(2),x_f(2),x_t(3),x_f(3),x_t(4),x_f(4));[x_it2,y_it2] = intersection_point(x_t(1),x_f(1),x_t(3),x_f(3),x_t(2),x_f(2),x_t(4),x_f(4));[x_it3,y_it3] = intersection_point(x_t(1),x_f(1),x_t(4),x_f(4),x_t(2),x_f(2),x_t(3),x_f(3));if (x_it1 > 0 && x_it1 <xrge(1) && y_it1 > 0 && y_it1 <xrge(2))x_tc = floor(x_it1+0.5);x_fc = floor(y_it1+0.5);x_tf = x_it1;x_ff = y_it1;elseif (x_it2 > 0 && x_it2 <xrge(1) && y_it2 > 0 && y_it2 <xrge(2))x_tc = floor(x_it2+0.5);x_fc = floor(y_it2+0.5);x_tf = x_it2;x_ff = y_it2;elseif (x_it3 > 0 && x_it3 <xrge(1) && y_it3 > 0 && y_it3 <xrge(2))x_tc = floor(x_it3+0.5);x_fc = floor(y_it3+0.5);x_tf = x_it3;x_ff = y_it3;elsex_tc = xrge(1)/2;x_fc = xrge(2)/2;x_tf = x_tc;x_ff = x_fc;endifprintf('Finished.\n');if (x_ff - searchf/2 >0)payback_clk = cv_pay;elsepayback_clk = -cv_pay;endif%figure(2);pos = start_updownpos + x_tc - 1;syn_pos = start_updownpos + x_tc - 1;updown_tf = raw_tf;sumsig = raw_iq(pos + chirp_n2).*exp(1j.*2*pi*(x_ff-searchf/2)/2./(2^SF)./K./K.*chirp_n2) .*  chirp_updown(1+chirp_n2);sumv = abs(sum(sumsig));updown_tf(x_tc,x_fc) = sumv;mesh(updown_tf);endif

结果：

需要注意的是，第四步里估计的钟差是一个绝对值，要根据本方法提供的频差的正负号进行调整。当然，如果使用的设备时钟很精确，则可以直接使用 LoRaPHY使用的方法。

6 卡住位置补偿解调

一旦找到了精确的时间、频率，则可以对剩余部分进行解调。解调的方法就是先乘以 down-chip，做fft、折半到 B带宽后，读取峰值的位置。

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3. 解调if (status==3)draw_mark(pos) = -30;%搬移频率到0中心data = raw_iq(pos:pos+SYM-1).*exp(1j.*2*pi*(deltaf-searchf/2)/2./(2^SF)./K./K.*chirp_n).*chirp_dw_dem;fft_data_raw = abs(fftshift(fft(data,SYM)));code_len = SYM/K;if (K>1)wrp_N = SYM;wrp_P = wrp_N / 2;wrp_M = 2^(SF-1);fft_data_raw(wrp_P-wrp_M+1:wrp_P)= fft_data_raw(wrp_P-wrp_M+1:wrp_P) + fft_data_raw(wrp_P+wrp_M+1:wrp_P+2*wrp_M);fft_data_raw(wrp_P+1:wrp_P+wrp_M)= fft_data_raw(wrp_P+1:wrp_P+wrp_M) + fft_data_raw(wrp_P-2*wrp_M+1:wrp_P-wrp_M);fft_data1 = fft_data_raw(wrp_P-wrp_M+1:wrp_P+wrp_M);elsefft_data1 = fft_data_raw;endif[codev,code_d] = max(fft_data1);fft_data = fft_data1;draw_hdseq(pos:pos+code_len-1)=fft_data/K;figure(1);clf;hold on;plot(draw_mark,'k');plot(draw_hdseq*30/SYM,'b');plot(draw_freqs,'g');hold off;pause(0.001);raw_demsym(1,codi) = mod(code_d+2^(SF-1),2^SF);codi = codi + 1;%样点偏移补偿(timeing offset),这里存疑，因为 payback_clk 是靠频偏估计的clk_pay = clk_pay + SYM;PZ = payback_clk * codi * SYM;pos = round(pos_start + clk_pay + PZ);endif

输出:

>>>lora_decode
Hit up！！20406   24502
Hit up！！24502   28599
Hit up！！28599   32695
Hit up！！32695   36791
Hit up！！36791   40887
Hit up！！40887   44983
Hit up！！44983   49079
Hit up！！49079   531760   1.1271e+024.0960e+03   7.8453e+018.1920e+03  -1.6873e+011.2288e+04  -1.6863e+021.6385e+04  -4.3814e+012.0481e+04   6.2779e+012.4577e+04   1.0725e+022.8673e+04   8.4182e+013.2769e+04   4.2651e+003.6865e+04  -1.5306e+024.0962e+04  -6.7274e+01
xcorr pre-upchirp amp fft payback=3.97002e-05
Hit down！！65187   69283
Fine Time Freq begin...Finished.Columns 1 through 23:670    826    590    146   1818    326   1526   1170    966    950    858   2033    657    966   1101   1854    650   1563   1361   1820   1574    331   1120Columns 24 through 46:1761    133   1753    922   1027   1564    622    984   1461   1831   1472   1690    234     90    110   1004    325   1917   1740    908    142    261   1955Columns 47 through 69:140   1916    835   1426   1628    216   1417   1575   1021   1650    217     35    218    559    689   1294   1114   1590    304   1753    148   1437   1625Columns 70 through 75:543    369   1994    609    609   1015

解调后的符号，通过 LoRaPHY提供的解析方法，可以还原出实际的数据:

>>> lora_decode
Syn code start.start offset = 15, fix = 0
01ECAE9230011111111111111112222222222222222333333333333333344448F16
078357676000123456789ABCDEF0123456789ABCDEF0123456789ABCDEF0123CFF4

值得注意的是，上述解调后，从2.25符号往后，其实有15个符号（30样点）的偏差。这个偏差是固定的，可以直接根据头部校验搜索并找到。这个偏差的产生机制，并不是很清楚。

7 拉远距离（加白噪声）测试

我们添加一个比波形本身强度高10倍的噪声，相位已经没法看了，但通过相关，发现依旧能正确处理：
在10倍的噪声下，首部的相位也有抖动，但大的周期还在。偶尔估计会出现偏差，但依旧能够较大概率解析出正确的数据（会出错）。

8 代码获取和运行

参考代码库：a4_lora/octave

从本次摸索，我发现读别人的文章与自己摸索实现一下，认识程度是显著不同的。我想LoRa本身可能也采取了类似的时频差补偿算法，否则它那60多块钱的成本，以及巨大的参数公差，根本就收不好。

自己的Octave实现与网上的例子相比，既不要求安装matlab R2019以上这种20GB的工程数学软件，也不和GNU-Radio这种多层封装的实现发生关系。

要运行代码，请直接安装开源的计算工具 GNU-Octave （和Matlab很像，但没有版权问题，安装包在几百兆左右）。Octave也可以直接使用MSYS2的pacman -S来安装。本次的代码不依赖任何工具箱，包括communicaitons。