设计被动LiDAR检测器设备 - 固件篇

欢迎回来！在上篇文章中，我描述了检查iPhone Pro TrueDepth LiDAR的过程，并构思了检测它的硬件设计方案。本文将介绍该设备的固件、其中采用的概念以及在微型SAMD21微控制器上实现这些概念的方法。

让我们先回顾上一篇文章中的一个关键点：“LiDAR是一种闪烁的光，我们如何测量闪烁的光”

1. 将光信号视为散布在区域上的光束
2. 感应这些光束并将其转换为模拟信号
3. 将信号从模拟转换为数字
4. 测量它们

事实证明，不仅仅是检测闪烁的光，而是检测多个离散的闪烁光束并将其与其他噪声区分开来。上一篇文章中描述的硬件提供了来自4个离散光电二极管的超快数字信号。我们可以使用这些信号进行多种不同的测量，以区分iPhone Pro TrueDepth LiDAR与其他信号源或噪声的特性。我为该设备设计的固件采用了以下测量方法：

频率（Goertzel算法）：

• 目标60Hz，可选谐波30Hz和120Hz
• 注意：包含谐波会增加误报的敏感性

脉冲重复频率（PRF）：

• 当LiDAR信号存在时，我们期望它看起来稳定，边缘均匀，而像红外遥控器这样的信号更多地是突发性的

突发相似性（幅度四分位距）：

• 同样，我们期望信号在出现时看起来稳定，因此对最近的测量结果进行分析，以判断它们代表的是稳定信号还是突发信号
• 比较汉恩窗口平滑、去趋势后信号能量的顶部和底部四分位数
• 我对这个理解得不是很好，但通过网页搜索和StackExchange，我摸索出了方法

空间巧合性（离散光束与广域泛光对比）：

• 保持每个通道的指数移动平均峰值记录，并检查以确定是部分而非全部传感器同时触发

TSOP红外解码否决（有时暴力解法就是优雅的解决方案）：

• 当TSOP在测量窗口内成功解码遥控代码时，认为信号源极不可能是LiDAR

周期性闪烁（可检测的占空比）：

• iPhone Pro TrueDepth LiDAR似乎以一种可检测的方式闪烁开关，但这不太可靠，因为在测量窗口期间，手持信号源的人只要简单移动，信号就可能进入或退出检测范围

为了完成所有这些任务，固件利用了几个简单的概念来抽象和组织来自设计硬件的测量结果，并以有用的方式对其进行分析。

评分系统： 为了有意义地分析测量结果并以编程方式决定给定信号是否是所寻找的LiDAR设备，固件采用了加权评分系统。必须在测量窗口内达到超过给定阈值的分数才能触发检测。

• 对所需测量的积极检测，例如所需的主导频率、测量窗口内部分而非全部传感器触发等，会给予一定的正分。
• TSOP对红外遥控码的积极解码将给予显著的负分。

测量分桶（采样）： 测量事件被存储到固定时间（4毫秒）的样本或“桶”中。

环形缓冲区（快速读，快速写）： 测量样本或“桶”使用环形缓冲区进行写入和读取。这很快！它需要快，这样才不会浪费硬件的能力。

中断驱动信号采集（快速边沿捕获）： SAMD21可以在中断计时器上实现微秒级精度。这很快，而这也需要快！

锁存滞后（防抖动/防“Flappy Bird”效应）： 在早期模型中，我注意到检测要么会突然开启然后永不关闭，要么几乎不可能保持触发状态，只能通过串口输出观察到，而不能通过LED信号看到。这不好玩，它需要既可爱又能亮起来。为了解决这个问题，我使用计数器来测量：

• 在给定时间段内实现了多少次积极检测
• 在同一时间段内没有触发多少次

基于这些计数，当检测达到足够长的时间时，“检测”指示器被触发并锁存；当检测长时间未达到时，指示器解锁存。这防止了一些误报，并有助于平滑成功检测的指示。

这确实涉及很多内容，而且在开始这个项目之前，我对其中很多内容都不熟悉，所以请允许我更详细地介绍其中一些概念。

频率分析（FRE）

在我们讨论频率相关的事情之前，我们需要谈谈数学（抱歉，我也不喜欢这个）。为了在软件中获取我们的频率测量值，我们需要使用Goertzel算法。但在我们有效做到这一点之前，首先对测量信号进行去趋势处理和汉恩平滑处理。

去趋势 / 直流偏移消除 这到底是什么意思？我不得不查一下，这让我感到尴尬，因为我在数字音频工作站（DAW）中录制音频时也做这个，但不知为何在查资料之前我没想到这里的联系，而且竟然是从Audacity网站上了解到的。

