fortran90实现傅立叶变换的函数

以下是Fortran90实现傅里叶变换（FFT）的函数示例：

subroutine fft(x, y, n, sign)
  implicit none
  integer, intent(in) :: n, sign
  real, dimension(n), intent(inout) :: x, y
  real, dimension(n) :: xr, xi, yr, yi
  integer :: i, j, k, m, istep
  real :: theta, wtemp, wpr, wpi, wr, wi
  
  ! Bit-reverse the input array
  j = 1
  do i = 1, n
    if (j > i) then
      xr(j) = x(i)
      xi(j) = y(i)
      xr(i) = x(j)
      xi(i) = y(j)
    else
      xr(i) = x(i)
      xi(i) = y(i)
    end if
    m = n/2
    do while (m >= 1 .and. j > m)
      j = j - m
      m = m/2
    end do
    j = j + m
  end do
  
  ! Perform the FFT
  m = 1
  do while (m < n)
    istep = 2*m
    theta = sign*2.0*acos(-1.0)/m
    wtemp = sin(0.5*theta)
    wpr = -2.0*wtemp*wtemp
    wpi = sin(theta)
    wr = 1.0
    wi = 0.0
    do k = 1, m
      do i = k, n, istep
        j = i + m
        wtemp = wr*xr(j) - wi*xi(j)
        xi(j) = wr*xi(j) + wi*xr(j)
        xr(j) = wtemp
        yr(i) = xr(i) + xr(j)
        yi(i) = xi(i) + xi(j)
        yr(j) = xr(i) - xr(j)
        yi(j) = xi(i) - xi(j)
      end do
      wtemp = wr
      wr = wr*wpr - wi*wpi + wr
      wi = wi*wpr + wtemp*wpi + wi
    end do
    m = istep
  end do
  
  ! Normalize the output array
  if (sign == 1) then
    x = yr/n
    y = yi/n
  else
    x = xr/n
    y = xi/n
  end if
end subroutine fft

该函数接受一个长度为n的实数数组x和虚数数组y作为参数，并将其进行FFT变换。参数sign指定正向（1）或反向（-1）变换。请注意，FFT是一种高效的算法，可以在O(nlogn)时间内完成变换，但它要求输入数组的长度为2的幂。在实际使用中，通常会将数组长度补足为2的幂，或使用其他的FFT变种来处理非2的幂长度的数组。