Make your own free website on Tripod.com
Chirp z-Transform Algorithm Program

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: TEST PROGRAM FOR CHIRP Z-TRANSFORM SUBROUTINE
C               USING A DELAY IMPULSE
C AUTHOR:       L R RABINER
C               BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974
C
C INPUT:        NDATA IS THE NUMBER OF DATA POINTS
C                   1 <= NDATA <= 512
C               NOPTS IS THE NUMBER OF OUTPUT POINTS
C                   1 <= NOPTS <= 512
C               IDEL IS THE DELAY OF THE IMPULSE IN SAMPLES
C               FS IS THE SAMPLING FREQUENCY IN HZ
C               SIG0 IS THE INITIAL VALUE OF SIGMA IN HZ
C               OME0 IS THE INITIAL VALUE OF OMEGA IN HZ
C               DLTSIG IS THE INCREMENT IN SIGMA IN HZ
C               DLTOMG IS THE INCREMENT IN OMEGA IN HZ
C               IMD REQUESTS ADDITIONAL RUNS
C                   IMD = 1  NEW RUN
C                   IMD = 0  TERMINATES PROGRAM
C-----------------------------------------------------------------------
C
      COMMON WTR(1024), WTI(1024), XR(1024), XI(1024)
      INTEGER TTI, TTO
C
C DEFINE I/O DEVICE CODES
C INPUT: INPUT TO THIS PROGRAM IS USER-INTERACTIVE
C        THAT IS - A QUESTION IS WRITTEN ON THE USER
C        TERMINAL (TTO) AND THE USER TYPES IN THE ANSWER.
C
C OUTPUT: ALL OUTPUT IS WRITTEN ON THE STANDARD
C         OUTPUT UNIT (LPT)
C
      TTI = I1MACH(1)
      TTO = I1MACH(4)
      LPT = I1MACH(2)
C
  10  WRITE (TTO,9999)
9999  FORMAT (23H NO OF DATA POINTS(I4)=)
      READ (TTI,9998) NDATA
9998  FORMAT (I4)
      WRITE (LPT,9997) NDATA
9997  FORMAT (19H NO OF DATA POINTS=, I5)
      WRITE (TTO,9996)
9996  FORMAT (25H NO OF OUTPUT POINTS(I4)=)
      READ (TTI,9998) NOPTS
      WRITE (LPT,9995) NOPTS
9995  FORMAT (21H NO OF OUTPUT POINTS=, I5)
      WRITE (TTO,9994)
9994  FORMAT (33H DATA INPUT DELAY IN SAMPLES(I4)=)
      READ (TTI,9998) IDEL
      WRITE (LPT,9993) IDEL
9993  FORMAT (30H DELAYED IMPULSE AT SAMPLE NO=, I4)
      WRITE (TTO,9992)
9992  FORMAT (33H SAMPLING FREQUENCY IN HZ(F10.0)=)
      READ (TTI,9991) FS
9991  FORMAT (F10.0)
      WRITE (LPT,9990) FS
9990  FORMAT (20H SAMPLING FREQUENCY=, F10.3)
      WRITE (TTO,9989)
9989  FORMAT (28H INITIAL SIGMA VALUE(F10.0)=)
      READ (TTI,9991) SIG0
      WRITE (LPT,9988) SIG0
9988  FORMAT (24H INITIAL VALUE OF SIGMA=, F10.3)
      WRITE (TTO,9987)
9987  FORMAT (28H INITIAL OMEGA VALUE(F10.0)=)
      READ (TTI,9991) OME0
      WRITE (LPT,9986) OME0
9986  FORMAT (24H INITIAL VALUE OF OMEGA=, F10.3)
      WRITE (TTO,9985)
9985  FORMAT (20H DELTA SIGMA(F10.0)=)
      READ (TTI,9991) DLTSIG
      WRITE (LPT,9984) DLTSIG
9984  FORMAT (35H INCREMENT IN SIGMA ALONG CZT PATH=, F10.3)
      WRITE (TTO,9983)
9983  FORMAT (20H DELTA OMEGA(F10.0)=)
      READ (TTI,9991) DLTOMG
      WRITE (LPT,9982) DLTOMG
9982  FORMAT (35H INCREMENT IN OMEGA ALONG CZT PATH=, F10.3)
C
C CREATE DELAYED INPUT SAMPLE WITH DELAY IDEL SAMPLES
C
      DO 20 I=1,NDATA
        XR(I) = 0.
        XI(I) = 0.
  20  CONTINUE
      XR(IDEL+1) = 1.
      NFFT = NDATA + NOPTS
      DO 30 I=1,10
        NTEST = 2**I
        IF (NTEST.GE.NFFT) GO TO 40
  30  CONTINUE
      WRITE (TTO,9981)
