查表法实现FFT

只看该作者 · 2025-2-5 22:59

选择合适的 FFT 点数对于查表法的实现非常重要。不同的 FFT 点数对应着不同大小的查找表，点数越大，查找表的规模也越大，所需的内存空间也越多。同时，FFT 点数还会影响计算的复杂度和精度，需要根据实际应用需求进行权衡。

只看该作者 · 2025-2-6 11:05

查表法需要预先计算并存储一部分或全部可能的DFT结果。因此，查表的大小会直接影响内存的使用。对于64点的FFT，查表的大小相对较小，但对于更大的点数（如1024点），查表的存储需求会显著增加。
需要确保目标设备有足够的内存来存储查表数据。

只看该作者 · 2025-2-6 12:21

本帖最后由 andy520520 于 2025-2-6 12:32 编辑

Y1 = X1 + X2 * twiddle factor

Y2 = X1 - X2 * twiddle factor

X2 * twiddle factor 是不是反复乘法了？这样是不是没有减少乘法次数？

奇偶项，分成N/2 两部分

Y(k)          = even(k)  + W(N->K) * odd(k)

Y(k + N/2) = even(k)  -  W(N->K) * odd(k)
基于递归调用的FFT, deepseek 给出的

// 快速傅里叶变换
void fft(complex double* x, int N) {
if (N <= 1) return;

// 分割偶数点和奇数点
complex double even[N/2];
complex double odd[N/2];
for (int i = 0; i < N/2; i++) {
      even = x[2*i];
      odd = x[2*i + 1];
}
// 递归计算偶数点和奇数点的FFT
fft(even, N/2);
fft(odd, N/2);

// 合并结果
for (int k = 0; k < N/2; k++) {
      complex double t = cexp(-2.0 * I * PI * k / N) * odd[k];
      x[k] = even[k] + t;
      x[k + N/2] = even[k] - t;
} }

cexp(-2.0 * I * PI * k / N) 也可以改成查表法的

只看该作者 · 2025-2-6 15:11

考虑查找表本身可能存在的误差，以及查表过程中可能引入的误差。
根据应用需求确定可接受的误差范围，并在设计和实现过程中采取相应的措施来减小误差。

只看该作者 · 2025-2-6 15:34

使用标准FFT算法的结果作为参考，验证查表法实现的FFT的正确性。

只看该作者 · 2025-2-6 15:56

在FFT计算过程中，需要根据当前的蝶形运算需要，从查找表中检索相应的乘法结果进行组合。要确保正确地选择和读取查找表中的值，避免因索引错误或数据类型不匹配而导致的计算错误。

只看该作者 · 2025-2-6 17:21

查表操作应尽可能高效，避免在计算过程中引入过多的延迟。可以通过优化查表算法和使用快速查找技术来提高效率

2万 · 2025-2-6 17:41

查找表的大小与FFT的点数和数据类型有关。对于较大的FFT点数，查找表的存储需求也会相应增加。

1万 · 2025-2-6 18:00

查表法只能适用部分场景吧

只看该作者 · 2025-2-6 18:01

查表**增加内存的使用量，因此需要合理分配内存，避免内存不足的问题。

只看该作者 · 2025-2-6 18:21

查找表通常存储在内存中，频繁的内存访问可能会成为计算的瓶颈。为了提高内存访问效率，可以采用一些优化策略，如将查找表存储在高速缓存中，或者合理安排数据的存储顺序，减少内存访问冲突。

只看该作者 · 2025-2-6 18:40

对代码进行优化，减少循环开销，提高计算效率。

只看该作者 · 2025-2-6 18:45

andy520520 发表于 2025-2-6 12:21
Y1 = X1 + X2 * twiddle factor

Y2 = X1 - X2 * twiddle factor

