OpenCL Workshop 1 —— 数字音频滤波

Introduction

这两年深度学习大火，Cuda跟着吃红利，OpenCL发展也很快。虽然OpenCL不是事实上的标准，但是作为开放标准，适应性是很强的，除了显卡之外，CPU/FPGA上都可以执行。

第一个OpenCL Workshop的具体目标就是编写一个音频文件升频工具，用来给PCM编码的WAV文件四倍频，把升频结果写到一个新的WAV文件里面。

用于升频的主要方法，数字滤波，可以广泛用于数字音频的处理。

首先会用传统的C语言编写单线程升频算法，然后用OpenCL编写并行加速版本，最后用CodeXL比较两者性能差距。

OpenCL Fundamental

首先简单介绍一下OpenCL并行开发的基本思路。

• OpenCL总体设计

OpenCL程序是分为Host和Device两端的，Host一般来说就是CPU，而Device可以是CPU/GPU/FPGA/ASIC加速器等。

Host的功能主要是调用CL提供的API，对Device进行配置、命令，以及内存操作，本质上就是一个可执行程序；Device负责执行CL Kernel核函数，进行具体的计算，Kernel一般是在运行之前编译以保证设备无关性（但是也可以针对设备预先编译，提高执行效率）。

Kernel是用OpenCL Kernel C语言编写的（2.2版OpenCL带来了OpenCL C++ Kernel语言），基本上语法就和C语言一致，但是算法上要根据设备的不同有所调整（否则性能得不到优化），这块在Workshop实作中会有所体现。

• 并行计算体系

OpenCL的最小计算单元是Processing Element（PE），对应CPU的线程、N卡的Cuda Core、A卡的ALU；多个PE组成一个Compute Unit（CU），对应CPU的Core、N卡的Stream Multiprocessor、A卡的Compute Unit（对，A卡这一层就叫CU）。

PE和CU是OpenCL中相当重要的两个概念，CU是指令单位，而PE是数据单位，由于CU是由多个PE组成的，这也就构成了SIMD（单指令多数据）的并行架构。

在PE上执行的程序单元叫做work item，在CU上执行的叫做work group，不难看出work group是由work item组成的。

这个work item，其实就是执行起来的核函数，所以核函数其实就类似于多线程编程中的线程函数，对于有HT（超线程）能力的CPU来说，一个核心可以同时执行2个线程（但是我们可以分配很多个线程给核心，CPU自己会调度）；而对于GPU来说，一个CU可以同时执行N个work item，这个N取决于CU里有多少个PE（同样我们可以分配很多个work item给work group，CL自己会调度）。如同线程有线程号一样，work item也有自己的全局编号，从0到work group size - 1，实际上有的时候work item就叫做thread。

总结一下对应关系：

work group ->  CU -> CPU Core / AMD CU / NV SM

Kernel -> work item -> PE -> CPU Thread / AMD ALU / NV Cuda Core

可以看到CU和PE这两个抽象概念，把具体的软硬件概念粘合起来了。

DSP Background

音频升频的方法，属于数字信号处理（DSP）的范畴。升频的主要原因，是为了把DAC（数模转换）的量化噪音提高到更高的频率，这样对于DAC之后的模拟低通滤波器的要求就低得多（因为有更大的频率范围来滚降，滤波器级数可以比较小）。

假设我们有一段44.1kHz/16bit的音频数据，做4倍频的方法：

1. 给每个采样之前插入3个0，得到176.4kHz采样率的数据
2. 用一个转折频率20kHz的LPF（低通滤波器）滤除插0带来的高频信号

上面的方法可行是因为，插0不会在原始音频带内插入信号，故而可以用LPF滤掉插入的高频信号，即可得到原来的音频带信号（但是采样频率变成了4倍）。用傅立叶变换查看频谱即可证明，请看官老爷自行实验。

对信号做低通滤波，频域上描述很简单，就是把经过快速傅立叶变换（FFT）的频谱信号，转折频率以上的频谱系数直接乘0，再做傅立叶逆变换得到时域信号。但是实际上，我们不会在频域上变换（因为有傅立叶-傅立叶逆变换的开销）。