“直流偏移是信号相对于零点的偏移。在Audacity波形图中，这意味着默认视图中的波形似乎不以0.0水平线为中心，如本图像中的上轨道所示：”

这是为了从信号中去除直流分量，以防止0Hz基线干扰噪声中60Hz信号的测量。我对此理解不深，这需要大量递归的网页搜索来学习如何测量和分析信号。看起来它似乎有效，但我的理解相当肤浅。

为此，执行均值去除，这听起来是积极的！均值去除（去趋势）以最简单的方式完成，通过计算窗口中所有样本的平均值，然后从每个样本中减去该平均值，如下所示：

1
2
3
4
5
6
7
8

float mu = 0.0f;
for (int i=0; i<N; i++) {
 mu += prep[i];
}
mu /= N;
for (int i=0; i<N; i++) {
 prep[i] -= mu;
}

去趋势处理后，对测量结果应用汉恩平滑处理。

汉恩平滑 同样，我对此理解不深，但基于所有的网页搜索，汉恩平滑是一种平滑信号波形的方法，使其变得锥形平滑，而不是古怪且参差不齐。设备在分析当前窗口、均值去除之后、执行Goertzel之前进行汉恩平滑。这样，观察窗口的边缘具有平滑的锥度，没有尖峰，否则尖峰会泄漏到相邻频率并阻碍测量。没有汉恩平滑，检测会非常跳跃，很难成功检测。有了汉恩平滑，主导频率区间就清晰得多。

事情是这样的，我不是科学家，也不是软件开发人员，我只是个黑客。我不得不通过网页搜索来摸索这一切。这显然就是你在C语言中做这件事的方法。

1
2
3
4

for (int i=0; i<N; i++) {
  float w = 0.5f * (1.0f - cosf(2.0f * PI * i / (N - 1)));
  prep[i] *= w;
}

就这么简单！？显然如此！在执行了当前分析窗口的均值去除和汉恩平滑之后，设备使用Goertzel算法来定位哪些频率最强！

Goertzel算法 Goertzel滤波器是一种在软件中从信噪比中滤除信号的方法。在这种情况下，设备有四个光电二极管接收大量光子并产生大量信号，微控制器则在每次发生时进行记录。这是很多噪声。设备在信号中寻找特定频率，特别是60Hz及其可选谐波（30Hz, 120Hz等）。这将告诉我们观察到的信号中60Hz或其谐波的强度，这在C语言中相当容易完成，正如StackExchange上的这个例子所示。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

    float   scalingFactor = numSamples / 2.0;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q0 = coeff * q1 - q2 + data[i];
        q2 = q1;
        q1 = q0;
    }

    // calculate the real and imaginary results
    // scaling appropriately
    real = (q1 - q2 * cosine) / scalingFactor;
    imag = (q2 * sine) / scalingFactor;

    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}

在LiDAR检测器固件中，这种滤波用于检查特定目标Hz的功率 - Goertzel功率！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

static float goertzel_power(const float* x, int N, float targetHz, float fs){
  float kf = (N * targetHz) / fs; // bin index fraction
  int k = (int)(kf + 0.5f); // nearest
  if (k <= 0 || k >= N-1) {
    return 0.0f;
  }
  float w = 2.0f * PI * k / N;
  float c = 2.0f * cosf(w);
  float s0 = 0;
  float s1 = 0;
  float s2 = 0;
  for (int i = 0; i < N; i++) {
    s0 = x[i] + c*s1 - s2; s2=s1;
    s1=s0;
  }
  float P = s1*s1 + s2*s2 - c*s1*s2;  // power!
  return ( P > 0 ) ? P:0;
}

设备检查60Hz区间（以及可选的如前面提到的谐波等其他区间）是否强且稳定（主导）。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

// find in 30, 60, 24, 50, or 120Hz
float topP = 0
float runnerP = 0;
int topIdx = -1;
for (int i=0;i<FRE_BINS_N;i++) {
  float p = goertzel_power(prep, N, FRE_BINS_HZ[i], MICRO_FS); // Goertzel Power!
  if (p > topP) {
    runnerP=topP;
    topP=p;
    topIdx=i;
  } else if (p > runnerP) {
    runnerP=p;
  }
}
float dom = (runnerP>0)? topP/runnerP : (topP>0?1e9f:0.0f);
if (topIdx>=0 && dom>=FRE_DOM_RATIO) {
  freScore = SCORE_FRE;
}

