1

離散フーリエ変換を計算するサブルーチンを含む 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. のディラックが得られます。その上、スペクトルは奇妙に見えます...これが私が得たものです:

スペクトラム

4

1 に答える 1

0

実数値シーケンスの離散フーリエ変換を計算するルーチンでよくあることですが、結果として得られる複素数値スペクトルは、スペクトルの半分だけが返されます。値を元の 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、...) の部分を 2 つの別々のグラフとしてプロットする必要があるでしょう。または振幅と位相（ここでも2つの別々のグラフとして）。

于 2015-02-05T16:25:51.517 に答える