subroutine timf2(x0,k,nfft,nwindow,nb,peaklimit,faclim,x1, & slimit,lstrong,px,nzap) ! Sequential processing of time-domain I/Q data, using Linrad-like ! "first FFT" and "first backward FFT", treating frequencies with ! strong signals differently. Noise blanking is applied to weak ! signals only. ! x0 - real input data ! nfft - length of FFTs ! nwindow - 0 for no window, 2 for sin^2 window ! x1 - real output data ! Non-windowed processing means no overlap, so kstep=nfft. ! Sin^2 window has 50% overlap, kstep=nfft/2. ! Frequencies with strong signals are identified and separated. Back ! transforms are done separately for weak and strong signals, so that ! noise blanking can be applied to the weak-signal portion. Strong and ! weak are finally re-combined, in the time domain. parameter (MAXFFT=1024,MAXNH=MAXFFT/2) parameter (MAXSIGS=100) real x0(0:nfft-1),x1(0:nfft-1) real x(0:MAXFFT-1),xw(0:MAXFFT-1),xs(0:MAXFFT-1) real xwov(0:MAXNH-1),xsov(0:MAXNH-1) complex cx(0:MAXFFT-1),cxt(0:MAXFFT-1) complex cxs(0:MAXFFT-1) !Strong signals complex cxw(0:MAXFFT-1) !Weak signals real*4 w(0:MAXFFT-1) real*4 s(0:MAXNH),stmp(0:MAXNH) logical*1 lstrong(0:MAXNH),lprev integer ia(MAXSIGS),ib(MAXSIGS) logical first equivalence (x,cx),(xw,cxw),(xs,cxs) data first/.true./ data k0/99999999/ ! save w,xsov,xwov,s,ntc,ntot,nh,kstep,fac,first,k0 save if(first) then pi=4.0*atan(1.0) do i=0,nfft-1 w(i)=(sin(i*pi/nfft))**2 enddo s=0. ntc=0 ntot=0 nh=nfft/2 kstep=nfft if(nwindow.eq.2) kstep=nh fac=1.0/nfft slimit=1.e30 first=.false. endif if(k.lt.k0) then xsov=0. xwov=0. endif k0=k x(0:nfft-1)=x0 if(nwindow.eq.2) x(0:nfft-1)=w(0:nfft-1)*x(0:nfft-1) call four2a(x,nfft,1,-1,0) !First forward FFT, r2c cxt(0:nh)=cx(0:nh) ! Identify frequencies with strong signals. !### do i=0,nh p=real(cxt(i))**2 + aimag(cxt(i))**2 s(i)=p enddo ave=sum(s(0:nh))/nh lstrong(0:nh)=s(0:nh).gt.10.0*ave !### ! ntot=ntot+1 ! if(mod(ntot,128).eq.5) then ! call pctile(s,stmp,nh,50,xmedian) ! slimit=faclim*xmedian ! endif ! if(ntc.lt.12000/nfft) ntc=ntc+1 ! uu=1.0/ntc ! smax=0. ! do i=0,nh ! p=real(cxt(i))**2 + aimag(cxt(i))**2 ! s(i)=(1.0-uu)*s(i) + uu*p ! lstrong(i)=(s(i).gt.slimit) ! if(s(i).gt.smax) smax=s(i) ! enddo nsigs=0 lprev=.false. iwid=1 ib=-99 do i=0,nh if(lstrong(i) .and. (.not.lprev)) then if(nsigs.lt.MAXSIGS) nsigs=nsigs+1 ia(nsigs)=i-iwid if(ia(nsigs).lt.0) ia(nsigs)=0 endif if(.not.lstrong(i) .and. lprev) then ib(nsigs)=i-1+iwid if(ib(nsigs).gt.nh) ib(nsigs)=nh endif lprev=lstrong(i) enddo if(nsigs.gt.0) then do i=1,nsigs ja=ia(i) jb=ib(i) if(ja.lt.0 .or. ja.gt.nh .or. jb.lt.0 .or. jb.gt.nh) then cycle endif if(jb.eq.-99) jb=ja + min(2*iwid,nh) lstrong(ja:jb)=.true. enddo endif ! Copy frequency-domain data into array cs (strong) or cw (weak). do i=0,nh if(lstrong(i)) then cxs(i)=fac*cxt(i) cxw(i)=0. else cxw(i)=fac*cxt(i) cxs(i)=0. endif enddo call four2a(cxw,nfft,1,1,-1) !Transform weak and strong back call four2a(cxs,nfft,1,1,-1) !to time domain, separately (c2r) if(nwindow.eq.2) then xw(0:nh-1)=xw(0:nh-1)+xwov(0:nh-1) !Add previous segment's 2nd half xwov(0:nh-1)=xw(nh:nfft-1) !Save 2nd half xs(0:nh-1)=xs(0:nh-1)+xsov(0:nh-1) !Ditto for strong signals xsov(0:nh-1)=xs(nh:nfft-1) endif ! Apply noise blanking to weak data if(nb.ne.0) then do i=0,kstep-1 peak=abs(xw(i)) if(peak.gt.peaklimit) then xw(i)=0. nzap=nzap+1 endif enddo endif ! Compute power levels from weak data only do i=0,kstep-1 px=px + xw(i)**2 enddo x1(0:kstep-1)=xw(0:kstep-1) + xs(0:kstep-1) !Recombine weak + strong return end subroutine timf2