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