使用 fftpack5.1 的傅立叶变换问题

Question

我在 Fortran 90 中使用 FFTPACK5.1 时遇到问题，它包含计算离散傅立叶变换的子例程。我设法安装它并使用例程，但是当我检查频率是否为 A 的简单正弦波是否一切正常时，我得到一个不在 A 处的非零系数（在频率空间中，在光谱），但在 2A。频谱发生了变化，我不明白为什么。我几乎可以肯定（但我有疑问）我正确计算了频率轴步长：

N 是原始正弦波的点数，Fech 是我的采样频率，我将频率轴步长计算为 df(i)=Fech(i-1)/N。

我正在使用 rfft1f 例程，所以如果有人对此有经验并且知道我的问题，我将非常高兴了解这里的问题所在。

这是我的代码：

! n: number of samples in the discret signal
integer ( kind = 4 ), parameter :: n = 4096
real, parameter :: deuxpi=6.283185307
!frequence is the frequence of the signal
!fech is the frequence of sampling
real :: frequence,fech 
integer :: kk
! r is the signal i want to process
! t is the built time and f is the built frequency
real ( kind = 4 ) r(n),t(n),f(n)


!Arrays routines need to work (internal recipe):
real ( kind = 4 ), allocatable, dimension ( : ) :: work
real ( kind = 4 ), allocatable, dimension ( : ) :: wsave

!size of arrays wsave and work for internal recipe 
lensav = n + int ( log ( real ( n, kind = 4 ) ) / log ( 2.0E+00 ) ) + 4
lenwrk = n
allocate ( work(1:lenwrk) )
allocate ( wsave(1:lensav) )


! initializes rttft1f, wsave array
call rfft1i ( n, wsave, lensav, ier )


frequence=0.5
fech=20

! I built here the signal
do kk=1,n
t (kk) = (kk-1) / (fech)
f (kk) = fech*(kk-1) / n
r (kk) = sin( deuxpi * t(kk) * frequence  ) 
end do


!and I call the rfft1f to build the Discrete Fourier Transform:
call rfft1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

!I get back r which contains now the fourier coefficients and plot it

我期望一个简单的正弦波在频率为 0.5（cf 代码）处有一个狄拉克，但在频域中我得到了一个 1. 处的狄拉克。此外，光谱看起来很奇怪......这是我得到的：

Answer 1

作为计算实值序列的离散傅立叶变换的典型例程，所得到的复值频谱仅返回一半频谱。为了将这些值放入原始 N 元素数组，仅返回第一个值的实部（始终为实数）。同样，在n 的偶数值的情况下，返回n/2-th 复数值的实部（也总是实数）。对于所有其他复数值，实部和虚部是交错的。

所以即使是n，对应的频率由下式给出：

r(1)       -> 0
r(2), r(3) -> Delta
r(4), r(5) -> 2*Delta
...
r(n)       -> (n/2)*Delta

对于奇数n：

r(1)        -> 0
r(2), r(3)  -> Delta
r(4), r(5)  -> 2*Delta
...
r(n-1),r(n) -> ((n-1)/2) * Delta

其中 Delta 以fech/n 给出。

要绘制复数值，您可能希望将实部（索引 1,2,4,6,...）和虚部（索引 3,5,7,...）绘制为两个单独的图表，或幅度和相位（同样作为两个单独的图表）。