Geology Reference
In-Depth Information
40 J=J+M
10 CONTINUE
C If input sequence is of unit length, exit.
IF(LG.EQ.1) GO TO 70
C Start calculation of DFTs by doubling.
C Set initial DFT length as power of 2.
L=1
50 LT2=L*2
DO 60 M=1,L
MM1=M-1
C Find argument of W raised to required power.
ARGM=DFLOAT(MM1)*ARG/DFLOAT(L)
C Find W to required power.
W=DCMPLX(DCOS(ARGM),DSIN(ARGM))
C Find double length transforms.
DO 60 K=M,LG,LT2
HOLD=W*G(K+L)
G(K+L)=G(K)-HOLD
G(K)=G(K)+HOLD
60 CONTINUE
C Update DFT length.
L=LT2
C Test if transform is complete.
IF(L.LT.LG) GO TO 50
70
CONTINUE
RETURN
END
2.3.4 Generalisation of the FFT
The fast Fourier transform presented in the previous section applies only to
sequences with lengths equal to an integral power of two. The sequences could
be padded with zeros to the required length, but the resultant change in record
length would alter the sample interval in frequency for the DFT, and in time for
the inverse DFT. Instead, in this section, we present a form of the FFT valid for
sequences of arbitrary length.
Our approach is based on the method of Bluestein (1970), which transforms the
DFT into a convolution. Beginning with expression (2.169) for the DFT, we have
the identities,
N
1
N
1
x j W kj + ( j 2
k 2 ) /2
j 2
+ k 2
x j W kj
G
=
=
k
j
=
0
j
=
0
W k 2 /2 N 1
x j W j 2 /2
W ( k j ) 2 /2
=
·
.
(2.176)
j = 0
The last sum on the right-hand side may be interpreted as the convolution of the
N -length sequence
W j 2 /2 x j ,
y j
=
j
=
0,1,..., N
1,
Search WWH ::




Custom Search