diff --git a/Makefile.win b/Makefile.win index 0206ddac0..980ac0286 100644 --- a/Makefile.win +++ b/Makefile.win @@ -36,7 +36,7 @@ SRCS2F77 = wsjt1.f avesp2.f bzap.f spec441.f spec2d.f mtdecode.f \ limit.f lpf1.f deep65.f morse.f nchar.f packcall.f packgrid.f \ packmsg.f packtext.f setup65.f short65.f slope.f spec2d65.f \ sync65.f unpackcall.f unpackgrid.f unpackmsg.f unpacktext.f \ - xcor.f xfft.f wsjt65.f astro.f azdist.f coord.f dcoord.f \ + xcor.f xfft.f xfft2.f wsjt65.f astro.f azdist.f coord.f dcoord.f \ deg2grid.f dot.f ftsky.f geocentric.f GeoDist.f grid2deg.f \ moon2.f MoonDop.f sun.f toxyz.f pfxdump.f diff --git a/horizspec.f90 b/horizspec.f90 index e551f2a18..0ff302411 100644 --- a/horizspec.f90 +++ b/horizspec.f90 @@ -36,7 +36,7 @@ subroutine horizspec(x,brightness,contrast,a) do i=1,nfft y(i)=1.4*x(i+i0) enddo - call xfft(y,nfft) + call xfft2(y,nfft) nq=nfft/4 do i=1,nq ss(i)=real(c(i))**2 + imag(c(i))**2 diff --git a/spec.f90 b/spec.f90 index 0b24027e5..65e7a7e60 100644 --- a/spec.f90 +++ b/spec.f90 @@ -36,7 +36,9 @@ subroutine spec(brightness,contrast,logmap,ngain,nspeed,a) save if(first) then - call zero(ss,nq) + do i=1,nq + ss(i)=0. + enddo istep=2205 nfft=4096 nq=nfft/4 @@ -136,7 +138,7 @@ subroutine spec(brightness,contrast,logmap,ngain,nspeed,a) go to 900 endif - call xfft(x,nfft) + call xfft2(x,nfft) do i=1,nq !Accumulate power spectrum ss(i)=ss(i) + real(c(i))**2 + imag(c(i))**2 @@ -167,7 +169,9 @@ subroutine spec(brightness,contrast,logmap,ngain,nspeed,a) enddo nsum=0 newdat=1 !Flag for new spectrum available - call zero(ss,nq) !Zero the accumulating array + do i=1,nq !Zero the accumulating array + ss(i)=0. + enddo if(jz.lt.300) jz=jz+1 endif diff --git a/wsjt.py b/wsjt.py index 9f44e186c..7d2eb317d 100644 --- a/wsjt.py +++ b/wsjt.py @@ -1,4 +1,4 @@ -#---------------------------------------------------------------- WSJT +#----------------------------------------------------------------- WSJT from Tkinter import * from tkFileDialog import * import Pmw diff --git a/xfft2.f b/xfft2.f new file mode 100644 index 000000000..0eae674a5 --- /dev/null +++ b/xfft2.f @@ -0,0 +1,184 @@ + SUBROUTINE xfft2(DATA,NB) +c +c the cooley-tukey fast fourier transform in usasi basic fortran +c +C .. Scalar Arguments .. + INTEGER NB +C .. +C .. Array Arguments .. + REAL DATA(NB+2) +C .. +C .. Local Scalars .. + REAL DIFI,DIFR,RTHLF,SUMI,SUMR,T2I,T2R,T3I,T3R,T4I, + + T4R,TEMPI,TEMPR,THETA,TWOPI,U1I,U1R,U2I,U2R,U3I,U3R, + + U4I,U4R,W2I,W2R,W3I,W3R,WI,WR,WSTPI,WSTPR + INTEGER I,I2,IPAR,J,K1,K2,K3,K4,KDIF,KMIN, + + KSTEP,L,LMAX,M,MMAX,NH +C .. +C .. Intrinsic Functions .. + INTRINSIC COS,MAX0,REAL,SIN +C .. +C .. Data statements .. + DATA TWOPI/6.2831853071796/,RTHLF/0.70710678118655/ +c +c 1. real transform for the 1st dimension, n even. method-- +c transform a complex array of length n/2 whose real parts +c are the even numbered real values and whose imaginary parts +c are the odd numbered real values. separate and supply +c the second half by conjugate symmetry. +c + + NH = NB/2 +c +c shuffle data by bit reversal, since n=2**k. +c + J = 1 + DO 131 I2 = 1,NB,2 + IF (J-I2) 124,127,127 + 124 TEMPR = DATA(I2) + TEMPI = DATA(I2+1) + DATA(I2) = DATA(J) + DATA(I2+1) = DATA(J+1) + DATA(J) = TEMPR + DATA(J+1) = TEMPI + 127 M = NH + 128 IF (J-M) 130,130,129 + 129 J = J - M + M = M/2 + IF (M-2) 130,128,128 + 130 J = J + M + 131 CONTINUE + +c +c main loop for factors of two. perform fourier transforms of +c length four, with one of length two if needed. the twiddle factor +c w=exp(-2*pi*sqrt(-1)*m/(4*mmax)). check for w=-sqrt(-1) +c and repeat for w=w*(1-sqrt(-1))/sqrt(2). +c + IF (NB-2) 174,174,143 + 143 IPAR = NH + 144 IF (IPAR-2) 149,146,145 + 145 IPAR = IPAR/4 + GO TO 144 + + 146 DO 147 K1 = 1,NB,4 + K2 = K1 + 2 + TEMPR = DATA(K2) + TEMPI = DATA(K2+1) + DATA(K2) = DATA(K1) - TEMPR + DATA(K2+1) = DATA(K1+1) - TEMPI + DATA(K1) = DATA(K1) + TEMPR + DATA(K1+1) = DATA(K1+1) + TEMPI + 147 CONTINUE + 149 MMAX = 2 + 150 IF (MMAX-NH) 151,174,174 + 151 LMAX = MAX0(4,MMAX/2) + DO 173 L = 2,LMAX,4 + M = L + IF (MMAX-2) 156,156,152 + 152 THETA = -TWOPI*REAL(L)/REAL(4*MMAX) + WR = COS(THETA) + WI = SIN(THETA) + 155 W2R = WR*WR - WI*WI + W2I = 2.*WR*WI + W3R = W2R*WR - W2I*WI + W3I = W2R*WI + W2I*WR + 156 KMIN = 1 + IPAR*M + IF (MMAX-2) 157,157,158 + 157 KMIN = 1 + 158 KDIF = IPAR*MMAX + 159 KSTEP = 4*KDIF + IF (KSTEP-NB) 160,160,169 + 160 DO 168 K1 = KMIN,NB,KSTEP + K2 = K1 + KDIF + K3 = K2 + KDIF + K4 = K3 + KDIF + IF (MMAX-2) 161,161,164 + 161 U1R = DATA(K1) + DATA(K2) + U1I = DATA(K1+1) + DATA(K2+1) + U2R = DATA(K3) + DATA(K4) + U2I = DATA(K3+1) + DATA(K4+1) + U3R = DATA(K1) - DATA(K2) + U3I = DATA(K1+1) - DATA(K2+1) + U4R = DATA(K3+1) - DATA(K4+1) + U4I = DATA(K4) - DATA(K3) + GO TO 167 + + 164 T2R = W2R*DATA(K2) - W2I*DATA(K2+1) + T2I = W2R*DATA(K2+1) + W2I*DATA(K2) + T3R = WR*DATA(K3) - WI*DATA(K3+1) + T3I = WR*DATA(K3+1) + WI*DATA(K3) + T4R = W3R*DATA(K4) - W3I*DATA(K4+1) + T4I = W3R*DATA(K4+1) + W3I*DATA(K4) + U1R = DATA(K1) + T2R + U1I = DATA(K1+1) + T2I + U2R = T3R + T4R + U2I = T3I + T4I + U3R = DATA(K1) - T2R + U3I = DATA(K1+1) - T2I + U4R = T3I - T4I + U4I = T4R - T3R + + 167 DATA(K1) = U1R + U2R + DATA(K1+1) = U1I + U2I + DATA(K2) = U3R + U4R + DATA(K2+1) = U3I + U4I + DATA(K3) = U1R - U2R + DATA(K3+1) = U1I - U2I + DATA(K4) = U3R - U4R + DATA(K4+1) = U3I - U4I + 168 CONTINUE + KDIF = KSTEP + KMIN = 4*KMIN - 3 + GO TO 159 + + 169 M = M + LMAX + IF (M-MMAX) 170,170,173 + 170 TEMPR = WR + WR = (WR+WI)*RTHLF + WI = (WI-TEMPR)*RTHLF + GO TO 155 + + 173 CONTINUE + IPAR = 3 - IPAR + MMAX = MMAX + MMAX + GO TO 150 +c +c complete a real transform in the 1st dimension, n even, by con- +c jugate symmetries. +c + 174 THETA = -TWOPI/REAL(NB) + WSTPR = COS(THETA) + WSTPI = SIN(THETA) + WR = WSTPR + WI = WSTPI + I = 3 + J = NB - 1 + GO TO 207 + + 205 SUMR = (DATA(I)+DATA(J))/2. + SUMI = (DATA(I+1)+DATA(J+1))/2. + DIFR = (DATA(I)-DATA(J))/2. + DIFI = (DATA(I+1)-DATA(J+1))/2. + TEMPR = WR*SUMI + WI*DIFR + TEMPI = WI*SUMI - WR*DIFR + DATA(I) = SUMR + TEMPR + DATA(I+1) = DIFI + TEMPI + DATA(J) = SUMR - TEMPR + DATA(J+1) = -DIFI + TEMPI + I = I + 2 + J = J - 2 + TEMPR = WR + WR = WR*WSTPR - WI*WSTPI + WI = TEMPR*WSTPI + WI*WSTPR + 207 IF (I-J) 205,208,211 + 208 DATA(I+1) = -DATA(I+1) + + 211 DATA(NB+1) = DATA(1) - DATA(2) + DATA(NB+2) = 0. + + DATA(1) = DATA(1) + DATA(2) + DATA(2) = 0. + + RETURN + END