void NuFFT(int *pi32RealBuffer, int *pi32ImageBuffer, uint32_t u32Len) {
const uint32_t N = u32Len;
const uint32_t u32TotalStage = log(N) / log(2);

uint32_t u32Stage, u32Span, u32Node, u32Twiddle;
int32_t i32X1Real, i32X1Image, i32X2Real, i32X2Image;
double cosVal, sinVal;

// 进行位反转
bit_reverse(pi32RealBuffer, pi32ImageBuffer, N);

// Iterations for log2(N) FFT butterfly stages
for (u32Stage = 0; u32Stage < u32TotalStage; u32Stage++) {
      // 使用位运算替代 pow 函数
      u32Span = 1 << u32Stage;

      for (u32Twiddle = 0; u32Twiddle < u32Span; u32Twiddle++) {
         // 计算旋转因子的角度
         double angle = -2 * M_PI * u32Twiddle / (2 * u32Span);
         cosVal = cos(angle);
         sinVal = sin(angle);

         for (u32Node = u32Twiddle; u32Node < N; u32Node += 2 * u32Span)
         {
            // Get the real and image part of the input variable X1
            i32X1Real = pi32RealBuffer[u32Node];
            i32X1Image = pi32ImageBuffer[u32Node];

            // Get the real and image part of the input variable X2
            i32X2Real = pi32RealBuffer[u32Node + u32Span];
            i32X2Image = pi32ImageBuffer[u32Node + u32Span];

            // 计算 X2 乘以旋转因子的结果
            int32_t twiddleReal = RESCALE(i32X2Real * cosVal - i32X2Image * sinVal);
            int32_t twiddleImage = RESCALE(i32X2Real * sinVal + i32X2Image * cosVal);

            // Y1 = X1 + X2 * twiddle factor
            pi32RealBuffer[u32Node] = i32X1Real + twiddleReal;
            pi32ImageBuffer[u32Node] = i32X1Image + twiddleImage;

            // Y2 = X1 - X2 * twiddle factor
            pi32RealBuffer[u32Node + u32Span] = i32X1Real - twiddleReal;
            pi32ImageBuffer[u32Node + u32Span] = i32X1Image - twiddleImage;
         }
      }
}
}

X2 * twiddle factor 计算一次这个就可以了
用上面的函数实现了FFT , 我做了修改，把缩放改成了1倍，直接用cos,sin 代替查表值

只看该作者 · 2025-2-6 19:00

优化查表算法，提高查找速度。例如，可以通过排序、索引或哈希表等技术来加速查找过程。
避免在查表过程中进行不必要的计算或重复查找。

只看该作者 · 2025-2-6 23:36

在实际计算过程中，需要对查找表中的数据进行高效的索引和查找。查找表的索引方法应该尽可能简单有效，以减少查找时间。

只看该作者 · 2025-2-7 12:12

本帖最后由 andy520520 于 2025-2-7 12:19 编辑

zhuotuzi 发表于 2025-1-24 17:58

#if DEBUG

            /* Y1 = X1 + X2 * twiddle factor , 注意这里是  -angle , 为了与原文一致 double angle = 2 * M_PI * u32Twiddle / (2 * u32Span) 用正的角度
               * The real part is real * cos() + image * sin()
               * The image part is -real * sin() + image * cos()
               */

            pi32RealBuffer[u32Node]  = i32X1Real  + RESCALE(i32X2Real * cos(angle)) + RESCALE(i32X2Image * sin(angle));
            pi32ImageBuffer[u32Node] = i32X1Image - RESCALE(i32X2Real * sin(angle)) + RESCALE(i32X2Image * cos(angle));

            /* Y2 = X1 - X2 * twiddle factor
               * The real part is  -real * cos() - image * sin()
               * The image part is  real * sin() - image * cos()
               */
            pi32RealBuffer[u32Node + u32Span]  = i32X1Real  - RESCALE(i32X2Real * cos(angle)) - RESCALE(i32X2Image * sin(angle));
            pi32ImageBuffer[u32Node + u32Span] = i32X1Image + RESCALE(i32X2Real * sin(angle)) - RESCALE(i32X2Image * cos(angle));

            #else

            uint32_t x2_tw_real  = RESCALE(i32X2Real * cos(angle))  - RESCALE(i32X2Image * sin(angle));
            uint32_t x2_tw_imag  = RESCALE(i32X2Real * sin(angle))  + RESCALE(i32X2Image * cos(angle));

            pi32RealBuffer [u32Node] = i32X1Real  + x2_tw_real ;
            pi32ImageBuffer[u32Node] = i32X1Image + x2_tw_imag  ;

            pi32RealBuffer [u32Node + u32Span] = i32X1Real  - x2_tw_real  ;
            pi32ImageBuffer[u32Node + u32Span] = i32X1Image - x2_tw_imag  ;

            #endif