通过调整滑动窗口内检查各区间功率的“桶”的时间间隔，我进行了调谐，以从噪声中筛选出持续高于阈值的60Hz区间。在测试中，这似乎是在200毫秒滑动窗口上以4毫秒间隔进行的。

突发相似性（BUR）

红外遥控器会显得“突发性”（或者根据某个LLM的说法是“spurty”，我选择在这里包含这个词），而像LiDAR这样的信号在占空比内可能看起来更稳定。我想尽可能测量所有能将LiDAR信号与其他一切区分开来的特征，所以我也递归地搜索了如何测量这个以及如何在C语言中实现。

使用与频率分析相同的去趋势和汉恩平滑后的分析窗口，我们制作一个“直方图”（我原以为这只是在TS-100中才有），然后进行幅度四分位距分析。根据我对这个的理解，我们想做的是：

1. 确定第25和第75个百分位数是多少
2. 将大于第75个百分位数的样本平均在一起
3. 将小于第25个百分位数的样本平均在一起
4. 用前者除以后者。

在C语言中，这并不复杂，如下所示。

1

burstRatio = average(samples_greater_than_Q3) / average(samples_less_than_Q1)

哇，那很容易！带着对如何做这件事的“理解”，并且在进行了更多关于如何在C语言中做这件事的递归网页搜索之后，我搞出了一个混乱的代码！设备使用数学计算出一个突发比率，然后将其与一个值（BURST_RATIO_NEED）进行比较，我仍然不确定什么值是恰到好处的，然后根据它评分。看看这个穴居人代码 / StackExchange / LLM解决方案的大杂烩：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

// prep[] after detrend + Hann-window of current window

// min/max for histogram
float mn = prep[0];
float mx = prep[0];
for (int i=1; i<N; i++) {
  if (prep[i] < mn) {
    mn = prep[i];
  }
  if (prep[i] > mx) {
    mx = prep[i];
  }
}

// build histogram
const int HB = 16;
int hist[HB] = {0};
if (mx == mn) {
  return 0;
} else {
  for (int i=0; i<N; i++) {
    int b = (int)((prep[i] - mn) * (HB - 1) / (mx - mn));
    if (b < 0) { 
      b = 0;
    } else if (b >= HB) {
      b = HB - 1;
    }
    hist[b]++;
  }
}

// find q1 and q3 in historgram
int q1T = N / 4;
int q3T = (3 * N) / 4;
int acc = 0;
int q1Idx = 0;
int q3Idx = HB - 1;

for (int b = 0; b < HB; b++) {
  acc += hist[b];
  if (acc >= q1T) { q1Idx = b; break; }
}
acc = 0;
for (int b = 0; b < HB; b++) {
  acc += hist[b];
  if (acc >= q3T) { q3Idx = b; break; }
}

// convert bin index to threshold
float q1 = (mx > mn) ? (mn + (mx - mn) * (q1Idx / (float)(HB - 1))) : mn;
float q3 = (mx > mn) ? (mn + (mx - mn) * (q3Idx / (float)(HB - 1))) : mx;

// averages of quartiles
float sumLo = 0;
float sumHi = 0;
int cntLo = 0;
int cntHi = 0;
for (int i = 0; i < N; i++) {
  float a = fabsf(prep[i]);
  if (prep[i] <= q1) {
    sumLo += a;
    cntLo++;
  }
  if (prep[i] >= q3) {
    sumHi += a;
    cntHi++;
  }
}

// compare observed ratio to limit and score
float burstRatio = 0.0f;
if (cntLo > 0 && cntHi > 0) {
  float meanLo = sumLo / cntLo;
  float meanHi = sumHi / cntHi;
  burstRatio = (meanLo < 1e-6f) ? ((meanHi > 0) ? 1e9f : 0.0f) : (meanHi / meanLo);
}
int burstScore = (burstRatio >= BURST_RATIO_NEED) ? SCORE_BURST : 0;

脉冲重复频率（PRF）

“脉冲重复频率是重复信号在特定时间单位内的脉冲数。”

我原以为这会像占空比或旋转LiDAR可能产生的效果一样。令我懊恼的是，PRF显然不同。它发生在微小时的时间尺度上，而来自旋转LiDAR或占空比的闪烁接近秒级时间尺度，对于可怜的人眼来说是明显的。这将被观察为信号随时间变化的某种方差或缺乏方差，而不是用我们凡人的眼睛感知到。

设备收集给定分析窗口内所有边沿中断触发的时间戳。这允许它枚举该窗口内所有的脉冲间隔。为了这个目的，设备还将所有四个通道合并为一个信号，因为光的结构对于这部分并不重要。然后它又进行了一些数学计算（抱歉，这里面有太多数学，我不得不查所有这些，因为我上的是商学院），将间隔转换为频率并存储在缓冲区中。

