第3A部分:使用纯C实现BFP FIR滤波器
在前面的3个阶段中,FIR滤波器是使用定点算术实现的。在接下来的三个阶段中,将使用块浮点(BFP)算术实现滤波器。
在第3A部分中,我们将使用C语言实现BFP逻辑,以确保读者能够看到发生的一切。在接下来的两个阶段中,我们将使用lib_xcore_math库中的函数来辅助我们。
来自lib_xcore_math
本页面引用了lib_xcore_math中的以下类型和操作:
实现
// 计算int32向量的头空间。
static inline
headroom_t calc_headroom(
const int32_t vec[],
const unsigned length)
{
// int32不能有超过32位的头空间。
headroom_t hr = 32;
for(int k = 0; k < length; k++)
hr = MIN(hr, HR_S32(vec[k]));
return hr;
}
如在第3部分中提到的,尾数向量的头空间定义为其元素中的最少头空间。这里hr初始化为32,因为32位整数不能有超过32位的头空间,我们遍历元素,找到最小的头空间。
HR_S32()是在lib_xcore_math中定义的一个宏,用于高效地获取有符号32位整数的头空间。它使用了本地指令CLS,该指令报告一个整数的前导符号位数。
/**
* 这是实际应用FIR滤波器的硬件线程的线程入口点。
*
* `c_audio`是用于与tile[0]交换PCM音频数据的通道。
*/
void filter_task(
chanend_t c_audio)
{
// 用作BFP向量的样本历史记录
struct {
int32_t data[HISTORY_SIZE]; // 样本数据
exponent_t exp; // 指数
headroom_t hr; // 头空间
} sample_history = {{0},-200,0};
// 用作输出帧的BFP向量
struct {
int32_t data[FRAME_SIZE]; // 样本数据
exponent_t exp; // 指数
headroom_t hr; // 头空间
} frame_output;
// 无限循环
while(1) {
// 读入一个新的帧
rx_and_merge_frame(&sample_history.data[0],
&sample_history.exp,
&sample_history.hr,
c_audio);
// 计算输出帧
filter_frame(&frame_output.data[0],
&frame_output.exp,
&frame_output.hr,
&sample_history.data[0],
sample_history.exp,
sample_history.hr);
// 在向量的前面为新的样本腾出空间。
memmove(&sample_history.data[FRAME_SIZE],
&sample_history.data[0],
TAP_COUNT * sizeof(int32_t));
// 发送处理后的帧
tx_frame(c_audio,
&frame_output.data[0],
frame_output.exp,
frame_output.hr,
FRAME_SIZE);
}
}
在这里,filter_task()的功能与之前的阶段基本相同,只是现在sample_history和frame_output是跟随尾数向量一起跟踪指数(exp)和头空间(hr)的struct。
请注意,在调用的每个函数中,还传递了指向指数和头空间的指针。这在块浮点中很常见,因为没有一个公共的对象类型来表示BFP向量(例如lib_xcore_math的bfp_s32_t)。这可能会变得繁琐,比如在调用filter_frame()时,需要传递2个BFP向量作为参数。许多操作涉及3个或更多的BFP向量。
// 计算整个输出帧
void filter_frame(
int32_t frame_out[FRAME_SIZE],
exponent_t* frame_out_exp,
headroom_t* frame_out_hr,
const int32_t history_in[HISTORY_SIZE],
const exponent_t history_in_exp,
const headroom_t history_in_hr)
{
// log2() of max product of two signed 32-bit integers
// = log2(INT32_MIN * INT32_MIN) = log2(-2^31 * -2^31)
const exponent_t INT32_SQUARE_MAX_LOG2 = 62;
// log2(TAP_COUNT)
const exponent_t TAP_COUNT_LOG2 = 10;
// 首先,确定要使用的输出指数。
const headroom_t total_hr = history_in_hr + filter_bfp.hr;
const exponent_t acc_exp = history_in_exp + filter_bfp.exp;
const exponent_t result_scale = INT32_SQUARE_MAX_LOG2 - total_hr
+ TAP_COUNT_LOG2;
const exponent_t desired_scale = 31;
const right_shift_t acc_shr = result_scale - desired_scale;
*frame_out_exp = acc_exp + acc_shr;
// 现在,使用该指数计算输出样本值。
for(int s = 0; s < FRAME_SIZE; s++){
timer_start(TIMING_SAMPLE);
frame_out[s] = filter_sample(&history_in[FRAME_SIZE-s-1],
acc_shr);
timer_stop(TIMING_SAMPLE);
}
// 最后,计算输出帧的头空间。
*frame_out_hr = calc_headroom(frame_out, FRAME_SIZE);
}
filter_frame()在这个阶段必须做几件事情。frame_out[]是输出尾数向量,frame_out_exp是与输出样本向量关联的指数,frame_out_hr是与输出样本向量关联的头空间,它们都是该函数的输出。
在filter_frame()能够计算输出样本之前,它必须确定在调用filter_sample()时应用的累加器的适当位移。
这是可以由带符号32位乘法得到的最大幅度整数。请注意:
TAP_COUNT_LOG2是。这个值是我们需要的额外结果位数,因为我们要将1024个项相加。
如果我们有两个int32_t的向量,其中每个元素都是INT32_MIN,那么内积将是:
无论哪个输入向量的头空间中的每一位都会减少最坏情况下的一个位数,这就是为什么它们从INT32_SQUARE_MAX_LOG2和TAP_COUNT_LOG2的和中减去以获得result_scale。
然而,为了将结果存储在int32_t中,值不能大于,这是我们的desired_scale。这两个值(或者说它们的log2())之间的差异就是所需的位移量acc_shr。
有了acc_shr,就可以计算输出样本。
输出指数只是累加器指数acc_exp加上累加器位移acc_shr。这是因为右移acc_shr位等效于(最多有舍入误差)除以。递增累加器指数可以防止移位实际上改变所表示的值。
请注意,如果将acc_shr定义为左移而不是右移,我们将减去acc_shr。
最后,需要计算输出样本的头空间,我们通过调用calc_headroom()来实现。
// 对一个输出样本应用滤波器,产生一个输出样本
int32_t filter_sample(
const int32_t sample_history[TAP_COUNT],
const right_shift_t acc_shr)
{
// 用于累加部分结果的累加器。
int64_t acc = 0;
// 计算历史记录和系数向量之间的内积。
for(int k = 0; k < TAP_COUNT; k++){
int64_t b = sample_history[k];
int64_t c = filter_bfp.data[k];
acc += (b * c);
}
return sat32(ashr64(acc, acc_shr));
}
filter_sample()看起来与第2部分中的代码相似。它的功能与第2部分中的代码相同,只是它以函数参数的形式接收acc_shr,而不是根据输入、滤波器和输出指数来计算它。
// 发送一个新的音频数据帧
static inline
void tx_frame(
const chanend_t c_audio,
const int32_t frame_out[],
const exponent_t frame_out_exp,
const headroom_t frame_out_hr,
const unsigned frame_size)
{
const exponent_t output_exp = -31;
// 输出通道期望使用指数为output_exp的PCM样本,因此在发送之前,我们需要将样本转换为使用正确指数的样本。
const right_shift_t samp_shr = output_exp - frame_out_exp;
timer_stop(TIMING_FRAME);
for(int k = 0; k < frame_size; k++){
int32_t sample = frame_out[k];
sample = ashr32(sample, samp_shr);
chan_out_word(c_audio, sample);
}
}
tx_frame()看起来与之前的阶段非常相似,只是这一次,当调用tx_frame()时,与frame_out[]关联的指数可能不是所需的输出指数-31。因此,tx_frame()必须在通过通道发送之前将样本尾数强制转换为Q1.31定点格式。
使用右移将某个指数转换为所需指数的转换始终是相同的:
// 接收一个新的音频数据帧
static inline
void rx_frame(
int32_t frame_in[FRAME_SIZE],
exponent_t* frame_in_exp,
headroom_t* frame_in_hr,
const chanend_t c_audio)
{
// 我们事先知道传入的样本将具有固定指数input_exp,并且没有理由更改它,因此我们将使用它。
*frame_in_exp = -31;
for(int k = 0; k < FRAME_SIZE; k++)
frame_in[k] = chan_in_word(c_audio);
timer_start(TIMING_FRAME);
// 确保头空间是正确的
*frame_in_hr = calc_headroom(frame_in, FRAME_SIZE);
}
rx_frame()也与之前的阶段相似,只是这一次它必须将与输入帧相关联的指数报告给调用函数。
在我们的情况下,我们正在接收PCM样本的帧,所有样本都使用指数-31。
// 接收一个新的音频数据帧并将其合并到样本历史记录中
static inline
void rx_and_merge_frame(
int32_t sample_history[HISTORY_SIZE],
exponent_t* sample_history_exp,
headroom_t* sample_history_hr,
const chanend_t c_audio)
{
// 用于存放接收到的输入帧的BFP向量。
struct {
int32_t data[FRAME_SIZE]; // 样本数据
exponent_t exp; // 指数
headroom_t hr; // 头空间
} frame_in = {{0},0,0};
// 接收一个新的输入帧
rx_frame(frame_in.data,
&frame_in.exp,
&frame_in.hr,
c_audio);
// 如果需要,重新调整BFP向量的尺度以便进行合并
const exponent_t min_frame_in_exp = frame_in.exp - frame_in.hr;
const exponent_t min_history_exp = *sample_history_exp - *sample_history_hr;
const exponent_t new_exp = MAX(min_frame_in_exp, min_history_exp);
const right_shift_t hist_shr = new_exp - *sample_history_exp;
const right_shift_t frame_in_shr = new_exp - frame_in.exp;
if(hist_shr) {
for(int k = FRAME_SIZE; k < HISTORY_SIZE; k++)
sample_history[k] = ashr32(sample_history[k], hist_shr);
*sample_history_exp = new_exp;
}
if(frame_in_shr){
for(int k = 0; k < FRAME_SIZE; k++)
frame_in.data[k] = ashr32(frame_in.data[k], frame_in_shr);
}
// 现在我们可以合并新的帧(反向顺序)
for(int k = 0; k < FRAME_SIZE; k++)
sample_history[FRAME_SIZE-k-1] = frame_in.data[k];
// 确保头空间是正确的
*sample_history_hr = calc_headroom(sample_history, HISTORY_SIZE);
}
rx_and_merge_frame()做了一些复杂的事情。首先,它声明了一个新的临时BFP向量来保存接收到的输入帧,并调用rx_frame()来获取 一个。
我们事先知道(因为我们刚刚看到了rx_frame())与输入帧关联的指数将是-31,因此我们可以只是覆盖sample_history[]中的值,就像我们在之前的阶段中所做的那样,并完成。然而,在这里我们不假设我们知道这一点。
如果与sample_history[]关联的指数sample_history_exp和接收到的帧frame_in.exp的指数不相同,我们就有了一个问题。如果我们只是简单地覆盖sample_history[]中的值,我们将会得到样本历史的不同部分使用不同的指数,这是不允许的。
为了解决这个问题,我们需要将这两个向量都强制转换为使用相同的指数。
当_添加_或_减去_两个BFP向量的尾数时,这总是发生的情况。将具有不同指数的尾数相加或相减在算术上是没有意义的。为了理解原因,让我们看一个简单的例子。
我们有两个16位尾数,和,与指数和关联。回想一下,指数的负数也是小数部分的位数。因此,我们可以用适当的基数将和重写为二进制,并尝试将它们相加:
x = 0 0.1 0 0 1 1 1 1 1 1 1 1 0 1 0
- y = 0 1 0 1.0 1 1 0 1 1 1 0 0 0 1 1
0 1 1 1 1 1 1 0 1 1 0 1 1 1 0 1 ??
但是我们应该在哪里放小数点?要将它们相加,它们必须在小数点上对齐,这意味着使用相同的指数。
决定使用哪个新的指数是确定在向量中的最小指数之间的最大最小指数的问题。让我解释一下。
如果BFP向量的指数是且头空间为,那么我们也可以使用新的指数来表示。为了使用这个新的指数,尾数需要左移位。然后,新的指数是,头空间。在这种情况下,是可以用来表示向量的最小指数,因为任何较低的指数都会导致尾数溢出。
如果向量具有最小指数,向量具有最小指数,并且我们计划用一个向量的尾数替换另一个向量的尾数,那么我们可以将和的最大值作为它们都可以强制转换为的新指数,而不会发生溢出。
虽然我们可以保证可以 在不产生任何量化误差的情况下用于,但将应用于(或反之亦然)可能会导致量化误差。
这并不完美。在我们的情况下,我们正在用frame_in.data[]的值替换sample_history[]的前FRAME_SIZE个样本。如果sample_history[]的最大值在前FRAME_SIZE个样本中,sample_history_hr可能会低估可用的头空间。但是,由于头空间的工作方式,对头空间的低估总是安全的,尽管可能会产生比严格必要的量化误差更多的量化误差。
一旦为样本历史记录选择了新的指数,就会为每个向量计算移位值(通常的方式:desired_exponent - current_exponent),并移位向量的尾数。一旦新的样本合并到sample_history[]中,就会计算sample_history[]的新的头空间。
结果
时间
| 时间类型 | 测量时间 |
|---|---|
| 每个滤波器系数 | 58.82 ns |
| 每个输出样本 | 60235.73 ns |
| 每个帧 | 15604809.00 ns |
输出波形
