! Last change: DOS 3 Aug 2000 5:55 pm ! *** copyright 2000 *** ! *** filename chirpz.f95 *** John F. Monahan ** ! ********************** program pchirpz ! TEST PROGRAM FOR CHIRPZ -- FAST FOURIER TRANSFORM ! computes FFT for any length series is O(N log N) ! using only "Power of 2" algorithm FFT2N implicit none complex, dimension(64) :: a,b,wpower complex s real t real, parameter :: twopi = 6.2831853 integer j,jj,k,kk,length,lngth2 ! 21 format(4x,2f16.6) 22 format(' Test of Chirp-Z to compute FFT -- original vector' ) 23 format(' Transformed vector directly, real and imaginary parts' ) 24 format(' Transformed vector by chirpz, real and imaginary parts' ) 25 format(' Inverse of transform -- is it same as original?' ) 26 format(i4,2f16.6) ! open( unit=6, file='chirpz.out' ) ! ! work on a series of length 7 length = 7 lngth2 = 4 ! load original vector and fill table do j = 1,length a(j) = real(j) ! table powers of w ! wpower(j) = exp(-2*pi*i*(j-1)/length) t = twopi*real(j-1)/real(length) wpower(j) = cmplx( cos(t), -sin(t) ) end do ! loop on j ! write out original vector write(6,22) write(6,21) (a(j),j=1,length) write(6,23) ! compute DFT directly do jj = 1,length j = jj - 1 s = 0. do kk = 1,length k = kk - 1 s = s + a(kk)*wpower( mod(j*k,length) + 1 ) end do ! loop on kk write(6,26) j,s b(jj) = s end do ! loop on i write(6,23) write(6,21) (b(j),j=1,length) ! ! compute DFT using Chirp-Z call chirpz(a,length,lngth2,1.) write(6,24) write(6,21) (a(j),j=1,length) ! transform back with Chirp-Z call chirpz(a,length,lngth2,-1.) write(6,25) write(6,21) (a(j),j=1,length) stop end program pchirpz subroutine chirpz(a,length,lngth2,sgn) ! computes the discrete Fourier transform using chirpz transform ! ! a complex vector vector of length length to be transformed ! ! length integer length of vector to be transformed ! ! lngth2 integer integer satisfying (length le 2**(lngth2-1)) ! ! sgn real indicates forward or inverse transform ! sgn = 1. then forward transform ! sgn = -1. then inverse transform ! ! *** required subprogram *** ! subroutine fft2n(a,n,sgn) -- computes FFT by "Power of 2" algorithm ! implicit none integer, intent(in) :: length,lngth2 complex, dimension(length) :: a complex, dimension(2**lngth2) :: x,y real, parameter :: twopi = 6.2831853 integer nstar,j,jj real sgn,t,c,s,fl,f2l ! ! lngth2 is value of n such that ! need 2**(n-2) < length < 2**(n-1) nstar = 2 ** lngth2 ! so now 2*length < twon ! fl = real(length) f2l = real(2*length) ! set up both vectors for convolution do jj = 1,nstar j = jj - 1 t = -sgn*twopi * real( mod( j*j, 2*length ) )/ f2l c = cos(t) s = sin(t) ! x is a(j)*w(j*j/2) if( j .lt. length ) then x(jj) = a(jj) * cmplx( c, s ) else x(jj) = ( 0., 0. ) ! pad with zeros end if ! ( j .lt. length ) ! construct y = w(-j*j/2) if( j .lt. length) then y(jj) = cmplx( c, -s ) else y(jj) = ( 0., 0.) ! pad with zeros end if ! ( j .lt. length ) if( j .gt. nstar - length ) then ! circular convolution j = nstar - j t = -sgn*twopi * real( mod( j*j, 2*length ) )/ f2l c = cos(t) s = sin(t) y(jj) = cmplx( c, -s ) end if ! ( j .gt. nstar - length ) end do ! loop on jj ! ! ! transform both x and y call fft2n(x,lngth2,1.) call fft2n(y,lngth2,1.) ! ! multiply together -- elementwise vector product do jj = 1,nstar x(jj) = x(jj) * y(jj) end do ! loop on jj ! ! transform back to get convolution call fft2n(x,lngth2,-1.) ! ! modify to get final result do jj = 1,length j = jj - 1 t = -sgn*twopi * real( mod( j*j, 2*length ) )/ f2l c = cos(t) s = sin(t) a(jj) = x(jj) * cmplx( c, s) if( sgn .eq. -1. ) a(jj) = a(jj) / fl ! divide if inverse end do ! loop on jj ! all done return end subroutine chirpz SUBROUTINE FFT2N(A,N,SGN) ! ! SUBROUTINE FFT2N -- DISCRETE FOURIER TRANSFORM ! ! A COMPLEX VECTOR INPUT: VECTOR TO BE TRANSFORMED ! OUTPUT: TRANSFORMED VECTOR ! N INTEGER INPUT: LOG(BASE 2) OF LENGTH OF VECTOR A ! LENGTH OF A IS 2 ** N ! SGN REAL INPUT: INDICATES FORWARD OR INVERSE TRANSFO ! SGN = 1. THEN FORWARD TRANSFORM ! SGN = -1. THEN INVERSE TRANSFORM ! ! DISCRETE FOURIER TRANSFORM OF A VECTOR A(J) IS GIVEN BY ! NN-1 ! T( A ) = SUM A(J) EXP( - I SGN 2 PI J K / NN ) ! J K=0 ! ! WHERE I = SQRT( -1 ) ! NN = 2 ** N ! PI = 3.1415923565... ! ! NOTE THAT THE INDEXING OF A IS SHIFTED BY ONE -- ! A(0) STORED IN A(1), ... , A( NN-1 ) STORED IN A( NN ) ! ! J F MONAHAN (NOV,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695-8203 ! RECODED DECEMBER 1999, APRIL 2000 FOR FORTRAN 95 ! ! THIS IS A 'POWER OF 2' ALGORITHM BASED ON THE PRESENTATION IN ! ! A. V. AHO, J. E. HOPCROFT, AND J. D. ULLMAN (1974) THE DESIGN AND ! ANALYSIS OF COMPUTER ALGORITHMS, ADDISON-WESLEY, READING. ! ! ! BEGIN WITH DECLARATIONS IMPLICIT NONE INTEGER, INTENT(IN) :: N COMPLEX, DIMENSION(1:N), INTENT(IN OUT) :: A REAL, INTENT(IN) :: SGN REAL C,S COMPLEX, DIMENSION(25) :: WK COMPLEX WS,ALEFT,ARIGHT LOGICAL CARRY,NEW LOGICAL, DIMENSION(25) :: RBITS INTEGER I,J,L,LM,LL,M,MM,M2,NN ! ! COMPUTE WK(I) = EXP( - 2 PI IMAG SGN / 2**I ) ! ! SPECIAL VALUES FIRST WK(1) = CMPLX( -1., 0. ) WK(2) = CMPLX( 0., -SGN ) ! SKIP IF N IS ONLY 1 OR 2 IF( N .GT. 2 ) THEN ! AVOID SINES AND COSINES -- JUST USE SQUARE ROOTS C = 0. S = 1. DO I = 3,N C = SQRT( ( 1. + C ) / 2. ) S = S / ( 2.*C ) WK(I) = CMPLX( C, -SGN*S ) END DO ! LOOP ON I END IF ! ( N .GT. 2 ) ! ! THIS IS THE MAIN LOOP OF THE ALGORITHM ! DO MM = 1,N ! CONVERT TO REVERSE LOOP M = N - MM ! M STEPS DOWN FROM N-1 TO 0 ! ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I ! M2 M2 = 2**M ! LM CONTROLS THE NUMBER OF SUBLOOPS LM = 2**(MM-1) ! DO LL = 1,LM ! L IS THE MAJOR INDEX VARIABLE L = 2*M2*(LL-1) ! INITIALIZE TO GET W ** REV( L / M2 ) WS = CMPLX(1., 0.) CARRY = .TRUE. ! START AT 2 (SHIFT) SINCE L/M2 IS ALWAYS EVEN DO I = 2,N NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) WS = WS * WK(I) END DO ! LOOP ON I ! WS IS NOW W ** REV( L / M2 ) ! NOW DO THE REAL COMPUTATION DO J = 1,M2 ALEFT = A( L + J ) ARIGHT = WS*A( L + J + M2 ) A( L + J ) = ALEFT + ARIGHT A( L + J + M2 ) = ALEFT - ARIGHT END DO ! LOOP ON J ! DONE WITH LOOP ON L END DO ! LOOP ON LL ! DONE WITH MAIN LOOP ON M END DO ! LOOP ON MM ! ! SCRAMBLE THE NUMBERS -- SWITCH A(J) AND A( REV(J) ) ! BUT DON'T DO IT TWICE ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I NN = 2 ** N ! LOOP THROUGH ALL THE NUMBERS DO LL = 1,NN ! DO THE BIT REVERSAL LM = 0 M2 = NN CARRY = .TRUE. ! START AT 1 HERE -- USE ALL THE NUMBERS DO I = 1,N M2 = M2 / 2 NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) LM = LM + M2 END DO ! LOOP ON I ! DON'T SWITCH TWICE IF( LL-1 .LT. LM ) THEN ALEFT = A(LL) ARIGHT = A(LM+1) A(LL) = ARIGHT A(LM+1) = ALEFT END IF ! ( LL-1 .LT. LM ) END DO ! LOOP ON LL ! IF INVERSE TRANSFORM, REMEMBER TO DIVIDE BY NN IF( SGN .GT. 0. ) RETURN C = REAL(NN) DO I = 1,NN A(I) = A(I) / C END DO ! LOOP ON I RETURN END SUBROUTINE FFT2N ! *** end of filename chirpz.f95 *********************