在微控制器上近似两个平方和的平方根

Posted

技术标签:

【中文标题】在微控制器上近似两个平方和的平方根【英文标题】:Approximating the square root of sum of two squares on a microcontroller 【发布时间】:2011-07-28 12:28:57 【问题描述】:

为了好玩,我正在努力在 8 位微控制器 (HCS08) 上实现 FFT 算法。算法完成后,我将拥有一个 8 位实数/虚数对的数组,并且我想要找到每个值的大小。也就是如果x很复杂,我想找

|x| = sqrt(Rex^2 + Imx^2)  

现在我可以使用一个 16 位寄存器和一个 8 位寄存器。我想过只是将它们平方,相加,然后取结果的平方根,但这带来了一个问题:两个 8 位数字的平方和的最大可能值为 ~130k,大于一个 16 位寄存器可以保存的最大值(65.5k)。

我想出了一个计算 16 位数字的整数平方根的子程序,这似乎工作得很好,但显然我不能保证使用适合 16 位的值。我现在的想法是有一种算法可以直接近似我需要的东西,但我似乎找不到任何东西。任何想法将不胜感激。

总结一下:假设我有一个包含两个 8 位分量的向量,我想找到向量的长度。在不实际计算平方和平方根的情况下如何近似这个?

谢谢!

【问题讨论】:

CORDIC 算法 (en.wikipedia.org/wiki/CORDIC) 可用于将向量 <x,y> 旋转到某个新向量 <x1,0>(或等价于 <0,y1>x1(或 y1)给出原始向量的大小,CORDIC可以不用乘法来实现。不过我自己没做过,不知道有多难。 这是用于音频的吗?之后你会计算 log10 以获得 dB 值吗? 取决于目的:如果您需要长度,则没有其他方法可以计算,但是当您确实需要范数(通常是长度)时,您可以使用另一个范数而不是默认的 L2 范数,例如曼哈顿距离 (= |real|+|imag|)。 @Paul R:是的,这是针对我正在处理的音频项目的。我正在与之连接的硬件需要一个线性电压,并将其转换为对数刻度。 @user599599:好的,在这种情况下,您可能可以摆脱sqrt - 请参阅下面的答案。 【参考方案1】:

如果和大于 65535,则除以 4(右移 2 位),取平方根,再乘以 2。你会损失一位精度,结果自然不能保证适合 8 位。

【讨论】:

感谢您的回复。我唯一担心的是,如果总和大于 65535,它会溢出,我无法知道。 (我只有一个 16 位寄存器,所以添加两个 16 位数字可能会产生不可预知的结果。)我想我可以通过最初将 Rex 和 Imx 除以 2,然后将最终结果相乘来完成同样的事情回答2;这听起来和你的建议一样吗? 你已经接受了这个答案,所以我猜你想通了:将输入除以 4,然后将输出乘以 2。【参考方案2】:

嗯,你可以把 x 写成极坐标形式:

x = r[cos(w) + i sin(w)]

w = arctan(Im(x)/Re(x)),所以

|x| = r = Re(x)/cos(w)

这里没有大数字,但也许你会失去三角函数的精度(也就是说,如果你可以访问三角函数:-/)

【讨论】:

嗯,有趣的想法。不幸的是,我无法访问三角函数,而且微控制器也不支持浮点,所以我几乎仅限于基本的整数运算。不过,我计划有一个三角查找表,所以我会记住这一点。【参考方案3】:

有一个网页描述了Fast Magnitude Estimator。基本思想是让最小二乘(或其他高质量)拟合方程:

Mag ~= Alpha * max(|I|, |Q|) + Beta * min(|I|, |Q|)

对于系数 Alpha 和 Beta。列出了几个系数对,包括均方误差、最大误差等,包括适用于整数 ALU 的系数。

【讨论】:

看起来 61/64 选项中的一个很适合这个应用程序。【参考方案4】:

一种可能适合也可能不适合的廉价而肮脏的方法是使用

|x| ~ max(|Rex|,|Imx|) + min(|Rex|,|Imx)/2;

这往往会高估|x| 0 到 12% 之间。

【讨论】:

【参考方案5】:

如果您随后要将幅度转换为 dB,那么您将完全放弃 sqrt 操作。 IE。如果你的计算是:

magnitude = sqrt(re*re+im*im); // calculate magnitude of complex FFT output value
magnitude_dB = 20*log10(magnitude); // convert magnitude to dB

您可以将其重写为:

magnitude_sq = re*re+im*im; // calculate squared magnitude of complex FFT output value
magnitude_dB = 10*log10(magnitude_sq);  // convert squared magnitude to dB

【讨论】:

好点,但我的问题是 log10 也是一个计算量很大的操作。我仍然有找到最接近的整数或使用查找表的问题。 @user599599:是的,你仍然有log,但之前你有sqrt + log,现在你只有log【参考方案6】:

您可能仅限于 2 个寄存器,但您可以在 http://www.realitypixels.com/turk/opensource/index.html 查看此代码 定点平方根 使用 CORDIC 的定点三角函数

【讨论】:

【参考方案7】:

一种可能的替代方法是计算 sqrt((x*x+y*y)/2,它将所有可能的矢量幅度缩放到 0..255 范围内。

两种(快速)算法似乎可以提供近乎完美的结果,一种使用 Cordic,另一种使用最多点积。

void cordic_it(uint16 &x, uint16 &y, int n) 
    auto X = x + y >> n;  // vsraq_n_u16(x, y, n)  in arm neon
    y = abs(y - x >> n);  // vabdq_u16(y, x >> n)  in arm neon


uint16_t scaled_magnitude_cordic(uint8_t x, uint8_t y) 
     const int kRound = 1;
     if (x < y) std::swap(x,y);
     // multiply by factor of 256/sqrt(2) == 181.02
     // then reduce by the gain of the cordic iterations of 1.16
     // - with prescaling we also ensure, that the cordic iterations
     //   do not lose too much significant bits when shifting right
     uint16_t X = x * 156, Y = y * 156;
     // exactly 4 iterations. 3 is too little, 5 causes too much noise
     for (int j = 1; j <= 4; j++)  cordic_it(X,Y,j);
     return (X+kRound) >> 8;

通过改变 kRound,可以调整结果:

     Histogram of real - approx:   -1    0       1 
kRound == 0 -> smaller code         1    46617   18918
kRound == 1 -> approx >= real       0    46378   19158
kRound == -73 ->  balanced error    3695 58301   3540

在选择kRound == 1时,可以通过

修复所有结果
uint8_t fix_if_larger_by_one(uint8_t sqrt, uint8_t x, uint8_t y) 
    auto P = (x*x + y*y) / 2;
    auto Q = sqrt*sqrt;
    return sqrt - (P < Q);

也可以通过对多个角度近似 xa + yb 的点积来计算平方根,其中传统方法是使用单个角度a = 1, b = 1/2

有 5 个独特的角度,对于大约 [0 10 20 30 40][5 15 25 35 45] 的角度,可以得出任意一组系数,这两个系数都会产生近乎完美的结果,最多相差 1 个单位。

1) [181 0],  [178 31], [170 62], [157 91], [139 116]
2) [180 18], [175 46], [164 76], [148 104], [128 128]