在时域上完成滤波，用DSP设计软件得到LPF的系数数组（物理上相当于单位冲击响应的幅度），直接与插0值后的数据做卷积，这就在时域上完成了低通滤波。

做卷积的时候有个优化，就是每个原始数据都会插入3个0，这些0与LPF相乘还是0，所以可以不计算。优化的具体方法是把LPF系数重排成一个四行矩阵，假设原来LPF系数是[w0, w1, ... w4n-1]，重排之后变成：

[
  [w0, w4, ... w4n-4],
  [w1, w5, ... w4n-3],
  [w2, w6, ... w3n-2],
  [w3, w7, ... w4n-1]
]

然后用这个4行矩阵与原始（非插值）数据构成的列向量做卷积，每个原始数据会算出4个值，这就是4倍插值后的结果了。这种矩阵叫做coeficient banks。

这样优化总的来说节约了75%的计算量。

LPF的性能越好，它的单位冲击响应就越长，滤波器系数数组也就越长。对于通带纹波0.0001dB、阻带衰减-120dB，通带最高频率20kHz、阻带最低频率22.05kHz的LPF来说，需要的系数长度为596（正好可以被4整除，方便我们生成4倍升频的coeficient banks）：

这个滤波器指标比所有DAC和DF（数字滤波）芯片内置的滤波器都要好。

我们完全可以实现通带纹波0.00001dB、阻带衰减-139dB（32bit极限）的滤波器，只需要增加系数长度（相应的计算量也会增加）。

Implementation

首先是Wav的读取和写入，这里我们限定只处理44.1kHz/16bit的双声道PCM编码WAV文件。

需要一提的是8bit的PCM数据是用无符号整数编码的，而16bit则是用二进制反码编码的。我们主要关心16bit的情况，对应的CL Host数据类型是cl_short，在Kernel中直接使用short即可。编写CL程序的时候，用cl_*数据类型可以保证数据长度始终是固定的（与32/64bit平台无关，与设备无关）。

另外WAV文件是把左右声道的数据交错放置的（左1 右1 左2 右2... 左N 右N），处理的时候需要单独取出一个声道的数据，最后保存的时候再次交错放置。

• 单线程实现

下面来看看C语言单线程的升频代码：

int16_t* overSampling(int16_t* data, int sampleCount, float** banks, int tapCount, int osRatio) {
    // add tapCount-1 0s before data
    int tempSize = sampleCount + tapCount - 1;
    int16_t* tempData = (int16_t*)malloc(tempSize * sizeof(int16_t));
    memset(tempData, 0, tempSize * sizeof(int16_t));
    memcpy(&tempData[4], data, sampleCount * sizeof(int16_t));

    int osSampleCount = sampleCount * osRatio;
    int16_t* osData = (int16_t*)malloc(osSampleCount * sizeof(int16_t));
    memset(osData, 0, osSampleCount * sizeof(int16_t));

    int bankSize = tapCount / osRatio;
    int osIndex = 0;
    for (int i = 0; i < sampleCount; i++) {
        for (int j = 0; j < osRatio; j++) {
            float* bank = banks[j];
            float sum = 0;

            for (int k = 0; k < bankSize; k++) {
                sum += bank[k] * tempData[i + k];
            }

            osData[osIndex] = sum;
            osIndex++;
        }
    }

    free(tempData);

    return osData;
}

banks就是前述的四行LPF系数矩阵，tapCount是LPF系数的个数。这段代码首先给原始数据前面补充了tapCount - 1个0，这是为了刚好可以计算得到4倍采样。然后为计算结果分配了内存。

接下来的for循环就是计算卷积的具体过程。对于每个采样，我们都与系数矩阵做了卷积，求出了4个值。

• 并行实现思路

OpenCL改造的关键，就在于如何把运算分配到PE上。我们观察上面的单线程算法，实际上对于不同的i，计算结果是互不影响的。

