用于低成本接收机的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 .实验学习结论摘要
- 1. 波形样本来源
- 2. 粗略显示LoRa物理层突发的结构
- 2.1 chirp的数学描述
- 2.2 chirp的归一化离散形式
- 3. 在频差条件下检测头部
- 3.1 钟差频差对头部检测几乎没有影响
- 3.2 检测结果
- 4. 基于相关峰起伏周期的钟差估计算法
- 4.1 低成本接收机的困境
- 4.2 相关峰起伏周期的成因
- 4.3 计算结果
- 5 基于 up-down chirp 边缘峰值交汇的时频差精确估计算法
- 5.1 笨办法观察规律
- 5.2 基于边缘最大值直线交点快速交汇的时频差估计
- 6 卡住位置补偿解调
- 7 拉远距离（加白噪声）测试
- 8 代码获取和运行
0 .实验学习结论摘要
虽然是站在LoRa-PHY、GNU-Radio巨人的肩膀上，但是在近1个月的不断踩坑、尝试中，还是非常的不顺利。最终得到正确结果时，感觉就像攀过了一座陡峭的小山。
- 我们发现对于Lora这种使用Chirp的起始频率携带信息的低成本硬件，其固有的时钟差、频率差对结果影响极大。尤其是使用另一个不靠谱的硬件来接收（比如我的山寨 B205mini），两者的钟差是不靠谱+不靠谱=离谱。能在Low SNR下解析出正确的数据，需要对波形的深刻理解和认识，每一个步骤都值得思考。
- 对我学习的清华大学项目 LoRa-PHY里，用于同步时钟误差的方法，是假设主要的误差都来自LoRa设备，用频差去折算钟差（时差）。在本实验里，我们的接收设备也很不靠谱，使得这个举措变得不稳定。我们通过观察相关峰值的最大值的周期变化，发现了隐含在头部14组up-chirp里的低频周期，提出了基于相关峰起伏周期的钟差估计算法，并使用它对钟差进行补偿，效果很好。
- 对收发时钟、频率稳定度都很低的情形，常见的二维时频相关算法耗时太大。提出基于 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 秒，他们满足:
B T = 2 S F B T = 2^{SF} BT=2SF
在T秒内，基带IQ向量旋转的频率从 -B/2 变化到 B/2。
由于频率
f ( t ) = 2 π ( f 0 + B T t ) f(t)=2\pi(f_0+ \frac{B}{T}t ) f(t)=2π(f0+TBt)
而向量的相位是频率的积分，计算：
p ( t ) = ∫ f ( t ) d t = 2 π ( f 0 + B 2 T t ) t , t ⊂ [ 0 , T ) p(t)=\int f(t)dt=2\pi(f_0+ \frac{B}{2T}t )t, t\subset [0,T) p(t)=∫f(t)dt=2π(f0+2TBt)t,t⊂[0,T)
此时，样点的复平面向量可以表示为：
s ( t ) = e 2 π j ( f 0 + B 2 T t ) t s(t)=e^{2\pi j(f_0+ \frac{B}{2T}t )t} s(t)=e2πj(f0+2TBt)t
这个t的取值是 0 ~ T秒。
2.2 chirp的归一化离散形式
使用上述模拟波形的样式来计算chirp, 不容易出错。但如果换算为离散形式，则更方便写程序。假设我们采用高于LoRa带宽B的K倍速率采样，即 Sr = BK，样点为 n，则参考本文开头的第一个式子里SF、B、T的关系，一个chirp的样点数为：
S Y M = N ( S F , K ) = 2 S F B S r = 2 S F ⋅ K SYM=N(SF,K)=\frac {2^{SF}}{B} Sr=2^{SF}\cdot K SYM=N(SF,K)=B2SFSr=2SF⋅K
比如SF=11, K=2时，N就是4096点。根据离散、连续的关系，
t = n S r = n B K t=\frac{n}{Sr}=\frac{n}{BK} t=Srn=BKn
p ( t ) = 2 π ( f 0 + B 2 T n B K ) n B K = 2 π ( f 0 + n 2 T ⋅ K ) n B K 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} p(t)=2π(f0+2TBBKn)BKn=2π(f0+2T⋅Kn)BKn
由于 B T = 2 S F BT=2^{SF} BT=2SF, 有：
p ( t ) = 2 π ( f 0 B + n 2 ⋅ 2 S F ⋅ K ) n K p(t)=2\pi (\frac{f_0}{B}+ \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K} p(t)=2π(Bf0+2⋅2SF⋅Kn)Kn
剩下一个 f0/B，也可以兑换为比例，其实f0取-B/2时，他就是 -1/2。
经过上面的推导，我们把 chirp 表示为和具体采样率、频率无关的值，只和SF,K有关。
p u p − c h i r p ( t ) = 2 π ( − 1 2 + n 2 ⋅ 2 S F ⋅ K ) n K p_{up-chirp}(t)=2\pi (-\frac{1}{2}+ \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K} pup−chirp(t)=2π(−21+2⋅2SF⋅Kn)Kn
S u p − c h i r p ( n ) = e 2 π j ( − 1 2 + n 2 ⋅ 2 S F ⋅ K ) n K S_{up-chirp}(n)=e^{2\pi j (-\frac{1}{2}+ \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K}} Sup−chirp(n)=e2πj(−21+2⋅2SF⋅Kn)Kn
我们可以使用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查看有无单音来检测，那样速度更快。
检测这个步骤，使用共轭，能够和原始波形直接合成为很强的峰。
S u p − c h i r p ‾ ( n ) = e − 2 π j ( − 1 2 + n 2 ⋅ 2 S F ⋅ K ) n K \overline{S_{up-chirp}}(n)=e^{-2\pi j (-\frac{1}{2}+ \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K}} Sup−chirp(n)=e−2πj(−21+2⋅2SF⋅Kn)Kn
这个波形的瞬时频率为：
f u p − c h i r p ‾ ( n ) = − 2 π ( − 1 2 + n 2 S F ⋅ K ) 1 K \overline{f_{up-chirp}}(n)=-{2\pi (-\frac{1}{2}+ \frac{n}{2^{SF}\cdot K})\frac{1}{K}} fup−chirp(n)=−2π(−21+2SF⋅Kn)K1
我们在第一阶段，主要的工作就是用上面这个相关波形去和原始波形做相关，两两相乘而后求和，检测峰值。
3.1 钟差频差对头部检测几乎没有影响
即使在有频偏的情况下，这种检测依然是有效的。只要频率偏移不是很离谱，通过样点的移动，总能用时间延迟的损失来抵消频偏：
设存在一个时间偏移 n ′ = δ + n n'=\delta+n n′=δ+n, 以及一个固有的频率偏移 Δ \Delta Δ：
f u p − c h i r p ‾ ( n ′ ) = f u p − c h i r p ‾ ( δ + n ) = − 2 π j ( Δ − 1 2 + n + δ 2 ⋅ 2 S F ⋅ K ) 1 K \overline{f_{up-chirp}}(n')=\overline{f_{up-chirp}}(\delta+n)=-{2\pi j (\Delta-\frac{1}{2}+ \frac{n+\delta}{2\cdot2^{SF}\cdot K})\frac{1}{K}} fup−chirp(n′)=fup−chirp(δ+n)=−2πj(Δ−21+2⋅2SF⋅Kn+δ)K1
当这个 δ \delta δ满足下式，使得频偏 Δ \Delta Δ其恰好能被时间推移抵消：
Δ + δ 2 ⋅ 2 S F ⋅ K = 0 \Delta+\frac{\delta}{2\cdot2^{SF}\cdot K}=0 Δ+2⋅2SF⋅Kδ=0
则结果变为
f u p − c h i r p ‾ ( n ′ ) = − 2 π ( − 1 2 + n ′ − δ 2 ⋅ 2 S F ⋅ K ) 1 K \overline{f_{up-chirp}}(n')=-{2\pi (-\frac{1}{2}+ \frac{n'-\delta}{2\cdot2^{SF}\cdot K})\frac{1}{K}} fup−chirp(n′)=−2π(−21+2⋅2SF⋅Kn′−δ)K1
相关峰与无频偏时相比，只是提前或者推迟了 δ \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 π j Δ n e^{2\pi j \Delta n} e2πjΔn,而后，再平移 δ \delta δ样点，计算共轭相乘的和。
在窗口内，遍历所有的 Δ , δ \Delta,\delta Δ,δ,
D ( m , n ) = e 2 π j ( − 1 2 + Δ + m + δ 2 ⋅ 2 S F ⋅ K ) m + δ K ⋅ e − 2 π j ( − 1 2 + n 2 ⋅ 2 S F ⋅ K ) n K D(m,n)=e^{2\pi j (-\frac{1}{2}+\Delta+ \frac{m+\delta}{2\cdot2^{SF}\cdot K})\frac{m+\delta}{K}}\cdot e^{-2\pi j (-\frac{1}{2}+ \frac{n}{2\cdot2^{SF}\cdot K})\frac{n}{K}} D(m,n)=e2πj(−21+Δ+2⋅2SF⋅Km+δ)Km+δ⋅e−2πj(−21+2⋅2SF⋅Kn)Kn
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。