FFT

概述

FFT 实现了应用于 2D 图像的傅里叶变换，支持实值和复值输入。它返回频域中的图像表示。它对于基于信号频率的内容分析非常有用，并且能够进行在频域中更有效执行的图像操作，以及其他用途。

空间域中的输入	幅度谱

实现

FFT - 快速傅里叶变换是一种分而治之的算法，用于有效地计算复数或实值数据集的离散傅里叶变换，时间复杂度为 \(O(n \log n)\)。它是计算物理和通用信号处理中最重要和应用最广泛的数值算法之一。

离散傅里叶变换 (DFT) 通常将空间域上的复值向量 \(x_k\) 映射到其频域表示，如下所示

\[ I'[u,v] = \sum^{M-1}_{m=0} \sum^{N-1}_{n=0} I[m,n] e^{-2\pi i (\frac{um}{M}+\frac{vn}{N})} \]

其中

• \(I\) 是空间域中的输入图像
• \(I'\) 是其频域表示
• \(M\times N\) 是输入的维度

请注意，\(I'\) 保持未归一化，因为频谱最终会在管道的后续阶段（如果存在）中进行归一化。

根据图像维度，采用不同的技术以获得最佳性能

• CPU 后端
  • 当 \(M\) 或 \(N\) 可以分解为 \(2^a \times 3^b \times 5^c\) 时，使用快速路径
• CUDA 后端
  • 当 \(M\) 或 \(N\) 可以分解为 \(2^a \times 3^b \times 5^c \times 7^d\) 时，使用快速路径

一般来说，素因子越小，性能越好，即，2 的幂是最快的。

FFT 支持以下变换类型

• 复数到复数，C2C
• 实数到复数，R2C。

数据布局

数据布局严格取决于变换类型。在一般 C2C 变换的情况下，输入和输出数据都应为 VPI_IMAGE_FORMAT_2F32 类型，并且具有相同的大小。在 R2C 模式下，每个类型为 VPI_IMAGE_FORMAT_F32 的输入图像行 \((x_1,x_2,\dots,x_N)\) 会生成一个类型为 VPI_IMAGE_FORMAT_2F32 的非冗余复数元素输出行 \((X_1,X_2,\dots,X_{\lfloor\frac{N}{2}\rfloor+1})\)。在这两种情况下，输入和输出图像的高度相同。

FFT 类型	输入		输出
	类型	大小	类型	大小
C2C	VPI_IMAGE_FORMAT_2F32	\(W \times H\)	VPI_IMAGE_FORMAT_2F32	\(W \times H\)
R2C	VPI_IMAGE_FORMAT_F32	\(W \times H\)	VPI_IMAGE_FORMAT_2F32	\(\big( \lfloor\frac{W}{2}\rfloor+1 \big) \times H\)

对于实数输入，输出仅包含非冗余元素，因此仅返回频谱的左半部分。为了检索完整的厄米特（共轭对称）输出，必须完成右半部分，以便满足以下条件

\begin{align*} F_{\textbf{re}}(u,v) &= F_{\textbf{re}}(-u,-v) \\ F_{\textbf{re}}(-u,v) &= F_{\textbf{re}}(u,-v) \\ F_{\textbf{im}}(u,v) &= -F_{\textbf{im}}(-u,-v) \\ F_{\textbf{im}}(-u,v) &= -F_{\textbf{im}}(u,-v) \end{align*}

下图显示了如何完成此操作

完成输出频谱

频谱输出中的第零个元素表示 DC（0 Hz）分量。通常，为了可视化，最好将 DC 放在中心。实现此目的的一种方法是像下图所示那样移动 4 个象限

\(Q_0\) \(Q_1\) \(Q_2\) \(Q_3\)

\(\rightarrow\)

\(Q_3\) \(Q_2\) \(Q_1\) \(Q_0\)

C API 函数

有关实现该算法的限制、约束和后端的列表，请查阅以下函数的参考文档

函数	描述
vpiCreateFFT	为直接快速傅里叶变换算法创建有效载荷。
vpiSubmitFFT	在单个图像上运行直接快速傅里叶变换。

注意

此算法要求系统安装库 libcufft.so.10。

用法

语言 C/C++ Python

1. 导入 VPI 模块

   import vpi

2. （可选）将输入转换为 32 位浮点数并应用 FFT 操作。结果图像将具有复数分量，由 vpi.Format._2F32 格式表示。所有操作将由 CUDA 后端执行。

   with vpi.Backend.CUDA
   output = input.convert(vpi.Format.F32) \
   .fft()

3. （可选）将复数频率数据转换为幅度，以便显示。
   1. 将输出复制到 numpy 数组并将其转换为 np.complex64

      hfreq = output.cpu().view(dtype=np.complex64).squeeze(2)

   2. 用缺失值填充复数数据的右半部分。左半部分直接从 VPI 的输出复制。结果是一个完整的厄米特矩阵。

      if input.width%2==0
      wpad = input.width//2-1
      padmode = 'reflect'
      else:
      wpad = input.width//2
      padmode='symmetric'
      freq = np.pad(hfreq, ((0,0),(0,wpad)), mode=padmode)
      freq[:,hfreq.shape[1]:] = np.conj(freq[:,hfreq.shape[1]:])
      freq[1:,hfreq.shape[1]:] = freq[1:,hfreq.shape[1]:][::-1]

   3. 将复数频率转换为归一化的对数幅度

      lmag = np.log(1+np.absolute(freq))

   4. 移动频谱，使 DC 位于中心

      spectrum = np.fft.fftshift(lmag)

1. 初始化阶段
   1. 包含定义图像 FFT 函数和所需图像格式转换功能的头文件。

      #include <vpi/algo/ConvertImageFormat.h>
      #include <vpi/algo/FFT.h>

   2. 定义输入图像对象。

          VPIImage input = /*...*/;

   3. 创建将保存频谱的输出图像。由于输入是实数，因此仅会计算非冗余值，并且输出宽度为 \(\lfloor \frac{w}{2} \rfloor+1\)。

          int width, height;
          vpiImageGetSize(input, &width, &height);
       
          VPIImage spectrum;
          vpiImageCreate(width / 2 + 1, height, VPI_IMAGE_FORMAT_2F32, 0, &spectrum);

   4. 创建临时图像，该图像存储输入图像的浮点表示。还分配主机内存以保存从 VPI 输出计算出的完整厄米特表示。

          VPIImage inputF32;
          vpiImageCreate(width, height, VPI_IMAGE_FORMAT_F32, 0, &inputF32);
       
          float *tmpBuffer = (float *)malloc(width * 2 * height * sizeof(float));

   5. 在 CUDA 后端上创建 FFT 有效载荷。输入是实数，输出是复数。问题大小指的是输入大小。

          VPIPayload fft;
          vpiCreateFFT(VPI_BACKEND_CUDA, width, height, VPI_IMAGE_FORMAT_F32, VPI_IMAGE_FORMAT_2F32, &fft);

   6. 创建将提交算法以供执行的流。

          VPIStream stream;
          vpiStreamCreate(0, &stream);

2. 处理阶段
   1. 使用 CUDA 将输入转换为浮点型，保持输入范围不变。

          vpiSubmitConvertImageFormat(stream, VPI_BACKEND_CUDA, input, inputF32, NULL);

   2. 将 FFT 算法提交到流，传递浮点输入和输出缓冲区。由于有效载荷是在 CUDA 后端上创建的，因此它将使用 CUDA 进行处理。

          vpiSubmitFFT(stream, VPI_BACKEND_CUDA, fft, inputF32, spectrum, 0);

   3. 等待处理完成。

          vpiStreamSync(stream);

   4. 现在来看结果。首先锁定频谱数据。

          VPIImageData data;
          vpiImageLockData(spectrum, VPI_LOCK_READ, VPI_IMAGE_BUFFER_HOST_PITCH_LINEAR, &data);

   5. 将复数频率数据转换为幅度，以便显示。
      1. 用缺失值填充复数数据的右半部分。左半部分直接从 VPI 的输出复制。结果是一个完整的厄米特矩阵。

             /* 使宽度/高度为偶数 */
         width = width & -2;
         height = height & -2;
          
             /* 中心像素坐标 */
             int cx = width / 2;
             int cy = height / 2;
          
             /* 图像数据采用可主机访问的 pitch-linear 布局。 */
         assert(data.bufferType == VPI_IMAGE_BUFFER_HOST_PITCH_LINEAR);
             VPIImageBufferPitchLinear *dataPitch = &data.buffer.pitch;
          
             int i, j;
             for (i = 0; i < (int)height; ++i)
             {
                 for (j = 0; j < (int)width; ++j)
                 {
                     float re, im;
                     if (j < cx)
                     {
                         const float *pix =
         (const float *)((const char *)dataPitch->planes[0].data + i * dataPitch->planes[0].pitchBytes) +
         j * 2;
         re = pix[0];
         im = pix[1];
                     }
                     else
                     {
                         const float *pix = (const float *)((const char *)dataPitch->planes[0].data +
         ((height - i) % height) * dataPitch->planes[0].pitchBytes) +
         ((width - j) % width) * 2;
                         /* 复共轭 */
         re = pix[0];
         im = -pix[1];
                     }
          
         tmpBuffer[i * (width * 2) + j * 2] = re;
         tmpBuffer[i * (width * 2) + j * 2 + 1] = im;
                 }
             }

      2. 将复数频率转换为归一化的对数幅度

             float min = FLT_MAX, max = -FLT_MAX;
          
             for (i = 0; i < (int)height; ++i)
             {
                 for (j = 0; j < (int)width; ++j)
                 {
                     float re, im;
         re = tmpBuffer[i * (width * 2) + j * 2];
         im = tmpBuffer[i * (width * 2) + j * 2 + 1];
          
                     float mag = logf(sqrtf(re * re + im * im) + 1);
         tmpBuffer[i * width + j] = mag;
          
         min = mag < min ? mag : min;
         max = mag > max ? mag : max;
                 }
             }
          
             for (i = 0; i < (int)height; ++i)
             {
                 for (j = 0; j < (int)width; ++j)
                 {
         tmpBuffer[i * width + j] = (tmpBuffer[i * width + j] - min) / (max - min);
                 }
             }

      3. 移动频谱，使 DC 位于中心

             for (i = 0; i < (int)height; ++i)
             {
                 for (j = 0; j < (int)cx; ++j)
                 {
                     float a = tmpBuffer[i * width + j];
          
                     /* 象限 0？ */
                     if (i < cy)
                     {
                         /* 将其与象限 3 交换 */
         tmpBuffer[i * width + j] = tmpBuffer[(i + cy) * width + (j + cx)];
         tmpBuffer[(i + cy) * width + (j + cx)] = a;
                     }
                     /* 象限 2？ */
                     else
                     {
                         /* 将其与象限 1 交换 */
         tmpBuffer[i * width + j] = tmpBuffer[(i - cy) * width + (j + cx)];
         tmpBuffer[(i - cy) * width + (j + cx)] = a;
                     }
                 }
             }

   6. 由于不再需要 VPI 输出图像，只需解锁即可

          vpiImageUnlock(spectrum);

   7. 现在使用您喜欢的方法显示 tmpBuffer。
3. 清理阶段
   1. 释放流、有效载荷、输入和输出图像以及临时主机缓冲区所持有的资源。

          vpiImageDestroy(input);
          vpiImageDestroy(inputF32);
          vpiImageDestroy(spectrum);
      free(tmpBuffer);

有关更多信息，请参阅 VPI - Vision Programming Interface 的 “C API 参考” 部分中的 快速傅里叶变换。

性能

有关如何使用下表中的性能表的信息，请参阅 算法性能表。
在比较测量结果之前，请查阅 比较算法运行时间。
有关性能基准测试方式的更多信息，请参阅 性能基准。

清除过滤器 设备 Jetson AGX Orin  - 流 1 2 4 8