WSJT-X/spec2d.f

	subroutine spec2d(data,jz,nstep,s2,nchan,nz,psavg0,sigma)

C  Computes 2d spectrogram for FSK441 single-tone search and waterfall
C  display.

	parameter (NFFT=256)
	parameter (NR=NFFT+2)
	parameter (NH=NFFT/2)
	parameter (NQ=NFFT/4)

	real data(jz)
	real s2(nchan,nz)
	real x(NR)
	real w1(7),w2(7)
	real psavg(128)
	real psavg0(128)
	real ps2(128)
	complex c(0:NH)
	common/acom/a1,a2,a3,a4
	common/fcom/s(3100),indx(3100)
	equivalence (x,c)
	save

	df=11025.0/NFFT

C  Compute the 2d spectrogram s2(nchan,nz).  Note that in s2 the frequency
C  bins are shifted down 5 bins from their natural positions.

	call set(0.0,psavg,NH)
	do n=1,nz
	   j=1 + (n-1)*nstep
	   call move(data(j),x,NFFT)
	   call xfft(x,NFFT)

	   sum=0.
	   do i=1,NQ
	      s2(i,n)=real(c(5+i))**2 + imag(c(5+i))**2
	      sum=sum+s2(i,n)
	   enddo
	   s(n)=sum/NQ

C  Accumulate average spectrum for the whole file.
	   do i=1,nh
	      psavg0(i) = psavg0(i)+ real(c(i))**2 + imag(c(i))**2
	   enddo
	enddo

C  Normalize and save a copy of psavg0 for plotting.  Roll off the
C  spectrum at 300 and 3000 Hz.
	do i=1,nh
	   psavg0(i)=3.e-5*psavg0(i)/nz
	   f=df*i
	   fac=1.0
	   if(f.lt.300.0) fac=f/300.0
	   if(f.gt.3000.0) fac=max(0.00333,(3300.0-f)/300.0)
	   psavg0(i)=(fac**2)*psavg0(i)
	enddo

C  Compute an average spectrum from the weakest 25% of time slices.
	call indexx(nz,s,indx)
	call zero(ps2,NQ)
	do j=1,nz/4
	   k=indx(j)
	   do i=1,NQ
	      ps2(i+5)=ps2(i+5)+s2(i,k)
	   enddo
	enddo
	ps2(1)=ps2(5)
	ps2(2)=ps2(5)
	ps2(3)=ps2(5)
	ps2(4)=ps2(5)
	sum=0.
	do i=6,59
	   sum=sum+ps2(i)
	enddo
	if(sum.eq.0.0) then
	   sigma=-999.
	   go to 999
	endif

C  Compute a smoothed spectrum without local peaks, and find its max.
	smaxx=0.
	do i=4,NQ
	   sum=0.
	   do k=1,7
	      w1(k)=ps2(i+k-4)
	      sum=sum+w1(k)
	   enddo
	   ave=sum/7.0
	   if(i.ge.14 .and. i.le.58) then
	      call pctile(w1,w2,7,50,base)
	      ave=0.25*(w2(1)+w2(2)+w2(3)+w2(4))
	   endif
	   psavg(i)=ave
	   smaxx=max(psavg(i),smaxx)
	enddo

C  Save scale factors for flattening spectra of pings.
	a1=1.0
	a2=psavg(nint(2*441/df))/psavg(nint(3*441/df))
	a3=psavg(nint(2*441/df))/psavg(nint(4*441/df))
	a4=psavg(nint(2*441/df))/psavg(nint(5*441/df))
	afac=4.0/(a1+a2+a3+a4)
	a1=afac*a1
	a2=afac*a2
	a3=afac*a3
	a4=afac*a4

C  Normalize 2D spectrum by the average based on weakest 25% of time
C  slices, smoothed, and with local peaks removed.

	do i=1,NQ
	   do j=1,nz
	      s2(i,j)=s2(i,j)/max(psavg(i+5),0.01*smaxx)
	   enddo
	enddo

C  Find average of active spectral region, over the whole file.
	sum=0.
	do i=9,52
	   do j=1,nz
	      sum=sum+s2(i,j)
	   enddo
	enddo
	sigma=sum/(44*nz)

 999	return
	end