1subroutine spec64(c0,npts,nsps,mode_q65,jpk,s3,LL,NN) 2 3 parameter (MAXFFT=20736) 4 complex c0(0:npts-1) !Complex spectrum of dd() 5 complex cs(0:MAXFFT-1) !Complex symbol spectrum 6 real s3(LL,NN) !Synchronized symbol spectra 7 real xbase0(LL),xbase(LL) 8! integer ipk1(1) 9 integer isync(22) !Indices of sync symbols 10 data isync/1,9,12,13,15,22,23,26,27,33,35,38,46,50,55,60,62,66,69,74,76,85/ 11 12 nfft=nsps 13 j=0 14 n=1 15 do k=1,84 16 if(k.eq.isync(n)) then 17 n=n+1 18 cycle 19 endif 20 j=j+1 21 ja=(k-1)*nsps + jpk 22 jb=ja+nsps-1 23 if(ja.lt.0) ja=0 24 if(jb.gt.npts-1) jb=npts-1 25 nz=jb-ja 26 cs(0:nz)=c0(ja:jb) 27 if(nz.lt.nfft-1) cs(nz+1:)=0. 28 call four2a(cs,nsps,1,-1,1) !c2c FFT to frequency 29 do ii=1,LL 30 i=ii-65+mode_q65 !mode_q65 = 1 2 4 8 16 for Q65 A B C D E 31 if(i.lt.0) i=i+nsps 32 s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2 33 enddo 34 enddo 35 36 df=6000.0/nfft 37 do i=1,LL 38 call pctile(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape 39 enddo 40 41 nh=25 42 xbase(1:nh-1)=sum(xbase0(1:nh-1))/(nh-1.0) 43 xbase(LL-nh+1:LL)=sum(xbase0(LL-nh+1:LL))/(nh-1.0) 44 do i=nh,LL-nh 45 xbase(i)=sum(xbase0(i-nh+1:i+nh))/(2*nh+1) !Smoothed passband shape 46 enddo 47 48 do i=1,LL 49 s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization 50 enddo 51 52 return 53end subroutine spec64 54