我正在努力在 8 位微控制器 (HCS08) 上以汇编方式实现 FFT 算法，以获得乐趣。算法完成后，我将拥有一个 8 位实数/虚数对的数组，并且我想要找到每个值的大小。也就是说，如果x很复杂，我想找到
|x| = sqrt(Re{x}^2 + Im{x}^2)
现在我可以使用一个 16 位寄存器和一个 8 位寄存器。我想只是将它们平方，相加，然后取结果的平方根，但这带来了一个问题：两个 8 位数字的平方和的最大可能值为 ~130k，大于一个 16 位寄存器可以保存的最大值（65.5k）。
我想出了一个计算 16 位数字的整数平方根的子程序，这似乎工作得很好，但显然我不能保证使用适合 16 位的值。我现在的想法是有一种算法可以直接近似我需要的东西，但我似乎找不到任何东西。任何想法将不胜感激。
总结一下：假设我有一个包含两个 8 位分量的向量，我想找到向量的长度。在不实际计算平方和平方根的情况下如何近似这个？
谢谢！
有一个描述Fast Magnitude Estimator的网页。基本思想是让最小二乘（或其他高质量）拟合方程：
Mag ~= Alpha * max(|I|, |Q|) + Beta * min(|I|, |Q|)
对于系数 Alpha 和 Beta。列出了几个系数对，包括均方误差、最大误差等，包括适用于整数 ALU 的系数。
如果和大于 65535，则除以 4（右移 2 位），取平方根，再乘以 2。你会损失一位精度，结果自然不能保证适合 8位。
一种可能的替代方法是计算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);
}
也可以通过对多个角度近似 x a + y b的点积来计算平方根，其中传统方法是使用单个角度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_epu16instrinsics 的 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 之前的英特尔中很难矢量化，因为缺乏可变的每通道移位。
如果您随后要将幅度转换为 dB，那么您将sqrt完全放弃该操作。即，如果您的计算是：
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
好吧，你可以把 x 写成极坐标形式：
x = r[cos(w) + i sin(w)]
哪里w = arctan(Im(x)/Re(x))，所以
|x| = r = Re(x)/cos(w)
这里没有大数字，但也许你会失去三角函数的精度（也就是说，如果你可以访问三角函数：-/）
一种可能适合也可能不适合的廉价而肮脏的方法是使用
|x| ~ max(|Re{x}|,|Im{x}|) + min(|Re{x}|,|Im{x})/2;
这将倾向于高估 |x| 介于 0 到 12% 之间。
您可能仅限于 2 个寄存器，但您可以在http://www.realitypixels.com/turk/opensource/index.html 使用 CORDIC 的定点平方根定点三角函数查看此代码