X64汇编实现FFT代码

以下是一个简单的 X64 汇编实现 FFT 的代码，用于计算长度为 2^n 的复数序列的快速傅里叶变换。

section .data
w_real dq 1.0, 0.0, 0.7071067811865475, 0.7071067811865475, 0.0, 1.0, -0.7071067811865475, 0.7071067811865475
w_imag dq 0.0, -1.0, -0.7071067811865475, 0.7071067811865475, -1.0, 0.0, -0.7071067811865475, -0.7071067811865475

section .text
global fft
fft:
    mov r12, rdi ; r12 = ptr to input
    mov r13, rsi ; r13 = ptr to output
    mov r14, rdx ; r14 = n
    mov r15, rcx ; r15 = stride

    mov rax, 1
    mov rbx, r14
    shl rax, cl
    cmp rax, r14
    jne .error

    mov rcx, 0
.loop1:
    mov rdi, r12
    mov rsi, r13
    mov rdx, r14
    mov r8, r15
    mov r9, rcx
    call fft_recursion

    add rcx, 1
    cmp rcx, rax
    jl .loop1

    ret

.error:
    xor rax, rax
    ret

fft_recursion:
    push rbx
    push rbp
    push r12
    push r13
    push r14
    push r15
    push rsi

    mov r14, rdx
    cmp r14, 2
    jle .base_case

    mov rax, 1
    shl rax, cl
    mov r15, rax

    sub rsp, 56
    mov rdi, rsp
    mov rsi, rsp+16
    mov rdx, rsp+32
    mov rcx, rsp+48
    mov rbx, r14
    shr rbx, 1
    lea r12, [rdi]
    lea r13, [rsi]
    lea r10, [rdx]
    lea r11, [rcx]

    xor r9, r9
.loop2:
    mov rax, r9
    shl rax, cl
    mov r8, rax
    mov rax, r9
    add rax, rbx
    shl rax, cl
    mov rax, r12
    add rax, r8
    mov r8, r13
    add r8, r9
    mov rax, [rax]
    mov [r8], rax
    mov rax, r9
    add rax, rbx
    mov rax, r12
    add rax, r8
    mov r8, r13
    add r8, rax
    mov rax, [rax]
    mov [r8], rax
    inc r9
    cmp r9, rbx
    jl .loop2

    mov rax, 1
    mov rbx, rbx
    shl rax, cl
    shr rbx, 1
    mov rax, r15
    shr rax, 1
    mov rdx, rax
    shl rdx, cl

    mov r9, 0
.loop3:
    mov rax, r9
    shl rax, cl
    mov r8, rax
    mov rax, r9
    add rax, rbx
    shl rax, cl
    mov r10, r12
    add r10, r8
    mov r11, r10
    add r11, rdx
    mov rax, [r10]
    mov rcx, [r11]
    mov r10, r13
    add r10, r9
    mov [r10], rax
    mov r10, r13
    add r10, rax
    mov rax, rcx
    mov [r10], rax
    inc r9
    cmp r9, rbx
    jl .loop3

    mov r9, rbx
.loop4:
    mov rax, r9
    shl rax, cl
    mov r8, rax
    mov rax, r9
    add rax, rbx
    shl rax, cl
    mov r10, r12
    add r10, r8
    mov r11, r10
    add r11, rdx
    mov rax, [r10]
    mov rcx, [r11]
    mov r10, r13
    add r10, r8
    mov [r10], rax
    mov r10, r13
    add r10, rax
    mov rax, rcx
    mov [r10], rax

    mov rax, r9
    mov rcx, r14
    shr rcx, 1
    shl rax, cl
    mov r8, r15
    shr r8, 1
    xor rax, r8
    mov r8, rax
    mov rax, r12
    add rax, r11
    mov r9, r13
    add r9, r8
    mov r10, rdx
    add r10, r8
    mov r11, rsi
    add r11, r8
    call fft_recursion

    inc r9
    cmp r9, rbx
    jl .loop4

    add rsp, 56
    pop rsi
    pop r15
    pop r14
    pop r13
    pop r12
    pop rbp
    pop rbx
    ret

.base_case:
    mov rdx, r15
    sub rdx, 1
    mov rax, r12
    mov rcx, r13
    movsd xmm0, [w_real]
    movsd xmm1, [w_imag]
    xor r9, r9
.loop5:
    mov rax, r9
    shl rax, cl
    mov r8, rax
    mov rax, r12
    add rax, r8
    movsd xmm2, [rax]
    mov rax, r13
    add rax, r9
    movsd [rax], xmm2
    inc r9
    cmp r9, rdx
    jle .loop5

    movsd xmm2, [rsp+8]
    movsd xmm3, [rsp+16]
    mulsd xmm2, xmm0
    mulsd xmm3, xmm1
    addsd xmm2, xmm3
    movsd [rsp+8], xmm2

    movsd xmm3, [rsp+24]
    movsd xmm4, [rsp+32]
    mulsd xmm3, xmm1
    mulsd xmm4, xmm0
    subsd xmm3, xmm4
    movsd [rsp+24], xmm3

    mov rax, r13
    movsd xmm2, [rax]
    movsd xmm3, [rsp+8]
    addsd xmm2, xmm3
    movsd [rax], xmm2

    mov rax, r13
    add rax, rdx
    movsd xmm2, [rax]
    movsd xmm3, [rsp+24]
    addsd xmm2, xmm3
    movsd [rax], xmm2

    ret

该代码使用了递归算法，通过划分子问题并递归求解来完成 FFT。基本情况下，当输入序列长度为 2 时，直接计算 FFT 并返回结果；否则，将输入序列划分为两个长度为 n/2 的子序列，分别计算它们的 FFT，并将它们合并为长度为 n 的序列的 FFT。

该代码依赖于 w_real 和 w_imag 数组，它们包含了旋转因子的实部和虚部。它还需要一个额外的栈空间来保存递归调用的参数和临时变量。