C 999.FOR FOR DSP COMPLEX X(1024), Y(1024), A(1024), B(1024),W(1024) REAL P, PB, XX(1024), AIC(50),MDL(50),PSD(4096) INTEGER XXX(1024) OPEN(UNIT=8, FILE='929416.DAT', TYPE='OLD') 11 FORMAT(1E26.17,1F26.17) 22 FORMAT(I4,2E14.6) 33 FORMAT(I4,1E14.6) 44 FORMAT(2F14.3) N=1024 T=1.0/1000000.0 DO 10 J=1,N READ(8,*) YY, XX(J) XR=XX(J) XI=0 X(J)=CMPLX(XR,XI) 10 CONTINUE CLOSE(UNIT=8) WRITE(*, *) XX(10) CALL PREFFT(N,1,NEXP,W) CALL FFT(N,1,T,NEXP,W,X) NSIG=3 CALL EIGENF(1,N,5,NSIG,X,PSD,ISTAT) DO 20 I=1,NSIG WRITE(*,*) PSD(I),ISTAT 20 CONTINUE STOP END C LIST OF SUBROUTINE C These two programs set up the complex exponential table (PREFFT) and C compute the discrete-time Fourier series of an array of complex data C samples using a decimation-in-frequency fast Fourier transform (FFT) C algorithm. C SUBROUTINE PREFFT (N,MODE,NEXP,W) C C Input Parameters: C C N - Number of data samples to be processed (integer-must be a C power of two) C MODE - Set to 0 for discrete-time Fourier series (Eq. 2.C.1) or C 1 for inverse (Eq. 2.C.2) C C Output Parameters: C C NEXP - Indicates power-of-2 exponent such that N=2**NEXP . C Will be set to -1 to indicate error condition if N C is not a power of 2 (this integer used by sub. FFT) C W - Complex exponential array C C Notes: C C External array W must be dimensioned .GE. N by calling program. C COMPLEX W(1024),C1,C2 NEXP=1 5 NT=2**NEXP IF (NT .GE. N) GO TO 10 NEXP=NEXP+1 GO TO 5 10 IF (NT .EQ. N) GO TO 15 NEXP=-1 RETURN 15 S=8.*ATAN(1.)/FLOAT(NT) C1=CMPLX(COS(S),-SIN(S)) IF (MODE .NE. 0) C1=CONJG(C1) C2=(1.,0.) DO 20 K=1,NT W(K)=C2 20 C2=C2*C1 RETURN END SUBROUTINE FFT (N,MODE,T,NEXP,W,X) C C Input Parameters: C C N,MODE,NEXP,W - See parameter list for subroutine PREFFT C T - Sample interval in seconds C X - Array of N complex data samples, X(1) to X(N) C C Output Parameters: C C X - N complex transform values replace original data samples C indexed from k=1 to k=N, representing the frequencies C (k-1)/NT hertz C C Notes: C C External array X must be dimensioned .GE. N by calling program. C COMPLEX X(1024),W(1024),C1,C2 MM=1 LL=N DO 70 K=1,NEXP NN=LL/2 JJ=MM+1 DO 40 I=1,N,LL KK=I+NN C1=X(I)+X(KK) X(KK)=X(I)-X(KK) 40 X(I)=C1 IF (NN .EQ. 1) GO TO 70 DO 60 J=2,NN C2=W(JJ) DO 50 I=J,N,LL KK=I+NN C1=X(I)+X(KK) X(KK)=(X(I)-X(KK))*C2 50 X(I)=C1 60 JJ=JJ+MM LL=NN MM=MM*2 70 CONTINUE NV2=N/2 NM1=N-1 J=1 DO 90 I=1,NM1 IF (I .GE. J) GO TO 80 C1=X(J) X(J)=X(I) X(I)=C1 80 K=NV2 85 IF (K .GE. J) GO TO 90 J=J-K K=K/2 GO TO 85 90 J=J+K IF (MODE .EQ. 0) S=T IF (MODE .NE. 0) S=1./(T*FLOAT(N)) DO 100 I=1,N 100 X(I)=X(I)*S RETURN END SUBROUTINE EIGENF (METHOD,N,IP,NSIG,X,PSD,ISTAT) C C Input Parameters: C C METHOD - Set to 0 for MUSIC method, 1 for EV method (integer) C N - Number of data samples (integer) C IP - Dimension of data matrix (integer) C NSIG - Assumed number of "signal" space vectors (integer) C X - Array of complex data samples C C Output Parameters: C C PSD - Array of real frequency estimator values C ISTAT - Status indicator (integer): C 0 for normal exit C 1 if NSIG > IP-1 or NSIG < 0 C C Notes: C C External array X must be dimensioned .GE. N and array PSD must be C dimensioned .GE. NPSD in calling program. Internal array FB must C be dimensioned .GE. 2*(N-IP) x IP , array U must be dimensioned C .GE. 2*(N-IP) x 2*(N-IP), array V must be dimensioned .GE. IP x IP, C and array S must be dimensioned .GE. IP or 2*(N-IP) , whichever is C greater. C PARAMETER (NPSD=4096) PARAMETER (MAXU=100) PARAMETER (MAXV=50) COMPLEX X(1024),FB(MAXU,MAXV),U(MAXU,MAXU),V(MAXV,MAXV) COMPLEX W(NPSD),Z(NPSD) REAL PSD(1),S(MAXV) PI2=8.*ATAN(1.) ISTAT=0 IF (NSIG .LT. IP .AND. NSIG .GE. 0) GO TO 5 ISTAT=1 RETURN 5 NP=N-IP DO 20 I=1,NP DO 10 K=1,IP FB(I,K)=X(I-K+IP) 10 FB(I+NP,K)=CONJG(X(I+K)) 20 CONTINUE NP2=2*NP CALL CSVD (FB,MAXU,MAXV,NP2,IP,0,IP,IP,S,U,V) C AI or Expert Knowledge to choose "signal" singular values, or input C NSIG at this point CALL PREFFT (NPSD,0,NEXP,W) DO 30 K=1,NPSD 30 PSD(K)=0. DO 70 I=NSIG+1,IP DO 40 K=1,IP 40 Z(K)=V(K,I) DO 50 K=IP+1,NPSD 50 Z(K)=(0.,0.) CALL FFT (NPSD,0,1.,NEXP,W,Z) DO 60 K=1,NPSD TEMP=REAL(Z(K))**2+AIMAG(Z(K))**2 IF (METHOD .EQ. 0) PSD(K)=PSD(K)+TEMP 60 IF (METHOD .NE. 0) PSD(K)=PSD(K)+TEMP/S(I) 70 CONTINUE DO 80 K=1,NPSD 80 PSD(K)=1./PSD(K) RETURN END SUBROUTINE CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V) C C Singular value decomposition of an M by N complex matrix A, C where M .GT. N . The singular values are stored in the vector C S. The first NU columns of the M by M unitary matrix U and the C first NV columns of the N by N unitary matrix V that minimize C Det(A-USV*) are also computed. C C C P.A. Businger and G.H. Golub, "Singular Value Decomposition C of a Complex Matrix," Communications of the ACM, vol. 12, C pp. 564-565, October 1969. C C This algorithm is reprinted by permission, Association for C Computing Machinery; copyright 1969. C COMPLEX A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R REAL S(NMAX),B(100),C(100),T(100) DATA ETA,TOL/1.2E-7,2.4E-32/ NP=N+IP N1=N+1 C Householder reduction C(1)=0. K=1 10 K1=K+1 C Elimination of A(I,K) , I=K+1,...,M Z=0. DO 20 I=K,M 20 Z=Z+REAL(A(I,K))**2+AIMAG(A(I,K))**2 B(K)=0. IF (Z .LE. TOL) GO TO 70 Z=SQRT(Z) B(K)=Z W=CABS(A(K,K)) Q=(1.,0.) IF (W .NE. 0.) Q=A(K,K)/W A(K,K)=Q*(Z+W) IF (K .EQ. NP) GO TO 70 DO 50 J=K1,NP Q=(0.,0.) DO 30 I=K,M 30 Q=Q+CONJG(A(I,K))*A(I,J) Q=Q/(Z*(Z+W)) DO 40 I=K,M 40 A(I,J)=A(I,J)-Q*A(I,K) 50 CONTINUE C Phase transformation Q=-CONJG(A(K,K))/CABS(A(K,K)) DO 60 J=K1,NP 60 A(K,J)=Q*A(K,J) C Elimination of A(K,J) , J=K+2,...,N 70 IF (K .EQ. N) GO TO 140 Z=0. DO 80 J=K1,N 80 Z=Z+REAL(A(K,J))**2+AIMAG(A(K,J))**2 C(K1)=0. IF (Z .LE. TOL) GO TO 130 Z=SQRT(Z) C(K1)=Z W=CABS(A(K,K1)) Q=(1.,0.) IF (W .NE. 0.) Q=A(K,K1)/W A(K,K1)=Q*(Z+W) DO 110 I=K1,M Q=(0.,0.) DO 90 J=K1,N 90 Q=Q+CONJG(A(K,J))*A(I,J) Q=Q/(Z*(Z+W)) DO 100 J=K1,N 100 A(I,J)=A(I,J)-Q*A(K,J) 110 CONTINUE C Phase transformation Q=-CONJG(A(K,K1))/CABS(A(K,K1)) DO 120 I=K1,M 120 A(I,K1)=A(I,K1)*Q 130 K=K1 GO TO 10 C Tolerance for negligible elements 140 EPS=0. DO 150 K=1,N S(K)=B(K) T(K)=C(K) 150 EPS=AMAX1(EPS,S(K)+T(K)) EPS=EPS*ETA C Initialization of U and V IF (NU .EQ. 0) GO TO 180 DO 170 J=1,NU DO 160 I=1,M 160 U(I,J)=(0.,0.) 170 U(J,J)=(1.,0.) 180 IF (NV .EQ. 0) GO TO 210 DO 200 J=1,NV DO 190 I=1,N 190 V(I,J)=(0.,0.) 200 V(J,J)=(1.,0.) C QR diagonalization 210 DO 380 KK=1,N K=N1-KK C Test for split 220 DO 230 LL=1,K L=K+1-LL IF (ABS(T(L)) .LE. EPS) GO TO 290 IF (ABS(S(L-1)) .LE. EPS) GO TO 240 230 CONTINUE C Cancellation of B(L) 240 CS=0. SN=1. L1=L-1 DO 280 I=L,K F=SN*T(I) T(I)=CS*T(I) IF (ABS(F) .LE. EPS) GO TO 290 H=S(I) W=SQRT(F*F+H*H) S(I)=W CS=H/W SN=-F/W IF (NU .EQ. 0) GO TO 260 DO 250 J=1,N X=REAL(U(J,L1)) Y=REAL(U(J,I)) U(J,L1)=CMPLX(X*CS+Y*SN,0.) 250 U(J,I)=CMPLX(Y*CS-X*SN,0.) 260 IF (NP .EQ. N) GO TO 280 DO 270 J=N1,NP Q=A(L1,J) R=A(I,J) A(L1,J)=Q*CS+R*SN 270 A(I,J)=R*CS-Q*SN 280 CONTINUE C Test for convergence 290 W=S(K) IF (L .EQ. K) GO TO 360 C Origin shift X=S(L) Y=S(K-1) G=T(K-1) H=T(K) F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.*H*Y) G=SQRT(F*F+1.) IF (F .LT. 0.) G=-G F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X C QR step CS=1. SN=1. L1=L+1 DO 350 I=L1,K G=T(I) Y=S(I) H=SN*G G=CS*G W=SQRT(H*H+F*F) T(I-1)=W CS=F/W SN=H/W F=X*CS+G*SN G=G*CS-X*SN H=Y*SN Y=Y*CS IF (NV .EQ. 0) GO TO 310 DO 300 J=1,N X=REAL(V(J,I-1)) W=REAL(V(J,I)) V(J,I-1)=CMPLX(X*CS+W*SN,0.) 300 V(J,I)=CMPLX(W*CS-X*SN,0.) 310 W=SQRT(H*H+F*F) S(I-1)=W CS=F/W SN=H/W F=CS*G+SN*Y X=CS*Y-SN*G IF (NU .EQ. 0) GO TO 330 DO 320 J=1,N Y=REAL(U(J,I-1)) W=REAL(U(J,I)) U(J,I-1)=CMPLX(Y*CS+W*SN,0.) 320 U(J,I)=CMPLX(W*CS-Y*SN,0.) 330 IF (N .EQ. NP) GO TO 350 DO 340 J=N1,NP Q=A(I-1,J) R=A(I,J) A(I-1,J)=Q*CS+R*SN 340 A(I,J)=R*CS-Q*SN 350 CONTINUE T(L)=0. T(K)=F S(K)=X GO TO 220 C Convergence 360 IF (W .GE. 0.) GO TO 380 S(K)=-W IF (NV .EQ. 0) GO TO 380 DO 370 J=1,N 370 V(J,K)=-V(J,K) 380 CONTINUE C Sort singular values DO 450 K=1,N G=-1. J=K DO 390 I=K,N IF (S(I) .LE. G) GO TO 390 G=S(I) J=I 390 CONTINUE IF (J .EQ. K) GO TO 450 S(J)=S(K) S(K)=G IF (NV .EQ. 0) GO TO 410 DO 400 I=1,N Q=V(I,J) V(I,J)=V(I,K) 400 V(I,K)=Q 410 IF (NU .EQ. 0) GO TO 430 DO 420 I=1,N Q=U(I,J) U(I,J)=U(I,K) 420 U(I,K)=Q 430 IF (N .EQ. NP) GO TO 450 DO 440 I=N1,NP Q=A(J,I) A(J,I)=A(K,I) 440 A(K,I)=Q 450 CONTINUE C Back transformation IF (NU .EQ. 0) GO TO 510 DO 500 KK=1,N K=N1-KK IF (B(K) .EQ. 0.) GO TO 500 Q=-A(K,K)/CABS(A(K,K)) DO 460 J=1,NU 460 U(K,J)=Q*U(K,J) DO 490 J=1,NU Q=(0.,0.) DO 470 I=K,M 470 Q=Q+CONJG(A(I,K))*U(I,J) Q=Q/(CABS(A(K,K))*B(K)) DO 480 I=K,M 480 U(I,J)=U(I,J)-Q*A(I,K) 490 CONTINUE 500 CONTINUE 510 IF (NV .EQ. 0) GO TO 570 IF (N .LT. 2) GO TO 570 DO 560 KK=2,N K=N1-KK K1=K+1 IF (C(K1) .EQ. 0.) GO TO 560 Q=-CONJG(A(K,K1))/CABS(A(K,K1)) DO 520 J=1,NV 520 V(K1,J)=Q*V(K1,J) DO 550 J=1,NV Q=(0.,0.) DO 530 I=K1,N 530 Q=Q+A(K,I)*V(I,J) Q=Q/(CABS(A(K,K1))*C(K1)) DO 540 I=K1,N 540 V(I,J)=V(I,J)-Q*CONJG(A(K,I)) 550 CONTINUE 560 CONTINUE 570 RETURN END