用SIMD指令加速求复数的模

• 在使用fftw库傅里叶求频域的时候需要求复数的模来求频谱图，但是fftw库存储复数的类型使用的complex类型，其类型是复数的实部和虚部交替存储，如果将实部和虚部结果分别存到Mat里，再通过magnitude来求，这样的做法对于大图的耗时非常高，由于算法的特性，所以考虑直接用SIMD指令求模。

• 求摸的公式是：
  
  -算法思想：
  整体算法上没有太大的难度，最难的可能就是，在complex里实部和虚部是交替存放，如何把该部分正确地取出来并作运算，本人使用的方法是这样的：
  

• 代码如下

 parallel_for_(Range(0, h), [&](const Range& range) {
                for (int i = range.start; i < range.end; i++) {
              //   for (int i = 0; i < h; i++) {
                    int limit = 0;
                    int offset = i * new_w;
                    float* pReal_line = real_ptr + offset;
                    float* pImag_line = imag_ptr + offset;
                    float* pResult_line = result + offset;
                    for (int j = 0;;) {
                       // 剩下的元素处理j
                        for (; j < limit; j++) {
                            float  a = pow(output_array[i * new_w + j][0], 2);
                            float  b = pow(output_array[i * new_w + j][1], 2);
                            *(pResult_line + j) = sqrt(a+b);
                        }
                        if (limit == new_w) {
                            break;
                        }
                      //  AVX
                        if (VTC_HAS_SUPPORT_AVX) {
                            for (; j <= new_w - 8; j += 8) {

                               //  加载实部和虚部
                                __m256 real_parts = _mm256_loadu_ps(&output_array[i * new_w + j][0]);      // 加载实部
                                __m256 imag_parts = _mm256_loadu_ps(&output_array[i * new_w + j+4][0]);  // 加载虚部

                               // 计算实部和虚部的平方
                                __m256 real_squared = _mm256_mul_ps(real_parts, real_parts);
                                __m256 imag_squared = _mm256_mul_ps(imag_parts, imag_parts);

                                __m256 hadd = _mm256_hadd_ps(real_squared, imag_squared);

                                __m256i indices = _mm256_set_epi32(7, 6, 3, 2, 5, 4, 1, 0); // 设置索引
                                __m256  rearrange = _mm256_permutevar8x32_ps(hadd, indices); // 重排vec中的元素

                                // 计算模
                                __m256 magnitudes = _mm256_sqrt_ps(rearrange);

                                // 将结果存储到输出数组
                                _mm256_storeu_ps(pResult_line+j, magnitudes);

                            }

                        }
                     //   SSE
                        for (; j <= new_w - 4; j += 4) {
                            // 加载实部和虚部到寄存器
                            __m128 real = _mm_loadu_ps(pReal_line + j);
                            __m128 imag = _mm_loadu_ps(pImag_line + j);

                            // 计算
                            __m128 real_sq = _mm_mul_ps(real, real);
                            __m128 imag_sq = _mm_mul_ps(imag, imag);

                            // 计算平方和
                            __m128 sum_sq = _mm_add_ps(real_sq, imag_sq);

                            // 计算平方和的平方根（即模）
                            __m128 result = _mm_sqrt_ps(sum_sq);
                            // 将结果存储到结果图中
                            _mm_storeu_ps(pResult_line + j, result);
                           //加载实部和虚部
                            __m128 real_parts = _mm_loadu_ps(&output_array[i * new_w + j][0]);      // 加载实部
                            __m128 imag_parts = _mm_loadu_ps(&output_array[i * new_w + j + 2][0]);  // 加载虚部

                           // 计算实部和虚部的平方
                            __m128 real_squared = _mm_mul_ps(real_parts, real_parts);
                            __m128 imag_squared = _mm_mul_ps(imag_parts, imag_parts);

                            __m128 hadd = _mm_hadd_ps(real_squared, imag_squared);

                            __m128 rearrange = _mm_shuffle_ps(hadd, hadd, _MM_SHUFFLE(3, 1, 2, 0));

                            // 计算模
                            __m128 magnitudes = _mm_sqrt_ps(rearrange);

                            // 将结果存储到输出数组
                            _mm_storeu_ps(pResult_line + j, magnitudes);
                        }
                        limit = new_w;
                    }
                }
                });