我发现存在的问题点  ： ① 一个是偶数项乘以旋转因子是重复乘的，这个就不叫FFT了； ②、输入的数据没有倒序，这个造成输出完全错误 ; ③、实部和虚部都是整数，实用性几乎没有。上面验证的直接用cos ,sin 计算的

只看该作者 · 2025-2-7 12:49

andy520520 发表于 2025-2-7 12:12
#if DEBUG

/* Y1 = X1 + X2 * twiddle factor , 注意这里是 -angle , 为了与原文一致 ...

为了对照结果，我把整数改成了虚部和实部都改成了 float

void NuFFT(float *pi32RealBuffer, float *pi32ImageBuffer, uint32_t u32Len)
{
const uint32_t N = u32Len;
const uint32_t u32TotalStage = log(N) / log(2);

uint32_t u32Stage, u32Span, u32Node, u32Twiddle;

float i32X1Real, i32X1Image, i32X2Real, i32X2Image;

// 添加位反转重排调用

bitReverseReorder(pi32RealBuffer, pi32ImageBuffer, N);

/* Iterations for log2(N) FFT butterfly stages */

for (u32Stage = 0; u32Stage < u32TotalStage; u32Stage++)
{
      // Span indicates the buffer index for constituting a butterfly matrix

      u32Span = 1 << u32Stage;

      for (u32Twiddle = 0; u32Twiddle < u32Span; u32Twiddle++) {

         // 计算旋转因子的角度

   #if DEBUG

         double angle = 2 * M_PI * u32Twiddle / (2 * u32Span);
         #else

   double angle = -2 * M_PI * u32Twiddle / (2 * u32Span);

   #endif

         for (u32Node = u32Twiddle; u32Node < N; u32Node += 2 * u32Span) {
            /* Get the real and image part of the input variable X1 */
            i32X1Real = pi32RealBuffer[u32Node];
            i32X1Image = pi32ImageBuffer[u32Node];

            /* Get the real and image part of the input variable X2 */
            i32X2Real = pi32RealBuffer[u32Node + u32Span];
            i32X2Image = pi32ImageBuffer[u32Node + u32Span];

#if DEBUG

            /* Y1 = X1 + X2 * twiddle factor , 注意这里是  -angle
               * The real part is real * cos() + image * sin()
               * The image part is -real * sin() + image * cos()
               */

            pi32RealBuffer[u32Node]  = i32X1Real  + RESCALE(i32X2Real * cos(angle)) + RESCALE(i32X2Image * sin(angle));
            pi32ImageBuffer[u32Node] = i32X1Image - RESCALE(i32X2Real * sin(angle)) + RESCALE(i32X2Image * cos(angle));

            /* Y2 = X1 - X2 * twiddle factor
               * The real part is  -real * cos() - image * sin()
               * The image part is  real * sin() - image * cos()
               */
            pi32RealBuffer[u32Node + u32Span]  = i32X1Real  - RESCALE(i32X2Real * cos(angle)) - RESCALE(i32X2Image * sin(angle));
            pi32ImageBuffer[u32Node + u32Span] = i32X1Image + RESCALE(i32X2Real * sin(angle)) - RESCALE(i32X2Image * cos(angle));

#else

uint32_t x2_tw_real  = RESCALE(i32X2Real * cos(angle))  - RESCALE(i32X2Image * sin(angle));
uint32_t x2_tw_imag  = RESCALE(i32X2Real * sin(angle))  + RESCALE(i32X2Image * cos(angle));

            pi32RealBuffer [u32Node] = i32X1Real  + x2_tw_real ;
            pi32ImageBuffer[u32Node] = i32X1Image + x2_tw_imag  ;

            pi32RealBuffer [u32Node + u32Span] = i32X1Real  - x2_tw_real  ;
            pi32ImageBuffer[u32Node + u32Span] = i32X1Image - x2_tw_imag  ;

#endif
         }
      }
}
}

必须要倒序输入

只看该作者 · 2025-2-8 11:37

对于较大的fft点，查找表**会消耗大量内存。为了避免内存溢出，需要在内存使用和计算效率之间进行权衡

只看该作者 · 2025-2-8 20:45

实施有效的位反转算法来重新排列输入数据是fft算法中的关键步骤。

查表法实现FFT

核心会员奖章

坚毅之洋流

七世轮回

技术奇才奖章