WSJT-X/spec.f90

!---------------------------------------------------- End Module Audio1

!---------------------------------------------------- spec
subroutine spec(brightness,contrast,logmap,ngain,nspeed,a)

! Called by SpecJT in its TopLevel Python code.  
! Probably should use the "!f2py intent(...)" structure here.

! Input:
  integer brightness,contrast   !Display parameters
  integer ngain                 !Digital gain for input audio
  integer nspeed                !Scrolling speed index
! Output:
  integer*2 a(225000)           !Pixel values for 750 x 300 array

  real psa(750)                 !Grand average spectrum
  real ref(750)                 !Ref spect: smoothed ave of lower half
  real birdie(750)              !Spec (with birdies) for plot, in dB
  real variance(750)            !Variance in each spectral channel

  real a0(225000)               !Save the last 300 spectra
  integer*2 idat(11025)         !Sound data, read from file
  integer nstep(5)
  integer b0,c0
  real x(4096)                  !Data for FFT
  complex c(0:2048)             !Complex spectrum
  real ss(1024)                 !Bottom half of power spectrum
  logical first
  include 'gcom1.f90'
  include 'gcom2.f90'
  include 'gcom3.f90'
  include 'gcom4.f90'
  data jz/0/                    !Number of spectral lines available
  data nstep/15,10,5,2,1/       !Integration limits
  data first/.true./

  equivalence (x,c)
  save

  if(first) then
     call zero(ss,nq)
     istep=2205
     nfft=4096
     nq=nfft/4
     df=11025.0/nfft
     fac=2.0/10000.
     nsum=0
     iread=0
     first=.false.
     b0=-999
     c0=-999
     logmap0=-999
     nspeed0=-999
     nx=0
     ncall=0
     jza=0
     rms=0.
  endif

  nmode=1
  if(mode(1:4).eq.'JT65') nmode=2
  if(mode.eq.'Echo') nmode=3
  if(mode.eq.'JT6M') nmode=4

  nlines=0
  newdat=0
  npts=iwrite-iread
  if(ndiskdat.eq.1) then
     npts=jzc/2048
     npts=2048*npts
     kread=0
     if(nspeed.ge.6) then
        call hscroll(a,nx)
        nx=0
     endif
  endif
  if(npts.lt.0) npts=npts+nmax
  if(npts.lt.nfft) go to 900               !Not enough data available

10 continue
  if(ndiskdat.eq.1) then
! Data read from disk
     k=kread
     do i=1,nfft
        k=k+1
        x(i)=0.4*d2c(k)
     enddo
     kread=kread+istep                       !Update pointer
  else
! Real-time data
     dgain=2.0*10.0**(0.005*ngain)
     k=iread
     do i=1,nfft
        k=k+1
        if(k.gt.nmax) k=k-nmax
        x(i)=0.5*dgain*y1(k)
     enddo
     iread=iread+istep                       !Update pointer
     if(iread.gt.nmax) iread=iread-nmax
  endif

  sum=0.                                     !Get ave, rms of data
  do i=1,nfft
     sum=sum+x(i)
  enddo
  ave=sum/nfft
  sq=0.
  do i=1,nfft
     d=x(i)-ave
     sq=sq+d*d
     x(i)=fac*d
  enddo
  rms1=sqrt(sq/nfft)
  if(rms.eq.0) rms=rms1
  rms=0.25*rms1 + 0.75*rms
  
  if(ndiskdat.eq.0) then
     level=0                                    !Compute S-meter level
     if(rms.gt.0.0) then                        !Scale 0-100, steps = 0.4 dB
        dB=20.0*log10(rms/800.0)
        level=50 + 2.5*dB
        if(level.lt.0) level=0
        if(level.gt.100) level=100
     endif
  endif

  if(nspeed.ge.6) then
     call horizspec(x,brightness,contrast,a)
     ncall=Mod(ncall+1,5)
     if(ncall.eq.1 .or. nspeed.eq.7) newdat=1
     if(ndiskdat.eq.1) then
        npts=jzc-kread
     else
        npts=iwrite-iread
        if(npts.lt.0) npts=npts+nmax
     endif
     if(npts.ge.4096) go to 10
     go to 900
  endif

  call xfft(x,nfft)

  do i=1,nq                               !Accumulate power spectrum
     ss(i)=ss(i) + real(c(i))**2 + imag(c(i))**2
  enddo
  nsum=nsum+1

  if(nsum.ge.nstep(nspeed)) then      !Integrate for specified time
     nlines=nlines+1
     do i=225000,751,-1               !Move spectra up one row
        a0(i)=a0(i-750)               ! (will be "down" on display)
     enddo
     if(ndiskdat.eq.1 .and. nlines.eq.1) then
        do i=1,750
           a0(i)=255
        enddo
        do i=225000,751,-1
           a0(i)=a0(i-750)
        enddo
     endif

     if(nflat.gt.0) call flat2(ss,1024,nsum)

     do i=1,750                       !Insert new data in top row
        j=i+182                       ! ?? was 186 ??
        a0(i)=5*ss(j)/nsum
        xdb=-40.
        if(a0(i).gt.0.) xdb=10*log10(a0(i))
     enddo
     nsum=0
     newdat=1                          !Flag for new spectrum available
     call zero(ss,nq)                  !Zero the accumulating array
     if(jz.lt.300) jz=jz+1
  endif

  if(ndiskdat.eq.1) then
     npts=jzc-kread
  else
     npts=iwrite-iread
     if(npts.lt.0) npts=npts+nmax
  endif

  if(npts.ge.4096) go to 10

!  Compute pixel values 
  iz=750
  logmap=0
  if(brightness.ne.b0 .or. contrast.ne.c0 .or. logmap.ne.logmap0 .or.    &
          nspeed.ne.nspeed0 .or. nlines.gt.1) then
     iz=225000
     gain=40*sqrt(nstep(nspeed)/5.0) * 5.0**(0.01*contrast)
     gamma=1.3 + 0.01*contrast
     offset=(brightness+64.0)/2
     b0=brightness
     c0=contrast
     logmap0=logmap
     nspeed0=nspeed
  endif

!  print*,brightness,contrast,logmap,gain,gamma,offset
  do i=1,iz
     n=0
     if(a0(i).gt.0.0 .and. logmap.eq.1) n=gain*log10(0.001*a0(i)) + offset + 20
     if(a0(i).gt.0.0 .and. logmap.eq.0) n=(0.01*a0(i))**gamma + offset
     n=min(252,max(0,n))
     a(i)=n
  enddo

900 continue
  return
end subroutine spec