home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 326.lha / KFFT_v1.1 / fft1.for < prev    next >
Text File  |  1989-12-23  |  2KB  |  69 lines
  1. C
  2. C********************************************************************
  3.       SUBROUTINE FFT1(A,M,N)
  4. C***********************************************************************
  5. C
  6. C *** DECIMATION IN TIME FFT (COOLEY-TUKEY ALGORITHM) FOR NATURALLY
  7. C *** ORDERED INPUTS
  8. C
  9. C *** FFT1 COMPUTES FORWARD DFT.  INVERSE DFT CAN BE COMPUTED WITH 
  10. C *** FFT1 BY FOLLOWING SEQUENCE:
  11. C ***     * COMPLEX CONJUGATE INPUT A
  12. C ***     * CALL FFT1
  13. C ***     * COMPLEX CONJUGATE RESULT A
  14. C ***     * NORMALIZE BY 1/N IF DESIRED
  15. C
  16. C ***   INPUTS:
  17. C ***     A          COMPLEX VECTOR OF LENGTH N CONTAINING NATURALLY
  18. C ***                ORDERED VALUES
  19. C ***     N          FFT SIZE, MUST BE POWER OF 2, 2**M
  20. C ***     M          LOG2(N)
  21. C
  22. C ***   OUTPUTS:
  23. C ***     A          FORWARD N-POINT DISCRETE TRANSFORM
  24. C
  25.       COMPLEX A,U,W,T
  26. C
  27.       DIMENSION A(N)
  28. C
  29.       DATA PI/3.14159265/
  30. C
  31.       N=2**M
  32. C
  33. C *** SCRAMBLE INPUT ARRAY (BIT REVERSE) SO THAT OUTPUTS WILL BE
  34. C *** NATURALLY ORDERED
  35.       NV2=N/2
  36.       NM1=N-1
  37.       J=1
  38.       DO 40 I=1,NM1
  39.         IF(I.LT.J) THEN
  40.           T=A(J)
  41.           A(J)=A(I)
  42.           A(I)=T
  43.         ENDIF
  44.         K=NV2
  45.    20   CONTINUE
  46.         IF(K.LT.J) THEN
  47.           J=J-K
  48.           K=K/2
  49.           GO TO 20
  50.         ENDIF
  51.         J=J+K
  52.    40 CONTINUE
  53. C
  54. C *** PERFORM DECIMATION IN TIME FFT
  55.       DO 60 L=1,M
  56.         LE=2**L
  57.         LE1=LE/2
  58.         U=CMPLX(1.,0.)
  59.         W=CMPLX(COS(PI/LE1),-SIN(PI/LE1))
  60.         DO 60 J=1,LE1
  61.           DO 50 I=J,N,LE
  62.             IP=I+LE1
  63.             T=A(IP)*U
  64.             A(IP)=A(I)-T
  65.             A(I)=A(I)+T
  66.    50     CONTINUE
  67.           
  68.       END
  69.