1subroutine analytic(d,npts,nfft,c,pc,beq) 2 3! Convert real data to analytic signal 4 5 parameter (NFFTMAX=1024*1024) 6 7 real d(npts) ! passband signal 8 real h(NFFTMAX/2) ! real BPF magnitude 9 real*8 pc(5),pclast(5) ! static phase coeffs 10 real*8 ac(5),aclast(5) ! amp coeffs 11 real*8 fp 12 13 complex corr(NFFTMAX/2) ! complex frequency-dependent correction 14 complex c(NFFTMAX) ! analytic signal 15 16 logical*1 beq ! boolean static equalizer flag 17 18 data nfft0/0/ 19 data aclast/0.0,0.0,0.0,0.0,0.0/ 20 data pclast/0.0,0.0,0.0,0.0,0.0/ 21! data ac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000 22 data ac/1.0,0.0,0.0,0.0,0.0/ 23 24 save corr,nfft0,h,ac,aclast,pclast,pi,t,beta 25 26 df=12000.0/nfft 27 nh=nfft/2 28 if( nfft.ne.nfft0 ) then 29 pi=4.0*atan(1.0) 30 t=1.0/2000.0 31 beta=0.1 32 do i=1,nh+1 33 ff=(i-1)*df 34 f=ff-1500.0 35 h(i)=1.0 36 if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then 37 h(i)=h(i)*0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t)))) 38 elseif( abs(f) .gt. (1+beta)/(2*t) ) then 39 h(i)=0.0 40 endif 41 enddo 42 nfft0=nfft 43 endif 44 45 if( any(aclast .ne. ac) .or. any(pclast .ne. pc) ) then 46 aclast=ac 47 pclast=pc 48! write(*,3001) pc 49!3001 format('Phase coeffs:',5f12.6) 50 do i=1,nh+1 51 ff=(i-1)*df 52 f=ff-1500.0 53 fp=f/1000.0 54 corr(i)=ac(1)+fp*(ac(2)+fp*(ac(3)+fp*(ac(4)+fp*ac(5)))) 55 pd=fp*fp*(pc(3)+fp*(pc(4)+fp*pc(5))) ! ignore 1st two terms 56 corr(i)=corr(i)*cmplx(cos(pd),sin(pd)) 57 enddo 58 endif 59 60 fac=2.0/nfft 61 c(1:npts)=fac*d(1:npts) 62 c(npts+1:nfft)=0. 63 call four2a(c,nfft,1,-1,1) !Forward c2c FFT 64 65 if( beq ) then 66 c(1:nh+1)=h(1:nh+1)*corr(1:nh+1)*c(1:nh+1) 67 else 68 c(1:nh+1)=h(1:nh+1)*c(1:nh+1) 69 endif 70 71 c(1)=0.5*c(1) !Half of DC term 72 c(nh+2:nfft)=0. !Zero the negative frequencies 73 call four2a(c,nfft,1,1,1) !Inverse c2c FFT 74 return 75end subroutine analytic 76