home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Professional
/
OS2PRO194.ISO
/
os2
/
prgramer
/
adaptor
/
examples
/
f77
/
fftold.f
< prev
next >
Wrap
Text File
|
1993-06-28
|
4KB
|
181 lines
Program Fast_Fourier_Transform
C Programm zur Berechnung der Potenzen der Einheitswurzeln
C N Problemgroesse; 2er Potenz
C N2 Anzahl der zu berechnenden Potenzen; N2 = N/2
C Direkt Fouriersynthese bzw. Fourieranalyse (1/-1)
C W(N) Feld zur Aufnahme der Potenzen der Einheitswurzeln
integer N,N2,k
complex W(:)
cmf$ layout W(:serial)
complex C(:)
real F(:), maxf
real Direkt
print * ,'Input k : (N = 2**k)'
read * ,k
N = 2**k
N2 = 2**(k-1)
allocate(W(0:N2-1))
allocate(C(0:N-1))
allocate(F(0:N-1))
c
c random initialization
c
call cmf_random (f)
c
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i=0,N-1
c(i) = f(i)
end do
c
Direkt = -1
call Presort (N,k,C)
call Roots (W,N2,N,Direkt)
call Transform(C,W,N,N2,k)
c
c inverse fft
Direkt = 1
call Presort (N,k,C)
call Roots (W,N2,N,Direkt)
call Transform(C,W,N,N2,k)
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i=0,n-1
f(i) = abs (c(i)/n - f(i))
end do
maxf = 0.0
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, N-1
reduce (maxval, maxf, f(i))
end do
print *,'Max = ', maxf
deallocate(F)
deallocate(C)
deallocate(W)
end
C *******************************************************
C -------------------------------------------------------
C ------------ Umordnung der Speicherplaetze ------------
C ------------ nach Cooley und Tukey ------------
C -------------------------------------------------------
subroutine Presort (N,k,C)
complex C(0:N-1)
complex H(0:N-1)
Integer kl,kn
integer hi, hj, i, index(0:n-1), ix
c initialization of index
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, n-1
index(i) = i
end do
c executing bit reverse for all elements of index
do L=1,k
kl = 2**(L-1)
kn = N / (kl + kl)
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, n-1
ix = index(i)
c hi = ix / 2
hi = ishft (ix,-1)
c hj = hi / kn
hj = ishft (hi,L-k)
c hi = mod (hi,kn)
hi = iand (hi,kn-1)
index (i) = hi + 2*kn*hj + mod(ix,2)*kn
end do
end do
H = C(index)
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, n-1
C(i) = H(i)
end do
end
C ********************************************************
C -------------------------------------------------------
C ----- Berechnung der Potenzen der Einheitswurzeln -----
C ----- der Alg. spiegelt eine Rekurrenz erster -----
C ----- Ordnung wieder. -----
C -------------------------------------------------------
subroutine Roots (W,N2,N,Direkt)
complex W(0:N2-1),omega
cmf$ layout W(:serial)
real phi
real Direkt
phi = Direkt * 8.0 * Atan(1.0) / real(n)
omega = cmplx(cos(phi),sin(phi))
W(0) = 1
do i = 1,N2-1
W(i) = W(i-1) * omega
end do
end
C ********************************************************
C **********************************************************
C --------------------------------------------
C ----- Transformation nach Cooley-Tukey -----
C --------------------------------------------
subroutine Transform(C,W,N,N2,k)
C Feld c repraesentiert die Fourierkoeffizienten
C Feld w beinhaltet die Potenzen der Einheitswurzeln
complex C(0:N-1),W(0:N2-1),H(0:N-1)
cmf$ layout W(:serial)
complex t
integer kl,kn,kl1,hi,index(0:N-1),ni
do L=1,k
kl = 2**(L-1)
kn = N2/kl
kl1 = 2 * kl
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, n-1
index(i) = ieor (i,kl)
end do
H = C(index)
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, n-1
c hi = mod (i,kl1)
hi = iand (i,kl1-1)
c ni = mod (i,kl)
ni = iand (i,kl-1)
if (hi .ge. kl) then
C(i) = H(i) - C(i) * w(ni*kn)
else
C(i) = C(i) + w(ni*kn) * H(i)
end if
end do
end do
end