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