然后，甚至进行更多的数学计算！它根据这些缓冲的间隔频率值计算PRF和变异系数，并将其与代表所需变异紧密度或PRF阈值的预定义值进行比较。这个值越稳定，信号是LiDAR的可能性越高。

我查了很多资料来弄清楚如何在C语言中做到这一点，再说一次，我是个黑客，不是开发人员，所以不能保证这个代码写得高效甚至写得好：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

// edge interrupt timestamps for window
float freqs[MAX_PULSES];
float int nf=0;
for (int i=1; i<nEdges; i++){
    float dt = (edges[i].t - edges[i-1].t) * 1e-6; // seconds
    if (dt > 1e-6 && dt < 5e-3) { // 200khz to 200hz
        freqs[nf++] = 1.0f / dt; // hz
    }
}

// compute mean PRF and CV
float sum = 0;
float sum2 = 0;
for (int i = 0; i < nf; i++) {
  sum += freqs[i];
  sum2 += freqs[i] * freqs[i];
}
float mean = sum / nf;
float var = (sum2 / nf) - mean * mean;
float sd = sqrtf(fmaxf(var, 0));
float cv = ( mean > 1e-3 ) ? (sd / mean) : 1.0f;

// award score 
int prfScore=0;
if (nf >= MIN_PULSES && mean >= 5e3 && mean <= 40e3 && cv <= 0.10) {
    prfScore = +3;
}

空间巧合性（SPAT）

正常光源对于像光电二极管这样的传感器来说看起来像单一光源。当你将像电脑显示器或电视遥控器指向红外传感器，或者在本例中指向红外传感器阵列时，很可能所有传感器都会同时看到这个相同的光源。正如我之前描述的，iPhone Pro上的LiDAR光源投射的是点阵，更像是聚焦的光束，而不是广泛散布的光源。

因此，从适度距离直视LiDAR光源，更像是看向投影仪，正好处于进入和离开其发出的刺眼光线的位置。向左移动可能被强光覆盖，向右移动可能直接看向发射器而看不到信号。这种特性使得设备能够区分可能的LiDAR信号源与其他杂散信号（噪声）！它将单一传感器在时间帧内触发边沿视为可能性低，或者所有传感器在同一时间帧内触发边沿也视为可能性低，而2-3个传感器在时间帧内触发边沿则视为高可能性。

设备记录每个传感器的边沿时间戳，并针对给定的分析窗口比较是否有任何传感器在同一时间帧内相互触发。我们检查每个时间帧的巧合计数，当有2个或更多但并非全部传感器触发时，给予正分。设备可选地，当零个、一个或全部传感器在给定时间帧内触发时，给予负分。

TSOP红外解码（否决）

即使设备在寻找60Hz信号，并且我们在硬件层面过滤940nm红外线，我们仍然会遇到一些常见的误报情况。电视通常是60Hz或120Hz，美国电灯泡的交流电频率显然是60Hz左右，红外遥控信号似乎也通常在60Hz左右。至少，红外遥控码信号可以很容易地被否决，因为它们可以非常容易且高度自信地被检测到。对于这部分，我基本上没写什么代码，只是借用了一个为TSOP1738解码红外遥控信号的Arduino库：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#include <IRremote.hpp> // yolo
static uint16_t tsop_remote_frames = 0;
static uint32_t tsop_window_start_ms = 0;

static int scoreTSOP() {
  uint32_t now = millis();
  if (tsop_window_start_ms == 0) {
    tsop_window_start_ms = now;
  }
  if (now - tsop_window_start_ms >= TSOP_WINDOW_MS) {
    tsop_remote_frames = 0;
    tsop_window_start_ms = now;
  }
  return (tsop_remote_frames > 0) ? -3 : 0;
}
// In loop()
while (IrReceiver.decode()) {
  if (IrReceiver.decodedIRData.protocol != UNKNOWN) {
    tsop_remote_frames++;
  }
  IrReceiver.resume();
}
// at score time
int TSOP_Score = scoreTSOP();

周期性闪烁（可检测的占空比）

我仍在努力理解如何最好地实现通过网页搜索建议的两种检测信号这方面特征的方法：用于检测强周期信号的二进制自相关，以及用于检测一致占空比测量的游程分析。