9981  FORMAT (18H N TOO BIG FOR FFT)
      STOP
  40  NFFT = NTEST
      CALL CZT(XR, XI, NDATA, NOPTS, DLTSIG, DLTOMG, WTR, WTI, SIG0,
     *    OME0, 0, NFFT, FS)
      WRITE (LPT,9980)
9980  FORMAT (//)
      WRITE (LPT,9979) (XR(I),I=1,NOPTS)
9979  FORMAT (29H REAL PART OF OUTPUT SEQUENCE/(5E14.5))
      WRITE (LPT,9978) (XI(I),I=1,NOPTS)
9978  FORMAT (34H IMAGINARY PART OF OUTPUT SEQUENCE/(5E14.5))
      WRITE (TTO,9977)
9977  FORMAT (27H MORE DATA(1=YES, 0=NO)(I1))
      READ (TTI,9976) IMD
9976  FORMAT (I1)
      WRITE (LPT,9975)
9975  FORMAT (1H1)
      IF (IMD.EQ.1) GO TO 10
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CZT
C CHIRP Z-TRANSFORM
C REFERENCE--L R RABINER, R W SCHAFER, C M RADER--BELL SYSTEM
C TECHNICAL JOURNAL--MAY 1969, PP 1249-1291
C-----------------------------------------------------------------------
C
      SUBROUTINE CZT(XR, XI, NDATA, NOPTS, DLTSIG, DLTOMG, WTR, WTI,
     *    SIGO, OME, NTR, NFFT, FS)
      DIMENSION WTR(1), WTI(1), XR(1), XI(1)
C
C     XR = ARRAY OF SIZE NFFT
C          ON INPUT XR CONTAINS REAL PART OF INPUT DATA IN LOCATIONS
C          1 TO NDATA
C          ON OUTPUT XR CONTAINS REAL PART OF OUTPUT DATA IN LOCATIONS
C          1 TO NOPTS
C     XI = ARRAY OF SIZE NFFT
C          ON INPUT XI CONTAINS IMAGINARY PART OF INPUT DATA
C          ON OUTPUT XI CONTAINS IMAGINARY PART OF OUTPUT DATA
C  NDATA = NUMBER OF INPUT POINTS BEING TRANSFORMED
C  NOPTS = NUMBER OF OUTPUT VALUES BEING COMPUTED
C DLTSIG = INCREMENT IN SIGMA ALONG THE CZT CONTOUR--RELATIVE TO FS
C DLTOMG = INCREMENT IN OMEGA ALONG THE CZT CONTOUR--RELATIVE TO FS
C          NOTE THAT IF DLTSIG IS LESS THAN ZERO, THE CONTOUR SPIRALS
C          OUTSIDE THE UNIT CIRCLE IN THE Z-PLANE
C          WTR IS AN ARRAY OF SIZE NFFT WHICH HOLDS THE REAL PART OF THE
C          TRANSFORM OF W**(-N**2/2)
C          WTI IS AN ARRAY OF SIZE NFFT WHICH HOLDS THE IMAGINARY PART
C          OF THE TRANSFORM OF W**(-N**2/2)
C   SIG0 = INITIAL VALUE OF SIGMA OF CZT CONTOUR, RELATIVE TO FS
C   OME0 = INITIAL VALUE OF OMEGA OF CZT CONTOUR, RELATIVE TO FS
C    NTR = ROTATION FACTOR ON INPUT DATA TO MAKE INPUT LOOK LIKE IT
C          BEGINS AT N = -NTR INSTEAD OF N=0. SET NTR=0 TO AVOID
C          THIS OPTION
C   NFFT = SIZE OF FFT PERFORMED INTERNALLY IN CZT PROGRAM--GENERALLY
C          NFFT IS A POWER OF 2--NFFT MUST BE GREATER THAN (NDATA+NOPTS)
C     FS = SAMPLING FREQUENCY IN HZ
C
      DOUBLE PRECISION A, B, PHINC, RO, XIO, PI, W1INC, WR, WI, W2INC,
     *    WOHR, WOHI, W3INC, W4INC, WPR, WPI, WPOHR, WPOHI, WNR, WNI,
     *    WNSQR, WNSQI, RAD, ANG, WAR, WAI
      DOUBLE PRECISION XRT, XADRT, WOPR, WIPR, WINR, WINI, WJNR, WJNI,
     *    XASQ, WKNR, WKNI, WLNR, WLNI
      DOUBLE PRECISION XTR
C
C SET INITIAL CONSTANTS
C
      N3 = NFFT
      NFT = N3/2 + 1
      N4 = N3/2
      X9 = N3
C
C N0 FACTOR FOR SYMMETRY OF W**(-N**2/2)
C
      NADRT = (NDATA-NOPTS)/2
      XADRT = NADRT + 1
      NRT = NADRT - NTR
      XRT = NRT
      PI = 4.0D0*DATAN(1.0D0)
      RO = SIGO
      XIO = OME
C
C COMPUTE W**(N**2/2) USING RECURSION RELATIONS
C OBTAIN INITIAL CONSTANTS
C
      A = DLTSIG
      B = DLTOMG
      PHINC = -PI*B*2.D0/FS
      W1INC = DEXP(A*PI/FS)
C
C W STORED IN WR AND WI
C
      WR = W1INC*DCOS(PHINC)
      WI = W1INC*DSIN(PHINC)
      W2INC = DSQRT(W1INC)
C
C W**(1/2) STORED IN WOHR AND WOHI
C
      WOHR = W2INC*DCOS(PHINC/2.0D0)
      WOHI = W2INC*DSIN(PHINC/2.0D0)
      W3INC = 1.0D0/W1INC
      W4INC = 1.0D0/W2INC
C
C W**(-1) STORED IN WPR AND WPI
C
      WPR = W3INC*DCOS(-PHINC)
      WPI = W3INC*DSIN(-PHINC)
C
C W**(-1/2) STORED IN WPOHR AND WPOHI
C
      WPOHR = W4INC*DCOS(-PHINC/2.0D0)
      WPOHI = W4INC*DSIN(-PHINC/2.0D0)
      DO 10 I=1,N4
        WTR(I) = 0.
        WTI(I) = 0.
  10  CONTINUE
C
C D0=W**(-1/2)   D0 STORED IN WNR AND WNI
C
      WNR = WPOHR
      WNI = WPOHI
C
C C0=1
C
      WNSQR = 1.0D0
      WNSQI = 0.0D0
C
C SOLVE RECURSION RELATION FOR W**(-N**2/2)
C
      CALL RECUR(WTR, WTI, WNR, WNI, WNSQR, WNSQI, WPR, WPI, N3)
C
C COMPUTE TRANSFORM OF W**(-N**2/2)
C
      CALL FFT842(0, N3, WTR, WTI)
C
C COMPUTE A**(-N)*W**(N**2/2)
C INITIAL CONSTANTS
C
      RAD = DEXP(RO*PI/FS)
      ANG = -XIO*PI*2.D0/FS
C
C A STORED IN WAR AND WAI
C
      WAR = RAD*DCOS(ANG)
      WAI = RAD*DSIN(ANG)
      XASQ = XADRT*XADRT/2.0D0
      WOPR = W1INC**(XASQ)
      WOPI = W1INC**(.5D0-XADRT)
      WKNR = WOPR*DCOS(PHINC*XASQ)
      WKNI = WOPR*DSIN(PHINC*XASQ)
      WLNR = WOPI*DCOS(PHINC*(.5D0-XADRT))
      WLNI = WOPI*DSIN(PHINC*(.5D0-XADRT))
      WNSQR = (WAR*WKNR+WAI*WKNI)/(RAD*RAD)
      WNSQI = (WAR*WKNI-WAI*WKNR)/(RAD*RAD)
      WNR = WAR*WLNR - WAI*WLNI
      WNI = WAI*WLNR + WAR*WLNI
C
C WEIGHT INPUT DATA BY A**(-N)*W**(N**2/2)
C
      CALL DECUR(XR, XI, WNR, WNI, WNSQR, WNSQI, WR, WI, NDATA)
C
C PAD DATA WITH ZEROS
C
      N5 = NDATA + 1
      DO 20 I=N5,N3
        XR(I) = 0.
        XI(I) = 0.
  20  CONTINUE
C
C TRANSFORM INPUT USING FFT
C
      CALL FFT842(0, N3, XR, XI)
C
C MULTIPLY FFTS
C
      DO 30 I=1,N3
        J = I
        IF (I.GT.NFT) J = N3 + 2 - I
        XT = XR(I)*WTR(J) - XI(I)*WTI(J)
        XI(I) = -XR(I)*WTI(J) - XI(I)*WTR(J)
        XR(I) = XT
  30  CONTINUE
C
C INVERSE TRANSFORM TO GET TIME SEQUENCE
C
      CALL FFT842(0, N3, XR, XI)
C
C SHUFFLE DATA
C
      IF (NADRT.EQ.0) GO TO 70
      IF (NADRT.LT.0) GO TO 50
      DO 40 I=1,NOPTS
        J = I + NADRT
        XR(I) = XR(J)
        XI(I) = XI(J)
  40  CONTINUE
      GO TO 70
  50  DO 60 I=1,NOPTS
        J = NOPTS + 1 - I
        K = J + NADRT
        IF (K.LE.0) K = K + N3
        XR(J) = XR(K)
        XI(J) = XI(K)
  60  CONTINUE
  70  CONTINUE
      DO 80 I=1,NOPTS
        XR(I) = XR(I)/X9
        XI(I) = -XI(I)/X9
  80  CONTINUE
C
C COMPUTE POST WEIGHTING W**(-K**2/2) INITIAL CONSTANTS FOR RECURSION
C
      WOPR = W1INC**(.5D0-XRT)
      WOPI = 1.0D0/WOPR
      WINR = WOPR*DCOS(PHINC*(.5D0-XRT))
      WINI = WOPR*DSIN(PHINC*(.5D0-XRT))
      WJNR = WOPI*DCOS(PHINC*(XRT-.5D0))
      WJNI = WOPI*DSIN(PHINC*(XRT-.5D0))
      XTR = NTR
      WOPR = RAD**(-XTR)
      WAR = WOPR*DCOS(-ANG*XTR)
      WAI = WOPR*DSIN(-ANG*XTR)
      WNR = WJNR
      WNI = WJNI
      WNSQR = WINR*WAR - WINI*WAI
      WNSQI = WINR*WAI + WINI*WAR
C
C POST WEIGHT DATA
C
      CALL DECUR(XR, XI, WNR, WNI, WNSQR, WNSQI, WR, WI, NOPTS)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: RECUR
C RECURSION RELATION TO GIVE W**(-N**2/2)
C-----------------------------------------------------------------------
C
      SUBROUTINE RECUR(WR, WI, WNR, WNI, WNSQR, WNSQI, VR, VI, NFFT)
      DIMENSION WR(1), WI(1)
      DOUBLE PRECISION WNR, WNI, WNSQR, WNSQI, VR, VI
      DOUBLE PRECISION XT
C
C C(N+1) = C(N)*D(N)
C D(N+1) = D(N)*W**(-1)
C W**(-N**2/2) = D(N)
C WNR AND WNI ARE THE D TERMS
C WNSQR AND WNSQI ARE THE C TERMS
C WR AND WI STORE THE W**(-N**2/2) RESULTS
C
      WR(1) = 1.
      WI(1) = 0.
      N3 = NFFT/2 + 1
      DO 10 I=2,N3
C
C C(N+1)=C(N)*D(N)
C
        XT = WNSQR*WNR - WNSQI*WNI
        WNSQI = WNSQR*WNI + WNSQI*WNR
        WNSQR = XT
C
C D(N+1)=D(N)*W
C
        XT = WNR*VR - WNI*VI
        WNI = WNR*VI + WNI*VR
        WNR = XT
C
C STORE W**(-N**2/2)
C
        WR(I) = WNSQR
        WI(I) = WNSQI
        J = NFFT + 2 - I
        WR(J) = WR(I)
        WI(J) = WI(I)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DECUR
C RECURSION RELATION TO GIVE Y(N)=X(N)*A**(-N)*W**(N**2/2)
C-----------------------------------------------------------------------
C
      SUBROUTINE DECUR(XR, XI, WNR, WNI, WNSQR, WNSQI, VR, VI, N)
C
C Y(N) OVERWRITES X(N)
C X(N) STORED IN XR AND XI
C C(N+1) = C(N)*D(N)
C D(N+1) = D(N)*W
C WNR AND WNI STORE C(N)--INITIALLY SET TO C0
C WNSQR AND WNSQI STORE D(N)--INITIALLY SET TO D0
C
      DIMENSION XR(1), XI(1)
      DOUBLE PRECISION WNR, WNI, WNSQR, WNSQI, VR, VI
      DOUBLE PRECISION XT
      DO 10 J=1,N
C
C C(N+1)=C(N)*D(N)
C
        XT = WNSQR*WNR - WNSQI*WNI
        WNSQI = WNSQR*WNI + WNSQI*WNR
        WNSQR = XT
C
C D(N+1)=D(N)*W
C
        XT = WNR*VR - WNI*VI
        WNI = WNR*VI + WNI*VR
        WNR = XT
C
C WEIGHT INPUT DATA
C
        XS = XR(J)*WNSQR - XI(J)*WNSQI
        XI(J) = XR(J)*WNSQI + XI(J)*WNSQR
        XR(J) = XS
  10  CONTINUE
      RETURN
      END
Return to Main Software Page