选项 1 有 9 个非平凡系数(尽管 62 == 31*2)。 选项 2 有 8 个非平凡系数,可用于以下实现:

int approx(uint8_t x, uint8_t y) 
     if (x < y) std::swap(x,y);  // sort so that x >= y
     auto a4 = (x + y) / 2;      // vhaddq_u8(x,y) on Arm Neon
     auto a0 = (x * 180 + y * 18) >> 8;
     auto a1 = (x * 175 + y * 46) >> 8;
     auto a2 = (x * 164 + y * 76) >> 8;
     auto a3 = (x * 148 + y * 104) >> 8;
     return max_of_five_elements(a0,a1,a2,a3,a4);

这组大部分是偶数的系数可以很好地转换为带有_mm_maddubs_epi16_mm_max_epu16 的SSSE3 指令集:除了a1 之外的每个点积都可以通过交错x、y 和交错系数的一条指令轻松计算。当然,同时计算 16 个相邻近似值以消除延迟并且不浪费来自 _mm_packus_epi16 的任何计算、对 uint8_t 输入进行排序或平均更有意义。

auto a0 = _mm_maddubs_epi16(xy, coeffs0); // coeffs0 = 90 9 90 9 ...
auto a1 = _mm_maddubs_epi16(xy, coeffs1); // coeffs1 = 87 23 87 23 ...
auto a2 = _mm_maddubs_epi16(xy, coeffs2); // coeffs2 = 82 38 82 38 ...
auto a3 = _mm_maddubs_epi16(xy, coeffs3); // coeffs3 = 74 52 74 52 ...
auto a4 = _mm_maddubs_epi16(xy, coeffs4); // coeffs4 = 64 64 64 64 ...
a1 = _mm_add_epi16(a1, x_per_2); // LSB of the coefficient 87.5

// take the maximum, shift right by 7 and pack to uint8_t
a0 = _mm_max_epu16(a0, a1);
a0 = _mm_max_epu16(a0, a2);
a0 = _mm_max_epu16(a0, a3);
a0 = _mm_max_epu16(a0, a4);
a0 = _mm_srli_epi16(a0, 7);
a0 = _mm_packus_epi16(a0, a0);

仅使用 8 个系数也适用于 ARM Neon 实现,它现在可以使用 16 位乘 16 位标量乘法,将所有系数存储在一个全宽寄存器中。

为了获得完美的结果,点积算法必须补偿到另一个方向,因为它可能会给出值,这只是floor(sqrt((x*x+y*y)/2) 的参考实现下面的一个元素:

uint8_t fix_if_smaller_by_one(uint8_t sqrt, uint8_t x, uint8_t y) 
    auto P = (x*x + y*y) / 2;
    auto Q = (sqrt+1)*(sqrt+1);
    return sqrt + (Q <= P);

其他近似算法通常使用除法或缩放,这在 AVX2 之前的英特尔中很难矢量化,因为缺乏可变的每通道移位。

【讨论】:

以上是关于在微控制器上近似两个平方和的平方根的主要内容,如果未能解决你的问题,请参考以下文章

基于微控制器的日出/设置算法实现

在两个 arduino 微控制器上使用 i2c 通信发送字符串

如何使用 JTAG 链接调试两个或多个微控制器?

ARM Chromebook 上的微控制器开发

微控制器有哪些型别

在微控制器或嵌入式平台上执行自定义脚本