[Feature](mluOpExecFFT): add bluestein fft

docs/design_docs/fft/fft.md

@@ -105,6 +105,83 @@ $$
 
 可以将二维FFT看作是在输入的两个维度分别做一维FFT，因此其计算方式与一维FFT本质相同，这里不再对FFT计算细节进行赘述。
 
+#### 1.2.1 bluestein 算法简述
+算法是一种用于计算任意长度离散傅里叶变换（DFT）的快速算法，通过将DFT转换为循环卷积实现计算优化。以下是其核心步骤及公式推导
+
+##### 1.2.1.1 DFT二次项分解
+标准DFT表达式：
+
+$$
+X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j\frac{2\pi}{N}kn}
+$$
+
+通过引入二次项展开：
+
+$$
+kn = \frac{1}{2}[k^2 + n^2 - (k-n)^2]
+$$
+
+得到等效表达式：
+
+$$
+e^{-j\frac{2\pi}{N}kn} = e^{-j\frac{\pi}{N}k^2} \cdot e^{j\frac{\pi}{N}(k-n)^2} \cdot e^{-j\frac{\pi}{N}n^2}
+$$
+
+##### 1.2.1.2 卷积形式重构
+令：
+
+$$
+\begin{aligned}
+w[n] &= e^{-j\frac{\pi}{N}n^2} (chrizp 信号) \\
+h[n] &= e^{j\frac{\pi}{N}n^2} (辅助信号)
+\end{aligned}
+$$
+
+则DFT可表示为：
+
+$$
+X[k] = w[k] \cdot \sum_{n=0}^{N-1} (x[n]w[n]) \cdot h[k-n]
+$$
+
+##### 1.2.1.3 计算步骤
+
+###### 1.2.1.3.1 预处理
+1. **输入信号调制**：
+
+$$
+a[n] = x[n] \cdot w[n], \quad 0 \leq n < N
+$$
+
+2. **构造卷积核**：
+
+$$
+h[n] = 
+\begin{aligned} 
+e^{j\frac{\pi}{N}n^2}, & 0 \leq n < N \\
+\end{aligned}
+$$
+
+###### 1.2.1.3.2 快速卷积
+1. **补零扩展**：
+
+$$
+\begin{aligned}
+a_{pad} &= [a[0],...,a[N-1], \underbrace{0,...,0}_{M-N}]
+\end{aligned}
+$$
+
+$$
+\begin{aligned}
+h_{pad} &= [h[0],...,hj[N-1], \underbrace{0,...,0}_{M-N}]
+\end{aligned}
+$$
+
+2. **FFT加速计算**：
+
+$$
+X[k] = w[k] \cdot \text{IFFT}\left( \text{FFT}(a_{pad}) \odot \text{FFT}(h_{pad}) \right)[k]
+$$
+
 备注：
 
 需要说明对nan/inf的特殊处理，输入存在nan/inf时，输出为按照IEEE754标准根据FFT计算公式得到，其中任何数乘nan为nan，0乘inf为nan，非0乘正负inf根据符号为正负inf，任何数加nan为nan，任何非nan数加正负inf根据符号为正负inf。
@@ -415,11 +492,83 @@ fft2d可以理解为先做一个bacth=n0, n=[n1]的行主序fft1d, 再做一个b
 
 当前测试表明，列主序这样的分解以及对应的访存模式，可以大大减小stride过大带来的影响。
 
-
 ##### 3.1.3.2 行主序batch循环软流水实现
 
 对于fft2d而言，其对行主序fft1d的调用往往会有较多batch。我们拟增加对于batch之间的软流水优化，从而提升性能。
 
