第2章 嵌入式软件编程模式和优化
本章介绍嵌入式软件编程模式和通用软件优化方案。嵌入式软件编程模式关注的是底层支撑软件的架构，包括内存和CPU运行时间的分配方案，它服务于上层机器学习软件。嵌入式软件编程模式和传统的计算机编程模式在概念上有所不同，一般计算机软件编程模式侧重软件开发过程中的代码复用和架构标准化技术，以模块化和层次化的形式构建软件对象，而这一章所讨论的嵌入式编程模式是为了满足多任务的实时响应，在有限内存和有限运算能力的嵌入式CPU上高效地完成运算任务。本章讨论的内容除了编程模式外，还包括通用的嵌入式软件优化方案，它们和嵌入式软件编程模式相结合，保障了机器学习算法在嵌入式平台上高效运行。这一章的内容除了可用于提高机器学习应用软件效率外，也能够应用于其他具有类似特性的嵌入式软件。
2.1 嵌入式软件编程模式
嵌入式系统和应用紧密结合，不同应用有不同的软件运行要求。下面将讨论几种常见的嵌入式系统编程模式，每种模式有对应的应用场合，可以单独使用或者混合使用。这里讨论的编程模式主要针对没有操作系统的嵌入式软件运行环境，在这种情况下，CPU的全部算力可以分配到和应用相关的计算，不需要额外执行IO资源状态、内存清理、调度等软件操作系统的管理任务，因此运行效率和内存使用效率会更高，但付出的代价是需要手动管理任务并发、IO状态检查、资源共享等，对开发者有更高的要求。
下面首先给出常见的嵌入式系统运行模式示意图，如图2-1所示。
图2-1 常见的嵌入式系统运行模式
图中系统的功能是从传感器获得外部数据，进行分析运算后输出控制命令，运行期间还需要接收用户输入。这一架构模型反映了大多数嵌入式系统的运行模式。比如，在语音识别应用中，嵌入式系统周期性地从ADC读入语音的波形数据，经过处理后识别用户语音指令，并根据识别结果输出控制命令。对于运行机器人视频避障软件的嵌入式系统，软件需要通过摄像头周期性地获得视频数据，分析视频内容，识别障碍和目的地并输出机器人移动控制命令。
上面描述的嵌入式系统运行模式涉及了多个需要并行执行的任务，每个任务有不同的运行周期性和处理优先级要求。表2-1给出了一个例子，其中嵌入式系统中需要同时运行4个任务，分别是传感器采样数据读取、数据计算及控制输出、用户输入命令检查、状态显示更新。它们的执行周期如表2-1所示。
表2-1 任务类型和执行周期示例
这几个任务中，读入传感器数据的时间优先级最高，根据采样率所规定的时间间隔运行，不能错过数据；计算及控制输出同样有实时性要求，但它的执行频率低于数据采样；用户命令的检查和状态显示更新的优先级最低，只要满足人们对信息的感知速度即可。
上面几个任务的时间要求和运行频率要求也各不相同。设计时，需要根据优先级制定软件执行方案。当使用实时操作系统时，可以使用比如“单调速率调度算法”（Rate Monotonic Scheduling，RMS）等方案实现以上多任务并发执行，如图2-2所示。我们后面讨论的几种编程模式能够不依赖于操作系统实现，允许开发者根据应用特点实现更灵活的定制化多任务执行模式。
图2-2 多个需要周期重复的任务交替执行的示意图
2.1.1 基于周期调用的运行模式
周期运行模式的一种实现方式如代码清单2-1所示，该代码实现了表2-1给出的嵌入式系统的几个任务周期执行。其中while循环内部的几个任务按执行顺序排列，并且每轮while循环根据周期要求对各个任务执行1次或者多次。
代码清单2-1 基于周期调用的运行模式伪代码
while(1)
{
传感器采样数据读取();
用户输入命令检查();
传感器采样数据读取();
数据计算及控制输出();
传感器采样数据读取();
状态显示更新();
传感器采样数据读取();
数据计算及控制输出();
}
while循环内的代码对应所需要的任务的一个运行周期，如图2-2所示。但在实际情况中，由于每个模块运行时间存在差异，因此往往得到疏密不均匀的“周期”运行效果，如图2-3所示。这使得系统对传感器输入数据的“采样率”不再固定，并且对用户的输入响应间隔也时快时慢，系统运行缺少“确定性”。改进措施如下：
1）将运行时间过长的任务进行状态分割。
2）将采样固定在时间格点上运行。
图2-3 各个任务运行时长不一致造成运行周期疏密不均匀
下面详细介绍这两种改进方案。
·运行时间长的任务拆分
任务拆分是将一个长时间任务分成几个阶段实现，每次循环执行时仅完成其中一个阶段的任务。例如，假设“状态显示更新”任务由于具有LCD特性而运行很慢，我们将其分成几部分，比如假设原始的操作可以分成三部分，如代码清单2-2所示。
代码清单2-2 运行缓慢的显示更新任务的伪代码
void状态显示更新()
{
更新测量数据曲线图();
更新数据计算结果();
更新用户输入内容();
}
代码清单2-2所显示的三部分操作可以分成三个阶段执行，如代码清单2-3所示。
代码清单2-3 拆分成三个阶段执行的显示更新任务的伪代码
void状态显示更新()
{
static int state=STAGE0;
switch (state):
{
case STAGE0:
更新测量数据曲线图();
state=STAGE1;
break;
case STAGE1:
更新数据计算结果();
state=STAGE2;
break;
case STAGE2:
更新用户输入内容();
state=STAGE0;
break;
}
}
该函数内部有static类型状态变量state，它的取值在函数一次调用完成到下一次调用期间是保持不变的。我们第一次调用该函数时，由于state=STAGE0，因此执行“更新测量数据曲线图”；第二次调用时，由于state=STAGE1，因此执行“更新数据计算结果”；第三次调用时，由于state=STAGE2，因此执行“更新用户输入内容”；三次调用后，state回到STAGE0，下次调用再次执行“更新测量数据曲线图”，并重复上面的过程，如图2-4所示。
图2-4 将运行时间长的任务拆分为三个阶段的执行模式示意图
通过上面的例子可以看到，我们能够将需要一次运行时间长的函数拆分成多个阶段实现，使得周期执行循环中，执行该函数的时间降低到各个阶段的执行时间。该模式适用于可以分解为多个执行阶段的函数，并且各个执行阶段相对独立，相邻执行阶段的运行时间可以不连续。
·时间格点采样
对于传感器数据读取，当相邻两次调用读取数据的函数的间隔小于采样间隔时，通过代码清单2-4给出的伪代码实现采样时间调整，以确保每次数据读取在特定的时间格点上执行。
代码清单2-4 基于时间格点的采样过程的伪代码
void传感器采样数据读取()
{
while ((读入当前时间() % 采样时间间隔)>0)
延迟等待固定时间();
读取并保存数据();
}
上面的伪代码中，while循环所依赖的条件“读入当前时间()%采样时间间隔”用于确定当前时间是否在采样时间格点上，如果不在，则延迟等待固定时间。
值得注意的是，在很多情况下，采样时间间隔可以通过外部数据采样硬件确定，用户程序通过FIFO读取数据，当数据尚未到达时，FIFO为空，这使得我们可以基于FIFO是否为空来判断时间格点上的数据采样是否完成，如代码清单2-5所示。
代码清单2-5 基于FIFO读取的数据采样的伪代码
void传感器采样数据读取()
{
while (外部FIFO空())
延迟等待固定时间();
读取并保存数据();
}
上述模式只有当两次执行“传感器采样数据读取”的间隔小于采样间隔时才有效，如果大于采样间隔，则会导致采样点“遗漏”，如图2-5所示。
图2-5 调用“传感器采样数据读取”函数间隔大于采样周期导致的采样缺失问题示意图
图2-5中虚线是采样时间格点，在中间由于“计算机控制输出”任务运行时间过长，出现过一次采样缺失现象。下面将讨论另外几种编程模式以避免这一问题。
2.1.2 基于中断的前后台运行模式
前面基于循环调用的编程模式难以确保各个任务的等时间周期运行，并且把实时性要求高的任务和实时性要求低的任务混在一起，难以同时兼顾不同任务实时性要求的差异。为了能够确保特定任务的实时性，我们通常使用基于中断的前后台编程模式，如图2-6所示。
其中，时钟中断服务程序负责从外部传感器得到数据输入，由于时钟中断信号具有严格的周期性，因此可以确保数据间隔的稳定，并且由于中断服务程序随时能够打断前台的数据处理、用户输入、显示输出程序执行，因此不受它们的运行速度影响。
使用前后台编程模式时，需要考虑中断服务程序的运行效率，对于前面介绍的固定时间间隔的数据采样任务，要求中断服务程序能够在采样间隔内完成，因此需要尽可能地提高运行效率，在中断服务程序中只保留必不可少的代码，尽可能将那些对实时性要求不高的数据操作放到前台应用程序中执行。比如来自传感器的数据经过了压缩和纠错编码，在中断服务程序运行期间执行数据解压缩和纠错可能导致运行时间过长，进而导致错过下一个采样的时间点，为避免这一现象出现，可以将这些运行时间长但又没有必要在中断服务程序内完成的任务转移到“数据计算及控制输出”函数中执行。
图2-6 基于中断的前后台运行模式
使用前后台编程模式的一个极端示例是将所有操作移入中断服务程序。硬件提供不同重复间隔和优先级的时钟中断，这样主程序会很简单，只有一个无限循环，而所有的工作由中断服务程序完成，如图2-7所示。
图2-7 全部基于中断运行模式
基于中断服务程序实现所有操作的方案中需要编程者处理数据访问冲突问题，比如传感器数据读取中断服务程序有可能打断“数据计算及控制输出”中断服务程序，“篡改”还未处理完的数据，这就需要可靠的数据交换机制，以确保不同的代码能够安全地访问共享数据交换空间。比如通过图2-8所示的缓冲器交换数据。
图2-8 利用环形缓冲器实现数据交换
“传感器采样数据读取”程序负责将读取的数据加入环形缓冲器，并修改缓冲器的“写指针”，而“数据计算及控制输出”程序负责从环形缓冲器读取数据，每次读完数据就修改缓冲器的“读指针”，避免对相同的共享数据区同时读写。
2.1.3 基于事件队列的运行模式
这一运行模式是基于中断的前后台运行模式的进一步改进。中断服务程序可以是传感器数据读取程序，也可以是用户指令输入检测程序。后台程序由一系列中断服务程序构成，负责读取数据供前台程序处理，还负责根据输入内容生成“事件”，将待处理事件加上数据放入固定格式的“事件数据包”，最后将事件数据包挂到全局的“事件队列”中去。前台程序则不断检查事件队列，找到其中待处理的事件数据包，提取并执行，如图2-9所示。
事件队列可以以链表形式或者环形缓冲器形式实现，由于它是由中断服务程序和前台处理程序共享的，因此在编程时需要避免访问冲突。
基于事件列表运行模式的一种扩充是在系统中构建多个事件队列，每个事件队列具有对应的优先级，前台程序根据事件队列的优先级优先处理高优先级的事件，这样能够使对实时性要求高的事件（比如数据计算）能够及时得到处理，如图2-10所示。
这一编程模式在实际应用中要求程序员考虑“低优先级事件积压”问题，即高优先级事件不断出现，导致低优先级事件长时间得不到处理、处理延迟过长的问题。需要程序员采取特定措施防止低优先级事件长时间得不到处理，一个解决方案是规定一个低优先级事件队列长度门限，当积压的低优先级事件数量超过门限时，临时提升该队列的处理优先级，避免低优先级事件队列内事件被过度积压。
图2-9 基于事件队列的运行模式
图2-10 基于多重事件队列的运行模式，每个事件队列具有不同优先级
2.1.4 带时间信息的事件队列运行模式
对于需要在特定时间执行指定动作的应用，比如有特定时序要求的硬件设备控制和网络通信协议等，可以使用带时间信息的事件队列的设计模式（见图2-11）。这一设计模式和带优先级的事件队列运行模式相似，差别是这里的每个事件队列的优先级被执行时间所取代，系统根据当前时间检查事件队列，对于时间匹配的事件队列，执行队列中的所有事件，执行完删除该队列。
图2-11 带时间信息的事件队列的运行模式
事件队列可以由当前正在执行的事件生成，也可以由中断服务程序生成。代码清单2-6中显示了这一模式的运行过程。
代码清单2-6 带时间的事件队列的运行模式的伪代码
void main()
{
while (1)
{
if (事件队列链表中最近的执行时间小于等于当前时间)
{
while (获取当前时间待执行事件())
执行事件任务();
}
}
}
2.1.5 计算图运行模式
基于计算图的运行模式在机器学习领域得到广泛应用，多种主流的神经网络框架使用计算图来描述神经网络架构。比如图2-12a给出了某一神经网络的计算图。
图2-12中每个圆圈代表网络中的一个运算，而连接各个运算的箭头线上的数字代表运算输出的张量尺寸。该网络的输入是尺寸为100×100的RGB图片，输出128维的目标分类向量和一个10×10的目标位置信息图。为后面分析方便，我们用字母表示上面6个运算，如图2-12b所示。
图2-12 某一神经网络的计算图结构。图a为计算图结构和输入输出张量尺寸；图b为计算图中各个节点的字母表示
使用计算图进行计算时，需要根据各个运算之间的数据依赖关系制定执行顺序，即每个运算节点只有当它的所有输入数据都已经准备好之后才能执行。比如图中“拼接卷积”（D）运算依赖于两个输入端的卷积B和C的运算结果，需要在这前两个卷积全部完成后才能够执行。
对于一个计算图，可以通过下面的步骤获得运算顺序。
第一步：根据ASAP（As Soon As Possible，尽可能早）策略或者ALAP（As Late As Possible，尽可能晚）策略构建计算顺序，将运算分配到各个时间段。
第二步：根据处理器的并行运算能力将每个运算节点分配到不同的时间行执行。
我们首先说明第一步，这一步骤为每个运算分配可以执行的时间段，我们分别给出ASAP和ALAP两种时间分配策略。
·ASAP时间分配策略
这一策略自顶而下，寻找做好准备的运算，并将其分配到最早的可以运行的时间段。比如一开始运算A和B都能够执行，将它们分配到第一时段，接着发现运算C可以执行（由于D依赖于C，D不可以执行），将C分配到第二时段，然后看到D和E可以执行了，将它们分配到第三时段，最后剩下F，分配到第四时段。这一分配结果如图2-13a所示。
·ALAP时间分配策略
这一策略自下而上分配时段，首先设置最后一个时段的输出，并根据这些输出确定倒数第二个时段必须执行的运算，接着再倒退一个时段，根据倒数第二个时段的运算确定在倒数第三个时段必须完成的运算，依次类推，直到所有运算都被包括进来。这一分配结果如图2-13b所示。
图2-13 ASAP和ALAP方案对每个运算的时间分配结果
从上面两种方式得到的运算中，会出现某一个时段需要执行2个运算的情况，比如图2-13的ASAP的执行时间分配算法中，第一时段需要同时运行A和B两个运算，当我们只有一个处理器时，需要对这种情况进行拆分。将同时出现多个运算的时段进行拆分，可以使每个时段只包括一个运算，图2-14给出了ASAP时间分配算法的拆分结果。
完成上述运行时间分配之后，就能够获得运行这个图的执行流程了，参见代码清单2-7中的伪代码。
代码清单2-7 计算图运行模式的伪代码
while (1)
{
得到图像数据();
执行卷积A();
执行卷积B();
执行卷积C();
执行数据拼接和卷积D();
执行全连接层E();
执行卷积层F();
输出结果();
}
图2-14 将安排在同一时段的多个运算（基于ASAP策略）分到不同时段运行。图a为原先的ASAP分配方案，时段1和时段3分别有两个运算；图b为经过了时间拆分的运算分配方案，每个时段只有一个运算
类似地，ALAP时间分配的结果也需要分割，如图2-15所示。
图2-15 将安排在同一时段的多个运算（基于ALAP策略）分到不同时段运行。图a为原先的ALAP分配方案，时段2和时段4有两个运算；图b为经过了时间拆分的运算分配方案，每个时段只有一个运算
对于内存有限的嵌入式系统，上述运算过程可以根据内存进行优化，不同运算顺序所需的最大内存占用量不同，比如上面ALAP时间分配算法的原始结果中，第二时段的运算B和C所处的运算时间有两种“分裂”方式——先执行B运算或者先执行C运算，两种方式的选择是任意的。但不同的执行顺序可能带来不同的内存需求，因此执行顺序的选择是值得优化的。下面分别画出这两种时间分配方法，以及对应的系统最大数据存储量，如图2-16所示。
图2-16 两种运算时间分配方案对应的系统最大内存量。图a为时间段2执行运算C时，运算期间内存占用量；图b为时间段2执行运算B时，运算期间内存占用量
图2-16中写出了每个运算输出的数据尺寸，以及执行运算过程中最大的需要暂存的数据总量，可以看到两种方案对应不同的数据暂存总量。带下划线且用粗体标注出的为运算期间的最高数据暂存总量，其中图2-16a中为55 600，图2-16b中为51 200。可见不同的运算顺序带来不同的内存要求，对于存储空间有限的嵌入式系统，有必要仔细检查不同执行顺序的内存占用量，选择内存需求最低的执行顺序。
2.2 通用软件优化方法
在这一节我们介绍下面几类通用优化方法，分别是：
·循环结构优化——针对具有多层嵌套的循环执行结构的优化。
·时间空间平衡——增加内存空间换取运算时间降低或者反过来的优化方法。
·运算精度和性能平衡——通过降低运算精度来提升运算速度并降低存储空间的方法。
·底层运算的快速实现算法——对于底层的乘除运算用等效的位操作提速。
·内存使用优化——通过调整运算次序或者复用内存，降低程序对嵌入式系统的最大内存需求。
这几类优化方法除了应用于嵌入式机器学习领域外，也能用于其他嵌入式应用，属于相对通用的优化方法。
2.2.1 循环结构优化
循环结构优化可以通过对比代码清单2-8中给出的三段功能相同的代码来说明，其中第一段和第二段代码对应的循环结构最内层的运算次数相同，都是109次，但对于第一种写法calc1，最内层的循环变量k需要重复加载108次，而第二种写法calc2的循环变量k只需要重复加载4次，降低了循环变量k上的重复操作次数，提高了运算效率。一般而言，如果嵌套循环的循环变量次序可以更换的话，我们会倾向于将循环次数少的循环变量放置到外层。
在calc3中，最内层的循环k被展开，程序行数增加了，但完全消除了额外的循环变量k的运算，因此运行效率会有所提升。
代码清单2-8 不同循环结构的代码示例
void calc1()
{
for (unsigned i=0; i<10000; i++)
for (unsigned j=0; j<10000; j++)
for (unsigned k=0; k<4; k++)
data_proc(i,j,k);
}
void calc2()
{
for (unsigned k=0; k<4; k++)
for (unsigned i=0; i<10000; i++)
for (unsigned j=0; j<10000; j++)
data_proc(i,j,k);
}
void calc3()
{
for (unsigned i=0; i<10000; i++)
for (unsigned j=0; j<10000; j++)
{
data_proc(i,j,0);
data_proc(i,j,1);
data_proc(i,j,2);
data_proc(i,j,3);
}
}
下面再看一下代码清单2-9所列出的矩阵乘法代码中循环结构的优化，考虑不同循环结构的两种矩阵乘法实现代码。
代码清单2-9 不同循环结构的矩阵乘法实现代码
void mat_mul1(int **a, int **b, int **c)
{
for (int i=0; i<100; i++)
for (int j=0; j<100; j++)
for (int k=0; k<100; k++)
c[i][j] = c[i][j] + a[i][k] * b[k][j];
}
void mat_mul2(int **a, int **b, int **c)
{
for(int i=0; i<100; i++)
for(int k=0; k<100; k++)
for(int j=0; j<100; j++)
c[i][j] = c[i][j] + a[i][k] * b[k][j];
}
如果观察最内层循环对应的数组访问，可以看到函数mat_mul1有一个连续内存访问，而mat_mul2为两个。函数mat_mul1的最内层循环变量k使得对数据a[i][k]的访问是连续的，但对数据b[k][j]的访问是跳跃的，如图2-17所示。
图2-17 函数mat_mul1的最内层循环对数据b[k][j]的访问模式
函数mat_mul2改变了循环结构，使得最内层循环变量j涉及的两个数据访问都是连续的，如图2-18所示。
图2-18 函数mat_mul2的最内层循环对数据b[k][j]和c[i][j]的访问模式
即循环变量j使得对b[k][j]和c[i][j]的访问都是连续的。在高性能嵌入式系统中，通常对于跳跃的内存的访问效率比连续内存的低，这是因为CPU内核一般通过Cache访问数据。而Cache通过内存控制器间歇地从外部存储器读取地址连续的一批数据放入Cache。对于连续地址访问，CPU通过一次Cache填充就能够得到所需要访问的多个连续地址数据，相比之下，如果跳跃地访问，CPU需要频繁地启动数据加载硬件填充Cache，运行效率低。
2.2.2 时间空间平衡
在嵌入式系统中，由于成本和功耗约束，对软件的运算速度和内存也有约束，实际编程过程中可以通过时间换取空间或者空间换取时间的方式进行优化。考虑下面的代码，它计算输入数据二进制形式下的1的个数，可以使用代码清单2-10列出的几种方案实现。
代码清单2-10 统计二进制数据x中位1的个数的几种代码实现
int count1(unsigned long x)
{
int c=0;
for (unsigned i=0; i<32; i++)
if (x&(1UL<<i))
c++;
return c;
}
int count2(unsigned long x)
{
int c=0;
int lut[]={0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
for (unsigned i=0; i<8; i++)
{
c+=lut[x&0xf];
x>>=4;
}
return c;
}
int count3(unsigned long x)
{
x = (x&0x55555555) + ((x>>1)&0x55555555);
x = (x&0x33333333) + ((x>>2)&0x33333333);
x = (x&0x0f0f0f0f) + ((x>>4)&0x0f0f0f0f);
x = (x&0x00ff00ff) + ((x>>8)&0x00ff00ff);
x = (x&0x0000ffff) + ((x>>16)&0x0000ffff);
return (int)x;
}
其中算法count2通过32次循环，对1的个数累加得到答案，而count2通过8次循环得到结果，其中每一次分析x中的4个位，而这4个位中1的数量通过查表lut直接获得。函数count2的代码尺寸和存储空间增加了，但换取了运算速度的提升。
代码片段count3还给出了基于“并行位相加”的方案，该方案利用CPU内部的32位无符号整数加法器的特性进行计算，其中第一行
x = (x&0x55555555) + ((x>>1)&0x55555555);
实现x中序号为奇数的位和序号为偶数的位相加，结果依次存在连续2位空间，如图2-19所示。
图2-19 每两个相邻位相加，运算结果存放于相邻的两个位中
接下来的代码
x = (x&0x33333333) + ((x>>2)&0x33333333);
实现对相邻的两组2位数据相加，结果依次存于连续4位空间，如图2-20所示。
图2-20 每相邻两组2位数据相加，运算结果依次存放于连续4位空间中
后续的代码
x = (x&0x0f0f0f0f) + ((x>>4)&0x0f0f0f0f); x = (x&0x00ff00ff) + ((x>>8)&0x00ff00ff); x = (x&0x0000ffff) + ((x>>16)&0x0000ffff);
分别实现：相邻的两组4位数据相加，结果依次存于连续8位的空间；相邻的两组8位数据相加，结果依次存于连续16位的空间；相邻的两组16位数据相加，结果存于连续32位的空间，如图2-21所示。
图2-21 通过邻两组数据相加统计数据中非零位数目的示意图
上述运算结构巧妙地应用了CPU内部高数据位宽的加法器的特性实现并行加法，快速得到位1的计数结果。
2.2.3 运算精度和性能平衡
不少对速度要求很高的应用，允许牺牲运算精度来换取速度的提升。比如代码清单2-11中通过查表计算三角函数cos的近似结果。
代码清单2-11 基于查表得到余弦函数的实现代码
#define _USE_MATH_DEFINES
#include <math.h>
const float LUT[]=
{
1.0000f, 0.9848f, 0.9396f, 0.8660f,
0.7660f, 0.6427f, 0.5000f, 0.3420f,
0.1736f, 0.0000f, -0.1736f,-0.3420f,
-0.5000f,-0.6427f, -0.7660f,-0.8660f,
-0.9396f,-0.9848f, -1.0000f,-0.9848f,
-0.9396f,-0.8660f, -0.7660f,-0.6427f,
-0.5000f,-0.3420f, -0.1736f,-0.0000f,
0.1736f, 0.3420f, 0.5000f, 0.6427f,
0.7660f, 0.8660f, 0.9396f, 0.9848f, 1.0f
};
float lut_cos_deg(unsigned x)
{
x%=360;
unsigned u=x/10,v=x%10;
float a=LUT[u], b=LUT[u+1];
return (float)(a+(b-a)*((float)v)/10.0f);
}
函数lut_cos_deg计算输入数据x的余弦值，其中x是单位为度（°）的非负整数。将通过查表和线性插值结合得到的结果与进行余弦精确计算的结果的差别小于0.004。图2-22给出了精确计算余弦值和上述近似运算的结果的差别，两种运算方式对应的曲线几乎重合。
下面再给出双曲正切（tanh）运算的近似算法。tanh运算在神经网络的激活函数中会被用到，它的定义为
它和另一种常用的激活函数sigmoid关系紧密，sigmoid函数定义为
图2-22 精确计算的余弦值和使用查表近似算法计算的余弦值之间的差别比较。图a为两个运算结果的比较；图b为计算结果的差的曲线
tanh和sigmoid函数的关系为
传统的tanh运算需要进行至少一次指数运算，运算效率低，但我们可以用代码清单2-12给出的近似算法实现。
代码清单2-12 基于近似公式计算tanh的C程序代码
#include <math.h>
float tanh(float x)
{
if (x>(float) 3.4f) return 1;
if (x<(float)-3.4f) return -1;
x*= 1.0f/3.4f;
x*= fabsf(x)-2.0f;
return x*(fabsf(x)-2.0f);
}
这一运算需要进行3次乘法，运算量低于精确计算tanh，但付出的代价是运算结果有最大0.018的误差。图2-23是精确计算的tanh和上面的近似算法得到的tanh的曲线对比，可见两者相差无几。
最后我们介绍以2为底的浮点数对数运算近似算法[1]。由于浮点数本身是以指数加上尾数的形式存储的，对它取以2为底的对数可以直接用到浮点数表示形式的特点。下面是IEEE754格式的单精度浮点数的格式说明。
图2-23 精确计算的tanh和使用近似算法计算的tanh之间的差别比较。图a为两个运算结果的比较；图b为计算结果的差的曲线
图2-24中的指数位E可以看作无符号整数，上述格式的单精度浮点数代表的数值为
图2-24 IEEE754单精度浮点数格式示意图
对于正浮点数，s=0，代表的数值简化为
对该正浮点数求对数，并应用近似得到
根据上述近似可以给出求单精度正浮点数的近似算法，如代码清单2-13所示。
代码清单2-13 单精度浮点数的近似对数计算
float log2_approx(float x)
{
unsigned long v = *(unsigned long*)&x;
return (float)((v >> 23) & 0xFF)-127.0 + \
(float)(v & 0x7FFFFF) / (float)0x800000;
}
我们在应用中除了计算2为底的对数外，也会计算其他数为底的对数，比如计算，这可以通过实现，即使用给出的近似计算算法得到log2x后再乘以常数即可。
2.2.4 底层运算的快速实现算法
对于很多整数运算，可以通过其二进制的位结构的特点来优化计算。当嵌入式处理器没有高效乘法和除法硬件时，我们可以用加减法以及移位实现很多运算。比如考虑计算下面的乘法运算：
其中x和y都是整数（C语言的long类型），这一运算可以通过代码清单2-14给出的程序片段实现。
代码清单2-14 整数乘法运算的快速实现程序
long b=(x<<3)-x; // b=x*(0111)2 y=(b<<4)+b;
如上所示代码中，乘数(01110111)2被拆成相同的两部分(0111)2，因此计算出b=x×(0111)2后，通过让程序中的b左移4位再和它自己相加，就得到结果了。而b=x×(0111)2可以根据(0111)2=(1000)2-1简化，即
上述过程将乘法分解为移位相加运算。对于特定常数的运算，通过观察可以得到最优的分解方案，也可以利用软件搜索最优分解方案。这一部分内容会在第4章详细介绍。
对于除法运算，我们同样可以使用上述优化方案，这里的关键是“化除为乘”，即把计算x/y的运算改成。比如计算x/11，我们选取N=8，可以得到
当x为整数类型时，可以用
计算，进一步用移位取代乘除法得到下面的C代码片段：
((x<<4)+(x<<3)-x)>>8
上述运算的精度和选取的位宽N有关（比如在式（2-9）中N=8），N的取值越大，精度越高，但相应的运算越复杂（将乘法转成移位加减后，运算量可能增加）。
下面考虑浮点数和整数的乘法。由于数据移位运算不能直接应用在浮点数上，因此在浮点数和整数的乘法中，需要将移位操作用数据相加实现，比如计算y=x×(01110111)2，但这次x是浮点数，运算过程需要修改成代码清单2-15所示的形式。
代码清单2-15 浮点数和整数相乘的快速实现程序
float b,c; b=x+x; b+=b; b+=b; // b=x*(1000)2 b-=x; // b=x*(0111)2 c=b+b; c+=c; c+=c; c+=c; // c=b*(10000)2 y=c+b;
上述移位操作转成数据相加操作，对于左移n位需要数据反复和自己相加n次，当n的取值较大时，运算效率低。另一个可选的方案是根据浮点数格式的特点，直接将移位值作用到浮点数的指数域上。对于IEEE754标准定义的单精度浮点数（即C语言的float类型）由32位数据表示，其中第23～30位存放浮点数的指数部分，对这个区域加或者减n等效为乘以或者除以2n。根据这一特性可以得到代码清单2-16所示的计算浮点数乘以整数(01110111)2=119的代码。
代码清单2-16 基于IEEE754数据格式实现浮点数和整数快速相乘的代码
// 计算x*119，其中119的二进制形式为01110111
float mul_119(float x)
{
float y, b = x;
*(unsigned long*)&b += (3 << 23); // b=(x<<3)
b -= x; // b=(x<<3)-x
y = b;
*(unsigned long*)&y += (4 << 23); // y=(b<<4)
return y + b; // 返回(b<<4)+b;
}
在没有浮点乘法硬件的CPU平台上，上述代码能够快速实现浮点数乘法，但使用这一代码时需要注意数据溢出问题，实际应用时需要补充异常检测和处理过程。读者可以结合第4章内容理解上述代码的含义。
我们最后给出用加法实现浮点数乘法的近似算法。这一算法用了上一节介绍的浮点数近似对数计算的思想。首先，两个单精度浮点数相乘可以写成下面的形式：
x和y的浮点数表示格式中，指数和尾数域的内容分别用和表示，根据上一节近似对数计算方法，有
于是
令，并且z也用单精度浮点数表示，即，于是有
对比式（2-13）和式（2-14）得到
注意
上面的公式隐含假设，如果该条件不满足，则会产生额外的误差。
基于上面给出的z的指数和尾数，可以根据单精度浮点数的二进制格式重新拼接成为单精度浮点数z，整个运算过程通过代码清单2-17给出的程序实现。
代码清单2-17 单精度浮点数的近似乘法运算
float mul_approx(float x, float y)
{
long z=*(long*)&x + *(long*)&y -(127 << 23);
return *(float*)& z;
}
上述代码中隐含了多种运算技巧：首先，式（2-15）中浮点数x和y的指数和尾数相加可以用两个浮点数的二进制形式直接当成整数相加实现；其次，上述运算考虑了乘法结果的符号位可以直接通过两个被乘数的符号位直接相加并丢弃进位得到。代码中-(127<<23)对应式（2-15）中减去127的运算，考虑单精度浮点数的指数位在第23～30位，因此将其左移23位。需要注意的是，在实际应用上述代码时，还应补充额外代码来检测数据溢出问题并进行处理。
2.2.5 内存使用优化
1.缓冲区分配
嵌入式系统只有有限的内存空间，需要通过合理地复用内存来降低存储压力。比如考虑之前的计算图，我们之前计算了它在计算过程中的最大数据暂存量，在实现过程中，需要通过存储器复用策略达到最低内存占用量。
回顾图2-12，根据每个运算环节的输入和输出数据（张量）尺寸，可以看到总共有下面几种尺寸的张量：100×100×3，50×50×8，20×20×8，20×20×64，128，10×10。但在实际运行时，我们可以只分配3块数据缓存区：
·1号数据缓冲区——存储尺寸是25 600（20×20×64）。
·2号数据缓冲区——存储尺寸是30 000（100×100×3）。
·3号数据缓冲区——存储尺寸是3200（20×20×8）。
图2-25中给出了特定运算时间分配方案下的缓冲区的分配方案。
图2-25 计算图计算期间的数据缓冲区分配示意图
图2-25中灰色底色的矩形块代表数据缓冲区，里面数字分别对应1、2、3号缓冲区。虽然各个运算需要的输出缓冲区有大有小，但只要分配时确保输出数据小于缓冲区尺寸即可保证算法正常运行。上述分配方案需要对应的数据缓冲区的总存储量是58 800个数据。代码清单2-18给出了上述执行流程的伪代码。
代码清单2-18 计算图执行过程中缓冲区复用的伪代码
float buffer_1[25600];
float buffer_2[30000];
float buffer_3[3200];
while (1)
{
得到图像数据(buffer_2); // 输入数据存储于buffer_1
执行卷积A(buffer_2, buffer_1); // 输入buffer_2，输出buffer_1
执行卷积B(buffer_2, buffer_3); // 输入buffer_2，输出buffer_3
执行卷积C(buffer_1, buffer_2); // 输入buffer_1，输出buffer_2
执行数据拼接和卷积D (buffer_3, buffer_1); // 输入buffer_3，
// 输出buffer_1
执行全连接层E(buffer_2, buffer_3); // 输入buffer_2，输出buffer_3
执行卷积层F(buffer_1, buffer_2); // 输入buffer_1，输出buffer_2
输出结果();
}
2.原址运算
前面讨论的缓冲区内存分配是假设了每个运算执行期间，输入和输出数据缓存区不能相同，因此每个运算占用的空间是输入和输出数据缓存区的总和，对于内存受限的嵌入式系统，可以通过将部分运算改成“原址计算”的模式进一步降低内存占用量，比如比较代码清单2-19所列出的两个运算（神经网络的ReLU运算）的实现方案的内存使用情况。
代码清单2-19 不使用原址计算和使用原址计算的ReLU激活函数代码
void relu1(float *in, float *out, int size)
{ for (int n=0; n<size; n++) out[n]=max(in[n],0); }
void relu2(float *inout, int size)
{ for (int n=0; n<size; n++) inout[n]=max(inout[n],0); }
代码中，函数relu1的输入和输出是分离的两个缓存区in和out，而函数relu2的输入和输出使用相同的缓冲区inout。可见relu2在运算期间的内存占用总量降低到原先的50%。
3.内存“折叠”复用
原址运算并不能适用于所有的运算环节，但有一大类运算可以部分实现原址计算。下面以卷积神经网络中常用的深度可分离卷积（depthwise separable convolution）为例进行说明，该卷积运算作用于C个输入的二维特征图（矩阵），对输入特征图逐一进行卷积得到对应的C个输出特征图。比如考虑在输入3个特征图的情况下，图2-26a是输入和输出数据存放的内存分布，图2-26b是原始的内存访问模式，图2-26c是改进的内存访问模式。
图2-26 进行深度可分离卷积计算时内存访问模式示意图。图a为原始内存访问模式下占用的总内存；图b为运算的三个阶段实际占用的内存；图c为通过内存复用，将写数据“折叠”到低地址，降低总的内存占用
在图2-26b的内存访问模式中，灰色对应需要保留的内存空间，虚线对应暂时不需要的内存空间。可见在实际运算期间，低地址内存可以依次释放，而高地址内存被依次占用，每一时刻需要占用的内存量小于输入和输出占用的内存之和。根据这一特性，我们可以按图2-26c所示模式复用内存。图2-26c和图2-26b相比，可见由于我们在计算的第2阶段和第3阶段复用了之前释放的低地址内存，因此需要的总内存量小于图2-26a所示的内存访问方案。
上面给出的内存优化方式可以看作“地址折叠”访问方式的一个例子，即输出数据存储时，若写地址超出分配的内存空间上限，则折返到低位地址。代码清单2-20为地址折叠缓存区复用的代码片段，显示了基于地址折叠方式计算并存放输出数据的具体编程模式。
代码清单2-20 地址折叠缓存区复用的代码示例
// buffer是缓冲区，存放输入和输出数据， buffer[0]是第1个输入数据
// offset是运算结果的存放位置，即buffer[offset]对应第1个输出数据
// buffer_size是缓冲区buffer的总尺寸
void calc(float* buffer, int buffer_size)
{
float result;
// 循环处理每个输入数据
for (int n=0; n<N; n++)
{
// 运行数据处理
result=data_processing(buffer+n);
// 存放处理结果，注意数组下标用%buffer_size运算实现存储地址"折叠"
buffer[(offset+n)%buffer_size]=result;
}
}
上述程序运行前后的内存结构如图2-27所示。
图2-27 内存“折叠”访问模式下程序运行前后的缓冲区占用情况示意图。图a为程序运行前输入数据连续存放于缓存区“底部”；图b为程序运行后，输出数据存放位置
由图2-27可以看到缓冲区内部数据的存放位置。在程序运行期间，输出数据由offset依次向上存放，到达缓存区边界后“折返”到缓存器底部。这一访问模式的缺点是输出数据有可能被“切分”成两段，进行后续运算时要根据这一结构访问数据。
图2-28给出能够使用“内存折叠”的代码的内存访问模式的特点。如图2-28a所示，阴影部分代表程序运行期间使用的内存地址范围。程序执行期间，内存占用区间不断改变，呈现图中所示的“斜条带”状，对于这一类代码，内存按时间分配使用方式可以改成图2-28b所示模式，运算执行到后面阶段时内存访问被“折叠”到之前空闲下来的空间。
图2-28 内存“折叠”访问模式示意图，阴影部分代表程序运行期间使用的内存地址范围。图a为优化前的内存占用随时间变化的规律；图b为通过地址“折叠”降低内存占用
2.3 小结
本章首先介绍了几种编程模式，每种编程模式的实用性和具体的应用相关，并不一定总是选择最复杂的方案。嵌入式编程是对运行实时性、开发复杂度、CPU使用效率等多方面指标的平衡，并没有固定的或唯一的编程方式。另外，嵌入式系统的编程模式和系统复杂度有关，对于运行Linux或者Windows 10的系统，计算资源丰富，各种计算机编程框架和编程模式都能够使用。但对于资源受限、运算能力较弱的系统，需要我们精细地分配CPU运行时间，在没有操作系统的条件下完成任务。在实际应用中，还有另一种方案，就是通过现有各种实时操作系统（比如开源的FreeRTOS等）帮助我们实现多任务并发运行，并帮助管理外部设备实现共享。本书不重点关注这部分内容，对于嵌入式实时操作系统也不再展开介绍，感兴趣的读者可以阅读相关材料。
本章还介绍了通用软件优化方法，这些优化方法能够用于提高各种嵌入式应用软件的运行效率，但值得注意的是，这些优化方法在不同的硬件平台上的效果是不同的，比如对于支持硬件浮点指令的处理器，就没有必要用加减法简化浮点计算；对于使用数据Cache的处理器，不同的程序循环结构的性能受Cache行为的影响。对于每一种优化方法，还需要读者仔细对比使用前后软件性能的实际提升量，以确保所使用的优化方法有效果。
参考文献
[1] Mitchell J N.Computer Multiplication and Division Using Binary Logarithms[J].IRE Transactions on Electronic Computers,1962,EC-11（4）：512-517.