换个方向考虑， 这就说明每个采样的升频结果都是互不影响、互不依赖的——这其实并不一定，有DSP背景的看官老爷应该知道，这其实是FIR（有限冲击响应）型滤波器的特点。还有一种IIR（无限冲击响应）型滤波器，是要用前面计算结果累加到后面的。这里就不展开讨论，只要在设计滤波器时，选择FIR型，就可以运用以上算法实现，保证结果之间互不依赖。

总之，也就是说这里我们可以把第一层for循环展开成为并行处理的方式。用OpenCL的话说，就是把第一层循环内的代码写成Kernel，而循环里面这个i，就是前面说过的work item的global编号（Kernel里可以调用get_global_id(0)获取）。

OpenCL Kernel在GPU上的主要优化在于：

1. 全局内存的合并读写。global_id相邻的work item尽量访问相邻的全局数据，计算卷积的时候，根据global_id计算的输入数据偏移是相邻的，输出数据也是直接写入相邻的全局内存，所以CL会自动进行读写合并。
2. 局部内存避免bank conflict，计算卷积的时候是不需要LDS的，所以不存在bank conflict。
3. Kernel尽量避免分支，同样，计算卷积是不需要分支的。

所以FIR的计算是很适合在GPU上并行加速的。而IIR的计算，需要累加之前计算出来的结果，在GPU上并行实现会复杂很多。

• Host端实现

接下来我们看一下OpenCL实现中Host端的准备代码：

    cl_platform_id platform = getPlatform();

    cl_device_id device = getDevices(platform);

    cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL);

    cl_command_queue commandQueue = clCreateCommandQueue(context, device, 0, NULL);

    cl_program program = getProgram(clProgramPath, context, device);

    printClBuildingLog(program, device);

    cl_kernel kernel = clCreateKernel(program, "fir", NULL);

以上代码很清晰的展示了OpenCL Host端的准备流程：

1. 获取Platform。基本上Platform是按照各个公司划分的，Intel、NV、AMD等，对应他们各自的OpenCL实现（一般是在驱动内实现），一个系统内同时装多个Platform不会互相影响。AMD的Platform可以支持AMD、Intel的CPU作为Device（当然还有AMD自家的GPU），而Intel的Platform就只支持它自家的CPU和GPU，NV的只支持自家的GPU。
2. 获取Device。选定了Platform之后就可以获取到该平台下支持的所有设备，这里我们只选择一个设备（但是也可以选择多个设备）。
3. 创建Context。Context可以包含多个设备，并且CL会自动在Context中的Device之间分享数据，不需要用户操心。
4. 创建CommandQueue。CommandQueue是针对单个设备的，对于设备的指令就是通过CommandQueue发布的（换句话说从用户的角度看Kernel的执行是以单个Device为单位的）。
5. 获取CL Program。这里就是读取磁盘上的CL Kernel文本了，并且会把文本编译成可以在Device上执行的代码。后面我们用了一个工具函数来打印编译日志，方便查看编译过程中是否有错误等。
6. 从CL Program中抽取需要执行的Kernel。一个CL Program可能定义了多个Kernel，这里就是指定需要执行的Kernel名字。

以上代码中像获取Platform都是编写了工具函数，为了突出步骤，就不写出实现代码。具体的实现就请看官老爷自行查阅文后附带的项目源码。

前面有提到work group的概念，如果把所有的work item放进一个work group，那Device上就只会有一个CU执行，这明显效率不高。所以我们需要合理的设置work group大小，让work item尽量平均分配到CU上。同时，每个work group又不能太小，否则CU里面的PE又不能充分的并行起来。

如果对于硬件不是非常了解，或者不想针对某种硬件设置work group，那么可以不指定work group size，让CL自动分配，一般来说，CL是可以达到CU/PE的最大利用率。