+#### 3.1.2 任意长度bluestein fft实现
+#### 3.1.2.1 1d bluestein fft
+根据算法计算步骤 1.2.1.3 可知， 核心逻辑为对输入数据乘系数后得到a[n],进行pad后再进行fft，结果与另一个系数的fft 结果相乘再逆fft，结果再乘系数。
+其中关键步骤fft 及ifft 可直接调用前述fft 及ifft kernel 进行计算。需要实现的是1. x[n]*w[n], 复数矩阵每列乘以相应复数系数；2. fft(a_pad) * fft(h_pad), fft(h_pad)结果为1位复数向量, 其实质为也是一个复数矩每列阵乘以相应的复数系数;3. w[k]*ifft()也是计算复数矩阵每列乘以复数相应系数。这个三个计算步骤实质是相同的，所以实现完整的bluestein fft 算法除调用上述接口外还需要实现此功能, 命名接complex_coeff_matmul(),以及系数生成generate()
+方案如下：
+![bluesteinfft](bluestein_fft_01.jpg)
+
+- generate() 伪代码：
+```C++
+// |---------------nram size-----|
+// |nram_index|nram_n_2|nram_temp|
+// nram_size 均分为3等份 结果存在sram_buffer,供后续流程使用
+__mluop_get_stage_indices_tfuse(*nram_index, length);
+__bang_write(nram_n_2, length, value);
+__bang_square(nram_n_2, nram_n_2, length);
+__bang_mul(nram_index, nram_index, nram_n_2);
+if(chirpz){
+  __bang_mul_scalar(nram_index, nram_index, -Pi, length);
+} else {
+  __bang_mul_scalar(nram_index, nram_index, Pi, length);
+}
+__mluop_exp(nram_index, nram_temp, length);
+__bang_transpose(nram_index, nram_index, 2, length)
+memcpy(sram_buffer, nram_index, length, NRAM2SRAM);
+```
+
+- complex_coeff_matmul() 核心计算步骤伪代码
+  - 对x[b, length] 拆分b，当length 较长时，核内循环处理，核心计算步骤做3及流水 
+  - 设z1=a+bi，z2=c+di(a、b、c、d∈R)是任意两个复数，那么它们的积(a+bi)(c+di)=(ac-bd)+(bc+ad)i。
+```C++
+// input_gdram 为输入x_n的地址，output_gdram 为输出地址
+// 流水开始前对ouput空间刷零完成pad zero
+// pad_length 是bluestein fft 计算length 长度输入需要pad 到的长度 
+// pad_length >= 2* length - 1
+gdram_set(output_gdram, b*pad_length, 0);
+//
+// 3级流水内计算步骤
+// nram 内存划分：
+// |-----ping------------|------------pong-----|------|
+// |input_x|output_x|temp|input_x|output_x|temp|chirpZ|
+// 指定一次可以处理长度为 l = 8192;
+// input_x = ouput_x = temp = b*l
+// (6b+1)*l = nram_size;
+// b = (nram_size/l - 1)/6
+memcpy_async(chirpz, sram_buffer, SRAM2NRAM, 2*l);
+// z1_r, z2_r..., z1_i, z2_i....
+memcpy_async(input_x, input_gdram, GDRAM2NRAM, b*l);
+//transpoze
+// a_r b_i c_r d_i
+// e_r f_i g_r e_i
+//复数实数虚数分开
+//  转换成 a_r,c_r,e_r,g_r 
+//        b_i,d_i,f_i,e_i
+__bang_transpose(output_x, input_x, 2, b*l);
+// a_r*z1_r c_r*z2_r...e_r*z1_r,g_r*z2_r...
+// b_i*z1_i d_i*z2_i...
+__bang_cycle_mul(input_x, output_x, chirpz, b*l, l);
+__bang_cycle_mul(input_x + b*l, output_x + b*l, chirpz + l, b*l, l);
+// a_r*z1_r - b_i*z1_i ....
+__bang_sub(temp, input_x, input_x + b*l, b*l);
+
+// a_r*z1_i c_r*z2_i...e_r*z1_i,g_r*z2_i...
+// b_i*z1_r d_i*z2_r...
+__bang_cycle_mul(input_x, output_x + b*l, chirpz, b*l, l);
+__bang_cycle_mul(input_x + b*l, output_x, chirpz + l, b*l, l);
+__bang_add(temp+b*l, input_x, input_x + b*l, b*l);
+//tranpose转换回复数
+__bang_transpose(output_x, temp, b*l, 2);
+memcpy_async(output_gdram, output_x, NRAM2GDRAM, pad_length*sizeof(complex), l*sizeof(complex), b*sizeof(complex));
+```
+#### 3.1.2.1 2d bluestein fft
+用1d bluestein fft 方法做行主序 fft， 然后再做列主序fft
+
 ### 3.2调用方案设计
 
 目前实现的方案如下：