module FFT (
  input wire clk,
  input wire rst,
  input wire signed [15:0] x_real [0:15],  // 输入实部
  input wire signed [15:0] x_imag [0:15],  // 输入虚部
  output reg signed [15:0] y_real [0:15],  // 输出实部
  output reg signed [15:0] y_imag [0:15]   // 输出虚部
);

  parameter N = 16;  // FFT大小，必须是2的幂次方

  // FFT模块内部信号
  reg signed [15:0] twiddle_real [0:N/2-1];
  reg signed [15:0] twiddle_imag [0:N/2-1];
  reg signed [15:0] stage_real [0:N-1];
  reg signed [15:0] stage_imag [0:N-1];
  reg signed [15:0] butterfly_real [0:N-1];
  reg signed [15:0] butterfly_imag [0:N-1];

  // 初始化蝶形运算所需的旋转因子
  initial begin
    for (int i = 0; i < N/2; i = i + 1) begin
      twiddle_real[i] = cos(2 * $pi * i / N);
      twiddle_imag[i] = -sin(2 * $pi * i / N);
    end
  end

  // FFT计算
  always @(posedge clk or posedge rst) begin
    if (rst) begin
      // 复位时初始化状态
      for (int i = 0; i < N; i = i + 1) begin
        stage_real[i] <= 16'h0000;
        stage_imag[i] <= 16'h0000;
        butterfly_real[i] <= 16'h0000;
        butterfly_imag[i] <= 16'h0000;
        y_real[i] <= 16'h0000;
        y_imag[i] <= 16'h0000;
      end
    end else begin
      // 输入数据赋值到第一级蝶形运算
      stage_real[0] <= x_real[0];
      stage_imag[0] <= x_imag[0];

      // FFT计算的蝶形运算
      for (int stage = 1; stage < $clog2(N); stage = stage + 1) begin
        // 计算每个蝶形运算的输入
        for (int i = 0; i < N; i = i + 1) begin
          butterfly_real[i] = stage_real[i] + (stage_imag[i] * twiddle_real[i >> stage] - stage_imag[i] * twiddle_imag[i >> stage]);
          butterfly_imag[i] = stage_imag[i] + (stage_imag[i] * twiddle_real[i >> stage] + stage_real[i] * twiddle_imag[i >> stage]);
        end

        // 更新下一级蝶形运算的输入
        for (int i = 0; i < N; i = i + 1) begin
          stage_real[i] = butterfly_real[i];
          stage_imag[i] = butterfly_imag[i];
        end
      end

      // 输出结果
      for (int i = 0; i < N; i = i + 1) begin
        y_real[i] = butterfly_real[i];
        y_imag[i] = butterfly_imag[i];
      end
    end
  end

endmodule