设备在用于测量PRF的每个“桶”能量导出的组合二进制能量流上执行这些分析。主要问题是这种分析发生在更长的时间框架上。目标设备需要保持静止并被观察比空间巧合或频率分析所需长得多的时间，因此作为测量的价值较低。由于分析价值较低，花费处理能力来执行它可能并不值得。

环形缓冲区、窗口、桶

设备使用环形缓冲区非常快速地存储和读取原始、已处理和已分析的测量结果。来自特定时间分段（桶）的小测量值被收集用于给定的测量窗口。这种数据组织允许设备在保持已知过去测量的一定状态的同时，维持测量和分析的滑动窗口。

中断驱动

设备需要能够进行信号处理并且做得快，以利用传感器的速度。SAMD21有一个48MHz的CPU来处理其中断引脚。这非常快。设备依靠中断来测量来自传感器通道电路的信号，以便在软件端尽可能快。

滞后 / 锁存

设备在检测触发（分数超过阈值）和解除检测（分数未超过阈值）方面做了一些滞后处理（我为这个项目新学的这个词）。

如前所述，设备维护两个滑动窗口计数，记录给定时间段内有多少次测量结果导致分数高于阈值，以及有多少次没有。当高于阈值的正触发次数超过另一个阈值时，设备将指示。一旦在给定窗口内有给定次数的测量结果未能导致分数高于阈值，设备将停止指示。

演示

以下是硬件设计中运算放大器和施密特触发器版本的概念验证设备的几个视频！

在每个视频中，原型设备可以在画面中间看到，LiDAR光源（iPhone 15 Pro）在前景中，位于画面下方的三分之一处。LiDAR光源投射的点阵网格以闪烁的点状形式在可见表面上可见。可以看到一个蓝色LED快速闪烁，指示串行传输。来自设备的串行输出以滚动文本的形式覆盖在画面上。

该文本包括来自原型设备的实时测量和分数信息。在视频中，LiDAR光源被启动（iPhone相机应用被打开），不久之后另一个蓝色LED亮起，指示设备已检测到TrueDepth LiDAR信号。在此刻，屏幕上显示的文本表明设备在频率分析、空间巧合和突发性方面授予了足够高的分数，并且从TSOP获得了足够低的分数，从而触发了检测指示。

视频不可用或格式不受支持。请尝试其他浏览器。

运算放大器+比较器概念验证LiDAR检测器

视频不可用或格式不受支持。请尝试其他浏览器。

施密特触发器概念验证LiDAR检测器

还有很多可以改进的地方，无论是在可以使用哪些组件方面，还是在固件调优方面。总体而言，我对结果相当满意，尽管它们还很粗糙。这很有趣，我学到了很多对其它项目已经派上用场的东西。

感谢阅读！

向Samy Kamkar致敬，他激励我尝试这个项目。谢谢你，Samy！ JFIO 帮派 r%！学习很有趣！知识就是力量！黑掉这个星球！

附录 I

• https://www.st.com/resource/en/design_tip/dt0089-the-goertzel-algorithm-to-compute-individual-terms-of-the-discrete-fourier-transform-dft-stmicroelectronics.pdf
• https://en.wikipedia.org/wiki/Goertzel_algorithm
• https://en.wikipedia.org/wiki/Discrete_Fourier_transform
• https://www.reddit.com/r/AskElectronics/comments/g5vrio/how_to_calculate_rc_values_to_remove_dc_offset/
• https://en.wikipedia.org/wiki/Hann_function
• https://en.wikipedia.org/wiki/Interquartile_range

附录 II

现有技术 这些是我能找到的似乎与检测LiDAR系统有关的专利。但这些似乎都没有检测TrueDepth LiDAR。

• US 5,666,120 A / US 5,990,821 A – “检测器设备。” 汽车雷达/激光探测器，带有激光（LIDAR）光电二极管前端和透镜阵列
• US 6,670,905 B1 – “带有位置和导航信息的雷达预警接收器。” 雷达探测器明确包含用于“检测警察LIDAR信号”的激光探测器，并集成GPS。
• US 2020/0041622 A1 – “检测和跟踪LiDAR串扰。” 使用无源通道和时间相关性来感知和分类外部LiDAR源（例如，其他车辆）的方法——功能上
• US 2020/0158825 A1 – “LiDAR系统的干扰检测和缓解。” 检测外部LiDAR干扰并缓解它
• US 11,899,110 B2（及相关专利） – “用于检测LiDAR成像传感器中的串扰干扰的技术。” 通过分析时间模式/串扰来检测外部LiDAR
• US 2020/0073531 A1 – “用于深度感测的结构光检测。” 专注于用相机检测结构光图案（例如，点阵投影仪）