接下来Host端还需要负责内存分配和指定参数，以及向CommandQueue发布命令：

    // left channel
    cl_mem banksBuffer = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(cl_float) * tapCount, clBanks, NULL);
    cl_mem leftBuffer = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(cl_short) * sampleCountPerChannel, extendedLeftData, NULL);
    cl_mem osLeftBuffer = clCreateBuffer(context, CL_MEM_WRITE_ONLY, sizeof(cl_short) * osSampleCount, NULL, NULL);

    clSetKernelArg(kernel, 0, sizeof(cl_mem), &banksBuffer);
    clSetKernelArg(kernel, 1, sizeof(cl_int), &bankSize);
    clSetKernelArg(kernel, 2, sizeof(cl_mem), &leftBuffer);
    clSetKernelArg(kernel, 3, sizeof(cl_int), &sampleCountPerChannel);
    clSetKernelArg(kernel, 4, sizeof(cl_mem), &osLeftBuffer);
    clSetKernelArg(kernel, 5, sizeof(cl_int), &osRatio);

    size_t globalWorkSize[1] = { sampleCountPerChannel };
    cl_int status = clEnqueueNDRangeKernel(commandQueue, kernel, 1, NULL, globalWorkSize, NULL, 0, NULL, NULL);

    status = clEnqueueReadBuffer(commandQueue, osLeftBuffer, CL_TRUE, 0, sizeof(cl_short) * osSampleCount, osLeft, 0, NULL, NULL);

上面的代码思路也是非常清晰的：

1. 首先使用clCreateBuffer创建了Memory Object，Memory Object会自动将Host数据搬到Device（一般是Device第一次访问数据的时候才搬）。
2. 设置Kernel参数。
3. 发布命令，这里clEnqueueNDRangeKernel并没有指定local work group size，所以OpenCL会自动分配合适的work goup大小。命令加入CommandQueue之后是异步执行的，函数会立即返回。
4. 读取计算结果。指定了block参数，clEnqueueReadBuffer会阻塞直到完成数据读取。

最后Host还需要执行资源回收，限于篇幅就不示出代码了。

• Device端实现

前面已经介绍了并行实现思路，各位看官老爷可以自己试着写一下Kernel，有两个Tips：

1. Kernel C不允许指向指针的指针，所以LPF coefficient bank需要重新组织成一维向量的形式，在Kernel中根据global id计算bank的行偏移。
2. Kernel应该尽量减少写共享的全局变量，所以单线程代码中osIndex的处理方式需要改成在Kernel中计算结果index。

以下是Kernel参考代码，折叠起来先思考一下：

__kernel void fir(__global float* banks, int bankSize,
  __global const short* samples, int sampleCount,
  __global short* results, int osRatio) {

  int i, j, k;
  int bankOffset = 0;
  float sum = 0;

  i = get_global_id(0);

  for (j = 0; j < osRatio; j++) {
    bankOffset = j * bankSize;
    sum = 0;

    for (k = 0; k < bankSize; k++) {
      sum += banks[bankOffset + k] * samples[i + k];
    }

    results[i * osRatio + j] = (short)sum;
  }
}

Performance

Platform	Device	Kernel Exec Time(ms)
AMD Accelerated Parallel Processing	AMD FX(tm)-8350 Eight-Core Processor (FX8350 with 32GB DDR3-2133)	4160.62048
AMD Accelerated Parallel Processing	Pitcairn (R9 270X with 2GB GDDR5)	49.92667

CodeXL报告GPU的Kernel Occupancy是100%，Cache命中率68%左右，没有LDS Bank Conflict。

这里CPU也是执行OpenCL Kernel，而不是单线程代码，所以可以动用到CPU内部所有核心。使用OpenCL，即使是在CPU上，也可以轻易实现并行计算。

从Kernel执行时间来看，GPU取得了80多倍的效率提升，这是符合逻辑的。270X有1280个1080MHz的ALU，FX8350有4个4.0GHz的FPU（是的，就是4个）：(1280 * 1.08) / (4 * 4.0) =  86.4。

The End

本文的完整工程代码可以从Github上获取，开发环境需要安装VS2015，AMD APP SDK 3.0。