fortran90实现快速傅立叶变换的函数

下面是Fortran90实现快速傅立叶变换的函数：

subroutine fft(x, N, isign)
implicit none
integer, intent(in) :: N, isign
complex*16, intent(inout) :: x(N)
integer :: i, j, k, n1, n2, n3, n4, i1, i2, i3, i4
double precision :: c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, wtemp

! Bit-reverse permutation
j = 1
do i = 1, N-1
    if (i .lt. j) then
        t = x(j)
        x(j) = x(i)
        x(i) = t
    endif
    k = N/2
    do while (k .lt. j)
        j = j - k
        k = k/2
    enddo
    j = j + k
enddo

! FFT loop
n1 = 0
n2 = 1
do i = 0, log(N)/log(2)-1
    n3 = n2*2
    n4 = n1
    do j = 0, n2-1
        if (isign .lt. 0) then
            c1 = cos(2.0d0*3.141592653589793d0*j/n3)
            c2 = sin(2.0d0*3.141592653589793d0*j/n3)
        else
            c1 = cos(2.0d0*3.141592653589793d0*j/n3)
            c2 = -sin(2.0d0*3.141592653589793d0*j/n3)
        endif
        do k = j, N-1, n3
            i1 = k
            i2 = k + n2
            i3 = n4 + k
            i4 = i3 + n2
            h1r = dble(x(i1))
            h1i = dimag(x(i1))
            h2r = dble(c1*x(i2)) - dimag(c2*x(i2))
            h2i = dble(c2*x(i2)) + dimag(c1*x(i2))
            x(i1) = cmplx(h1r + h2r, h1i + h2i)
            x(i2) = cmplx(h1r - h2r, h1i - h2i)
            h1r = dble(x(i3))
            h1i = dimag(x(i3))
            h2r = dble(c1*x(i4)) - dimag(c2*x(i4))
            h2i = dble(c2*x(i4)) + dimag(c1*x(i4))
            x(i3) = cmplx(h1r + h2r, h1i + h2i)
            x(i4) = cmplx(h1r - h2r, h1i - h2i)
        enddo
    enddo
    n1 = n2
    n2 = n3
enddo

! Normalization
if (isign .lt. 0) then
    do i = 1, N
        x(i) = x(i)/dble(N)
    enddo
endif

end subroutine fft

此函数名为`fft`，有三个输入参数：`x`是一个长度为`N`的复数数组，`N`是数组的长度，`isign`是正负号标志，用于指定是进行正向变换还是反向变换。函数中使用了Cooley-Tukey快速傅立叶变换算法。