Geology Reference
In-Depth Information
with the (2
N
−
1)-length sequence
W
−
j
2
/2
z
j
=
,
j
=−
N
+
1,...,0,...,
N
−
1,
yielding the (3
N
−
2)-length sequence
h
k
,
k
=−
N
+
1,...,0,...,2
N
−
2.
The DFT is then expressible as
j
=
0
y
j
z
k
−
j
=
W
k
2
/2
N
−
1
W
k
2
/2
h
k
,
k
G
k
=
=
0,1,...,
N
−
1.
(2.177)
Thus, only the terms
h
0
,
h
1
,...,
h
N
−
1
of the convolved sequence are required in the
calculation of the transform.
The sequences y
j
and
z
j
are then zero padded to length an integral power of two,
equal to or larger than 3
N
2. The convolution is found, as outlined in Section 2.3.2,
by taking the inverse of the product of the DFTs of the padded sequences y
j
and
z
j
, using calls to the subroutine FFT2N.
The subroutine FFTN finds the DFT or its inverse for sequences of arbitrary
length
N
using the Bluestein convolution formulation;
N
can be any positive integer.
Again, the same subroutine can be used to find the inverse DFT by reversing the
sign of the exponent of
e
and dividing the input DFT sequence by
N
.
−
SUBROUTINE FFTN(LG,G,IS,L,WK1,WK2)
C
C FFTN finds the discrete Fourier transform of the sequence G of length LG
C for input IS=-1. For input IS=1, it finds the inverse of the DFT
C sequence G of length LG. LG can be any positive integer. The subroutine
C uses the method of Bluestein in which the transform is converted to
C a convolution. The sequences entering the convolution are padded with
C zeros to length L, an integral power of two, equal to or larger than 3LG-2.
C The convolution is then performed by finding the DFTs of each sequence
C using the subroutine FFT2N, multiplying them, and then finding the inverse
C of the resulting DFT. WK1 and WK2 are work area vectors, defined in
C the main programme, of length L or larger. On output, the transformed
C sequence replaces the input sequence.
C
IMPLICIT DOUBLE COMPLEX(A-H,O-Z)
DOUBLE PRECISION PI,ARG,ARGJ2,FACT,DC,DS
DIMENSION G(LG),WK1(L),WK2(L)
C Define pi.
PI=3.14159265358979D0
C Define common factor in argument of W.
ARG=PI*DFLOAT(IS)/DFLOAT(LG)
C Determine if DFT or inverse DFT.
FACT=1.D0
IF(IS.EQ.1) FACT=FACT/DFLOAT(LG)
C Find N, the power of two equal to, or next larger than, 3*LG-2,
C the length of the convolution.
LGT3M2=3*LG-2
N=0
I=1
Search WWH ::
Custom Search