home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / HARM.FOR < prev    next >
Text File  |  1985-11-29  |  13KB  |  520 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE HARM
  5. C
  6. C        PURPOSE
  7. C           PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX
  8. C           THREE DIMENSIONAL ARRAY
  9. C
  10. C        USAGE
  11. C           CALL HARM (A,M,INV,S,IFSET,IFERR)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           A     - AS INPUT, A CONTAINS THE COMPLEX, 3-DIMENSIONAL
  15. C                   ARRAY TO BE TRANSFORMED.  THE REAL PART OF
  16. C                   A(I1,I2,I3) IS STORED IN VECTOR FASHION IN A CELL
  17. C                   WITH INDEX 2*(I3*N1*N2 + I2*N1 + I1) + 1 WHERE
  18. C                   NI = 2**M(I), I=1,2,3 AND I1 = 0,1,...,N1-1 ETC.
  19. C                   THE IMAGINARY PART IS IN THE CELL IMMEDIATELY
  20. C                   FOLLOWING.  NOTE THAT THE SUBSCRIPT I1 INCREASES
  21. C                   MOST RAPIDLY AND I3 INCREASES LEAST RAPIDLY.
  22. C                   AS OUTPUT, A CONTAINS THE COMPLEX FOURIER
  23. C                   TRANSFORM.  THE NUMBER OF CORE LOCATIONS OF
  24. C                   ARRAY A IS 2*(N1*N2*N3)
  25. C           M     - A THREE CELL VECTOR WHICH DETERMINES THE SIZES
  26. C                   OF THE 3 DIMENSIONS OF THE ARRAY A.   THE SIZE,
  27. C                   NI, OF THE I DIMENSION OF A IS 2**M(I), I = 1,2,3
  28. C           INV   - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION
  29. C                   OF DIMENSION ONE FOURTH OF THE QUANTITY
  30. C                   MAX(N1,N2,N3)
  31. C           S     - A VECTOR WORK AREA FOR SINE TABLES WITH DIMENSION
  32. C                   THE SAME AS INV
  33. C           IFSET - AN OPTION PARAMETER WITH THE FOLLOWING SETTINGS
  34. C                      0    SET UP SINE AND INV TABLES ONLY
  35. C                      1    SET UP SINE AND INV TABLES ONLY AND
  36. C                           CALCULATE FOURIER TRANSFORM
  37. C                     -1    SET UP SINE AND INV TABLES ONLY AND
  38. C                           CALCULATE INVERSE FOURIER TRANSFORM (FOR
  39. C                           THE MEANING OF INVERSE SEE THE EQUATIONS
  40. C                           UNDER METHOD BELOW)
  41. C                      2    CALCULATE FOURIER TRANSFORM ONLY (ASSUME
  42. C                           SINE AND INV TABLES EXIST)
  43. C                     -2    CALCULATE INVERSE FOURIER TRANSFORM ONLY
  44. C                           (ASSUME SINE AND INV TABLES EXIST)
  45. C           IFERR - ERROR INDICATOR.   WHEN IFSET IS 0,+1,-1,
  46. C                   IFERR = 1 MEANS THE MAXIMUM M(I) IS GREATER THAN
  47. C                  20 , I=1,2,3   WHEN IFSET IS 2,-2 , IFERR = 1
  48. C                   MEANS THAT THE SINE AND INV TABLES ARE NOT LARGE
  49. C                   ENOUGH OR HAVE NOT BEEN COMPUTED .
  50. C                   IF ON RETURN IFERR = 0 THEN NONE OF THE ABOVE
  51. C                   CONDITIONS ARE PRESENT
  52. C
  53. C        REMARKS
  54. C           THIS SUBROUTINE IS TO BE USED FOR COMPLEX, 3-DIMENSIONAL
  55. C           ARRAYS IN WHICH EACH DIMENSION IS A POWER OF 2.  THE
  56. C           MAXIMUM M(I) MUST NOT BE LESS THAN 3 OR GREATER THAN 20,
  57. C           I = 1,2,3
  58. C
  59. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  60. C           NONE
  61. C
  62. C        METHOD
  63. C           FOR IFSET = +1, OR +2, THE FOURIER TRANSFORM OF COMPLEX
  64. C           ARRAY A IS OBTAINED.
  65. C
  66. C                  N1-1   N2-1   N3-1                L1   L2   L3
  67. C     X(J1,J2,J3)=SUM    SUM    SUM    A(K1,K2,K3)*W1  *W2  *W3
  68. C                  K1=0   K2=0   K3=0
  69. C
  70. C                  WHERE WI IS THE N(I) ROOT OF UNITY AND L1=K1*J1,
  71. C                        L2=K2*J2, L3=K3*J3
  72. C
  73. C
  74. C           FOR IFSET = -1, OR -2, THE INVERSE FOURIER TRANSFORM A OF
  75. C           COMPLEX ARRAY X IS OBTAINED.
  76. C
  77. C     A(K1,K2,K3)=
  78. C               1      N1-1   N2-1   N3-1                -L1  -L2  -L3
  79. C           -------- *SUM    SUM    SUM    X(J1,J2,J3)*W1  *W2  *W3
  80. C           N1*N2*N3   J1=0   J2=0   J3=0
  81. C
  82. C
  83. C           SEE J.W. COOLEY AND J.W. TUKEY, 'AN ALGORITHM FOR THE
  84. C           MACHINE CALCULATION OF COMPLEX FOURIER SERIES',
  85. C           MATHEMATICS OF COMPUTATIONS, VOL. 19 (APR. 1965), P. 297.
  86. C
  87. C     ..................................................................
  88. C
  89.       SUBROUTINE HARM(A,M,INV,S,IFSET, IFERR)
  90.       DIMENSION A(1),INV(1),S(1),N(3),M(3),NP(3),W(2),W2(2),W3(2)
  91.       EQUIVALENCE (N1,N(1)),(N2,N(2)),(N3,N(3))
  92.    10 IF( IABS(IFSET) - 1) 900,900,12
  93.    12 MTT=MAX0(M(1),M(2),M(3)) -2
  94.       ROOT2 = SQRT(2.)
  95.       IF (MTT-MT ) 14,14,13
  96.    13 IFERR=1
  97.       RETURN
  98.    14 IFERR=0
  99.       M1=M(1)
  100.       M2=M(2)
  101.       M3=M(3)
  102.       N1=2**M1
  103.       N2=2**M2
  104.       N3=2**M3
  105.    16 IF(IFSET) 18,18,20
  106.    18 NX= N1*N2*N3
  107.       FN = NX
  108.       DO 19 I = 1,NX
  109.       A(2*I-1) = A(2*I-1)/FN
  110.    19 A(2*I) = -A(2*I)/FN
  111.    20 NP(1)=N1*2
  112.       NP(2)= NP(1)*N2
  113.       NP(3)=NP(2)*N3
  114.       DO 250 ID=1,3
  115.       IL = NP(3)-NP(ID)
  116.       IL1 = IL+1
  117.       MI = M(ID)
  118.       IF (MI)250,250,30
  119.    30 IDIF=NP(ID)
  120.       KBIT=NP(ID)
  121.       MEV = 2*(MI/2)
  122.       IF (MI - MEV )60,60,40
  123. C
  124. C     M IS ODD. DO L=1 CASE
  125.    40 KBIT=KBIT/2
  126.       KL=KBIT-2
  127.       DO 50 I=1,IL1,IDIF
  128.       KLAST=KL+I
  129.       DO 50 K=I,KLAST,2
  130.       KD=K+KBIT
  131. C
  132. C     DO ONE STEP WITH L=1,J=0
  133. C     A(K)=A(K)+A(KD)
  134. C     A(KD)=A(K)-A(KD)
  135. C
  136.       T=A(KD)
  137.       A(KD)=A(K)-T
  138.       A(K)=A(K)+T
  139.       T=A(KD+1)
  140.       A(KD+1)=A(K+1)-T
  141.    50 A(K+1)=A(K+1)+T
  142.       IF (MI - 1)250,250,52
  143.    52 LFIRST =3
  144. C
  145. C     DEF - JLAST = 2**(L-2) -1
  146.       JLAST=1
  147.       GO TO 70
  148. C
  149. C     M IS EVEN
  150.    60 LFIRST = 2
  151.       JLAST=0
  152.    70 DO 240 L=LFIRST,MI,2
  153.       JJDIF=KBIT
  154.       KBIT=KBIT/4
  155.       KL=KBIT-2
  156. C
  157. C     DO FOR J=0
  158.       DO 80 I=1,IL1,IDIF
  159.       KLAST=I+KL
  160.       DO 80 K=I,KLAST,2
  161.       K1=K+KBIT
  162.       K2=K1+KBIT
  163.       K3=K2+KBIT
  164. C
  165. C     DO TWO STEPS WITH J=0
  166. C     A(K)=A(K)+A(K2)
  167. C     A(K2)=A(K)-A(K2)
  168. C     A(K1)=A(K1)+A(K3)
  169. C     A(K3)=A(K1)-A(K3)
  170. C
  171. C     A(K)=A(K)+A(K1)
  172. C     A(K1)=A(K)-A(K1)
  173. C     A(K2)=A(K2)+A(K3)*I
  174. C     A(K3)=A(K2)-A(K3)*I
  175. C
  176.       T=A(K2)
  177.       A(K2)=A(K)-T
  178.       A(K)=A(K)+T
  179.       T=A(K2+1)
  180.       A(K2+1)=A(K+1)-T
  181.       A(K+1)=A(K+1)+T
  182. C
  183.       T=A(K3)
  184.       A(K3)=A(K1)-T
  185.       A(K1)=A(K1)+T
  186.       T=A(K3+1)
  187.       A(K3+1)=A(K1+1)-T
  188.       A(K1+1)=A(K1+1)+T
  189. C
  190.       T=A(K1)
  191.       A(K1)=A(K)-T
  192.       A(K)=A(K)+T
  193.       T=A(K1+1)
  194.       A(K1+1)=A(K+1)-T
  195.       A(K+1)=A(K+1)+T
  196. C
  197.       R=-A(K3+1)
  198.       T = A(K3)
  199.       A(K3)=A(K2)-R
  200.       A(K2)=A(K2)+R
  201.       A(K3+1)=A(K2+1)-T
  202.    80 A(K2+1)=A(K2+1)+T
  203.       IF (JLAST) 235,235,82
  204.    82 JJ=JJDIF   +1
  205. C
  206. C     DO FOR J=1
  207.       ILAST= IL +JJ
  208.       DO 85 I = JJ,ILAST,IDIF
  209.       KLAST = KL+I
  210.       DO 85 K=I,KLAST,2
  211.       K1 = K+KBIT
  212.       K2 = K1+KBIT
  213.       K3 = K2+KBIT
  214. C
  215. C     LETTING W=(1+I)/ROOT2,W3=(-1+I)/ROOT2,W2=I,
  216. C     A(K)=A(K)+A(K2)*I
  217. C     A(K2)=A(K)-A(K2)*I
  218. C     A(K1)=A(K1)*W+A(K3)*W3
  219. C     A(K3)=A(K1)*W-A(K3)*W3
  220. C
  221. C     A(K)=A(K)+A(K1)
  222. C     A(K1)=A(K)-A(K1)
  223. C     A(K2)=A(K2)+A(K3)*I
  224. C     A(K3)=A(K2)-A(K3)*I
  225. C
  226.       R =-A(K2+1)
  227.       T = A(K2)
  228.       A(K2) = A(K)-R
  229.       A(K) = A(K)+R
  230.       A(K2+1)=A(K+1)-T
  231.       A(K+1)=A(K+1)+T
  232. C
  233.       AWR=A(K1)-A(K1+1)
  234.       AWI = A(K1+1)+A(K1)
  235.       R=-A(K3)-A(K3+1)
  236.       T=A(K3)-A(K3+1)
  237.       A(K3)=(AWR-R)/ROOT2
  238.       A(K3+1)=(AWI-T)/ROOT2
  239.       A(K1)=(AWR+R)/ROOT2
  240.       A(K1+1)=(AWI+T)/ROOT2
  241.       T= A(K1)
  242.       A(K1)=A(K)-T
  243.       A(K)=A(K)+T
  244.       T=A(K1+1)
  245.       A(K1+1)=A(K+1)-T
  246.       A(K+1)=A(K+1)+T
  247.       R=-A(K3+1)
  248.       T=A(K3)
  249.       A(K3)=A(K2)-R
  250.       A(K2)=A(K2)+R
  251.       A(K3+1)=A(K2+1)-T
  252.    85 A(K2+1)=A(K2+1)+T
  253.       IF(JLAST-1) 235,235,90
  254.    90 JJ= JJ + JJDIF
  255. C
  256. C     NOW DO THE REMAINING J'S
  257.       DO 230 J=2,JLAST
  258. C
  259. C     FETCH W'S
  260. C     DEF- W=W**INV(J), W2=W**2, W3=W**3
  261.    96 I=INV(J+1)
  262.    98 IC=NT-I
  263.       W(1)=S(IC)
  264.       W(2)=S(I)
  265.       I2=2*I
  266.       I2C=NT-I2
  267.       IF(I2C)120,110,100
  268. C
  269. C     2*I IS IN FIRST QUADRANT
  270.   100 W2(1)=S(I2C)
  271.       W2(2)=S(I2)
  272.       GO TO 130
  273.   110 W2(1)=0.
  274.       W2(2)=1.
  275.       GO TO 130
  276. C
  277. C     2*I IS IN SECOND QUADRANT
  278.   120 I2CC = I2C+NT
  279.       I2C=-I2C
  280.       W2(1)=-S(I2C)
  281.       W2(2)=S(I2CC)
  282.   130 I3=I+I2
  283.       I3C=NT-I3
  284.       IF(I3C)160,150,140
  285. C
  286. C     I3 IN FIRST QUADRANT
  287.   140 W3(1)=S(I3C)
  288.       W3(2)=S(I3)
  289.       GO TO 200
  290.   150 W3(1)=0.
  291.       W3(2)=1.
  292.       GO TO 200
  293. C
  294.   160 I3CC=I3C+NT
  295.       IF(I3CC)190,180,170
  296. C
  297. C     I3 IN SECOND QUADRANT
  298.   170 I3C=-I3C
  299.       W3(1)=-S(I3C)
  300.       W3(2)=S(I3CC)
  301.       GO TO 200
  302.   180 W3(1)=-1.
  303.       W3(2)=0.
  304.       GO TO 200
  305. C
  306. C     3*I IN THIRD QUADRANT
  307.   190 I3CCC=NT+I3CC
  308.       I3CC = -I3CC
  309.       W3(1)=-S(I3CCC)
  310.       W3(2)=-S(I3CC)
  311.   200 ILAST=IL+JJ
  312.       DO 220 I=JJ,ILAST,IDIF
  313.       KLAST=KL+I
  314.       DO 220 K=I,KLAST,2
  315.       K1=K+KBIT
  316.       K2=K1+KBIT
  317.       K3=K2+KBIT
  318. C
  319. C     DO TWO STEPS WITH J NOT 0
  320. C     A(K)=A(K)+A(K2)*W2
  321. C     A(K2)=A(K)-A(K2)*W2
  322. C     A(K1)=A(K1)*W+A(K3)*W3
  323. C     A(K3)=A(K1)*W-A(K3)*W3
  324. C
  325. C     A(K)=A(K)+A(K1)
  326. C     A(K1)=A(K)-A(K1)
  327. C     A(K2)=A(K2)+A(K3)*I
  328. C     A(K3)=A(K2)-A(K3)*I
  329. C
  330.       R=A(K2)*W2(1)-A(K2+1)*W2(2)
  331.       T=A(K2)*W2(2)+A(K2+1)*W2(1)
  332.       A(K2)=A(K)-R
  333.       A(K)=A(K)+R
  334.       A(K2+1)=A(K+1)-T
  335.       A(K+1)=A(K+1)+T
  336. C
  337.       R=A(K3)*W3(1)-A(K3+1)*W3(2)
  338.       T=A(K3)*W3(2)+A(K3+1)*W3(1)
  339.       AWR=A(K1)*W(1)-A(K1+1)*W(2)
  340.       AWI=A(K1)*W(2)+A(K1+1)*W(1)
  341.       A(K3)=AWR-R
  342.       A(K3+1)=AWI-T
  343.       A(K1)=AWR+R
  344.       A(K1+1)=AWI+T
  345.       T=A(K1)
  346.       A(K1)=A(K)-T
  347.       A(K)=A(K)+T
  348.       T=A(K1+1)
  349.       A(K1+1)=A(K+1)-T
  350.       A(K+1)=A(K+1)+T
  351.       R=-A(K3+1)
  352.       T=A(K3)
  353.       A(K3)=A(K2)-R
  354.       A(K2)=A(K2)+R
  355.       A(K3+1)=A(K2+1)-T
  356.   220 A(K2+1)=A(K2+1)+T
  357. C     END OF I AND K LOOPS
  358. C
  359.   230 JJ=JJDIF+JJ
  360. C     END OF J-LOOP
  361. C
  362.   235 JLAST=4*JLAST+3
  363.   240 CONTINUE
  364. C     END OF  L  LOOP
  365. C
  366.   250 CONTINUE
  367. C     END OF  ID  LOOP
  368. C
  369. C     WE NOW HAVE THE COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE
  370. C     BIT-REVERSED.  THE FOLLOWING ROUTINE PUTS THEM IN ORDER
  371.       NTSQ=NT*NT
  372.       M3MT=M3-MT
  373.   350 IF(M3MT) 370,360,360
  374. C
  375. C     M3 GR. OR EQ. MT
  376.   360 IGO3=1
  377.       N3VNT=N3/NT
  378.       MINN3=NT
  379.       GO TO 380
  380. C
  381. C     M3 LESS THAN MT
  382.   370 IGO3=2
  383.       N3VNT=1
  384.       NTVN3=NT/N3
  385.       MINN3=N3
  386.   380 JJD3 = NTSQ/N3
  387.       M2MT=M2-MT
  388.   450 IF (M2MT)470,460,460
  389. C
  390. C     M2 GR. OR EQ. MT
  391.   460 IGO2=1
  392.       N2VNT=N2/NT
  393.       MINN2=NT
  394.       GO TO 480
  395. C
  396. C     M2 LESS THAN MT
  397.   470 IGO2 = 2
  398.       N2VNT=1
  399.       NTVN2=NT/N2
  400.       MINN2=N2
  401.   480 JJD2=NTSQ/N2
  402.       M1MT=M1-MT
  403.   550 IF(M1MT)570,560,560
  404. C
  405. C     M1 GR. OR EQ. MT
  406.   560 IGO1=1
  407.       N1VNT=N1/NT
  408.       MINN1=NT
  409.       GO TO 580
  410. C
  411. C     M1 LESS THAN MT
  412.   570 IGO1=2
  413.       N1VNT=1
  414.       NTVN1=NT/N1
  415.       MINN1=N1
  416.   580 JJD1=NTSQ/N1
  417.   600 JJ3=1
  418.       J=1
  419.       DO 880 JPP3=1,N3VNT
  420.       IPP3=INV(JJ3)
  421.       DO 870 JP3=1,MINN3
  422.       GO TO (610,620),IGO3
  423.   610 IP3=INV(JP3)*N3VNT
  424.       GO TO 630
  425.   620 IP3=INV(JP3)/NTVN3
  426.   630 I3=(IPP3+IP3)*N2
  427.   700 JJ2=1
  428.       DO 870 JPP2=1,N2VNT
  429.       IPP2=INV(JJ2)+I3
  430.       DO 860 JP2=1,MINN2
  431.       GO TO (710,720),IGO2
  432.   710 IP2=INV(JP2)*N2VNT
  433.       GO TO 730
  434.   720 IP2=INV(JP2)/NTVN2
  435.   730 I2=(IPP2+IP2)*N1
  436.   800 JJ1=1
  437.       DO 860 JPP1=1,N1VNT
  438.       IPP1=INV(JJ1)+I2
  439.       DO 850 JP1=1,MINN1
  440.       GO TO (810,820),IGO1
  441.   810 IP1=INV(JP1)*N1VNT
  442.       GO TO 830
  443.   820 IP1=INV(JP1)/NTVN1
  444.   830 I=2*(IPP1+IP1)+1
  445.       IF (J-I) 840,850,850
  446.   840 T=A(I)
  447.       A(I)=A(J)
  448.       A(J)=T
  449.       T=A(I+1)
  450.       A(I+1)=A(J+1)
  451.       A(J+1)=T
  452.   850 J=J+2
  453.   860 JJ1=JJ1+JJD1
  454. C     END OF JPP1 AND JP2
  455. C
  456.   870 JJ2=JJ2+JJD2
  457. C     END OF JPP2 AND JP3 LOOPS
  458. C
  459.   880 JJ3 = JJ3+JJD3
  460. C     END OF JPP3 LOOP
  461. C
  462.   890 IF(IFSET)891,895,895
  463.   891 DO 892 I = 1,NX
  464.   892 A(2*I) = -A(2*I)
  465.   895 RETURN
  466. C
  467. C     THE FOLLOWING PROGRAM COMPUTES THE SIN AND INV TABLES.
  468. C
  469.   900 MT=MAX0(M(1),M(2),M(3)) -2
  470.       MT = MAX0(2,MT)
  471.   904 IF (MT-18) 906,906,13
  472.   906 IFERR=0
  473.       NT=2**MT
  474.       NTV2=NT/2
  475. C
  476. C     SET UP SIN TABLE
  477. C     THETA=PIE/2**(L+1) FOR L=1
  478.   910 THETA=.7853981634
  479. C
  480. C     JSTEP=2**(MT-L+1) FOR L=1
  481.       JSTEP=NT
  482. C
  483. C     JDIF=2**(MT-L) FOR L=1
  484.       JDIF=NTV2
  485.       S(JDIF)=SIN(THETA)
  486.       DO 950 L=2,MT
  487.       THETA=THETA/2.
  488.       JSTEP2=JSTEP
  489.       JSTEP=JDIF
  490.       JDIF=JSTEP/2
  491.       S(JDIF)=SIN(THETA)
  492.       JC1=NT-JDIF
  493.       S(JC1)=COS(THETA)
  494.       JLAST=NT-JSTEP2
  495.       IF(JLAST - JSTEP) 950,920,920
  496.   920 DO 940 J=JSTEP,JLAST,JSTEP
  497.       JC=NT-J
  498.       JD=J+JDIF
  499.   940 S(JD)=S(J)*S(JC1)+S(JDIF)*S(JC)
  500.   950 CONTINUE
  501. C
  502. C     SET UP INV(J) TABLE
  503. C
  504.   960 MTLEXP=NTV2
  505. C
  506. C     MTLEXP=2**(MT-L). FOR L=1
  507.       LM1EXP=1
  508. C
  509. C     LM1EXP=2**(L-1). FOR L=1
  510.       INV(1)=0
  511.       DO 980 L=1,MT
  512.       INV(LM1EXP+1) = MTLEXP
  513.       DO 970 J=2,LM1EXP
  514.       JJ=J+LM1EXP
  515.   970 INV(JJ)=INV(J)+MTLEXP
  516.       MTLEXP=MTLEXP/2
  517.   980 LM1EXP=LM1EXP*2
  518.   982 IF(IFSET)12,895,12
  519.       END
  520.