离散傅里叶变换似乎打印出错误的答案？

Question

我正在尝试编写一个程序来计算一组给定数据的离散傅立叶变换。我采样了一个正弦波，所以我的设置是 (pi/2,2*pi,3*pi/2,2*pi)。这是我的程序：

program DFT
implicit none
integer :: k, N, x, y, j, r, l
integer, parameter :: dp = selected_real_kind(15,300)
real, allocatable,dimension(:) :: h, rst
integer, dimension(:,:), allocatable :: W
real(kind=dp) :: pi


open(unit=100, file="dft.dat",status='replace')

N = 4

allocate(h(N))
allocate(rst(N))
allocate(W(-N/2:N/2,1:N))

pi = 3.14159265359

do k=1,N
  h(k) = k*(pi*0.5)
end do


do j = -N/2,N/2
do k = 1, N

    W(j,k) = EXP((2.0_dp*pi*cmplx(0.0_dp,1.0_dp)*j*k)/N)

end do
end do

rst = matmul(W,h)
!print *, h, w
write(100,*) rst

end program

这首先打印出数组：

0.00000000  0.00000000  15.7079639  0.00000000  0.00000000

使用在线计算器，结果应该是：

15.7+0j  -3.14+3.14j  -3.14+0j  -3.14-3.14j

我也不知道为什么第一个条目太长了。

谁能看出为什么 3/4 的结果打印出 0？我注意到 15.7 出现在实际答案和我的结果中。

谢谢

Answer 1

Fortran 为您计算复数，但您必须将适当的变量声明为复数。尝试

complex, allocatable,dimension(:) :: rst
complex, dimension(:,:), allocatable :: W

Answer 2

尽管问题已经被回答并接受了，但是给出的程序有很多问题，我不得不说...

给定的输入不是正弦波，它是时间的线性函数。有点像基于 1 的斜坡输入。

对于 DFT，索引通常被认为来自 0:N-1，而不是 1:N。

对于 W，奈奎斯特频率表示为两次，如 -N/2 和 N/2。同样，给行编号 0:N-1 是正常的，顺便说一句，这就是为什么你的 rst 向量中有一个额外的输出。

pi 是双精度但只初始化为 12 位有效数字。很难判断 pi 的值是否有错字，这就是为什么许多人会使用 4*atan(1.0_dp) 或 acos(-1.0_dp) 的原因。

请注意，h(N) 实际上将作为零时间输入结束，这是整个世界从零开始索引 DFT 向量的原因之一。

表达式 cmplx(0.0_dp,1.0_dp) 有点无用，因为如果第三个可选的 KIND= 参数不存在，CMPLX 内在函数总是 returns 单精度结果。作为复杂文字，(0.0_dp,1.0_dp) 将是双精度。但是，您也可以使用 (0,1)，因为它完全可以用单精度表示，并且在乘以其左侧不断增长的乘积时会转换为双精度。 2.0_dp 也可以成功地表示为 2，并且不会那么混乱。

表达式EXP((2.0_dp*pi*cmplx(0.0_dp,1.0_dp)*j*k)/N) 适用于逆 DFT，忽略归一化。因此，我会把整个事情写得更干净、更正确，如 EXP(-2*pi*(0,1)*j*k/N)。那么输出结果应该可以直接与在线计算器打印出来的结果进行比较。

离散傅里叶变换似乎打印出错误的答案？

Discrete Fourier Transform seems to be printing incorrect answers?

fortran

gfortran