使用 fftpack5.1 进行傅里叶变换的问题
Trouble with Fourier Transform using fftpack5.1
我在 Fortran 90 中使用 FFTPACK5.1 时遇到问题,它包含用于计算离散傅里叶变换的子例程。我设法安装它并使用例程,但是当我使用频率为 A 的简单正弦波检查是否一切正常时,我得到一个不在 A 处的非零系数(频率space,在频谱中)但在 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. 的狄拉克。此外,频谱看起来很奇怪......这是我得到的:
作为计算实值序列的离散傅里叶变换的典型例程,得到的复值谱仅返回一半的谱。为了将值放入原始的 N 元素数组中,仅返回第一个值(始终为实数)的实部。类似地,在 n
的偶数值的情况下,返回第 n
/2 个复数值(也始终为实数)的实部。对于所有其他复数值,实部和虚部交错。
所以偶数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
其中增量为 fech/n
。
要绘制复数值,您可能想要绘制实部(索引 1,2,4,6,...)和虚部(索引 3,5,7,...)作为两个单独的图表,或振幅和相位(再次作为两个单独的图表)。
我在 Fortran 90 中使用 FFTPACK5.1 时遇到问题,它包含用于计算离散傅里叶变换的子例程。我设法安装它并使用例程,但是当我使用频率为 A 的简单正弦波检查是否一切正常时,我得到一个不在 A 处的非零系数(频率space,在频谱中)但在 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. 的狄拉克。此外,频谱看起来很奇怪......这是我得到的:
作为计算实值序列的离散傅里叶变换的典型例程,得到的复值谱仅返回一半的谱。为了将值放入原始的 N 元素数组中,仅返回第一个值(始终为实数)的实部。类似地,在 n
的偶数值的情况下,返回第 n
/2 个复数值(也始终为实数)的实部。对于所有其他复数值,实部和虚部交错。
所以偶数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
其中增量为 fech/n
。
要绘制复数值,您可能想要绘制实部(索引 1,2,4,6,...)和虚部(索引 3,5,7,...)作为两个单独的图表,或振幅和相位(再次作为两个单独的图表)。