program spectrum

      program spectrum

      implicit integer(a-z)
      parameter (maxlen=512)
      integer over

      integer*2 status

      byte ss(maxlen)
      character*(maxlen) cdata

      logical average

      external get_wave

      equivalence (ss,cdata)
c
           
      over = 10
      average = .true.
    1 continue
      ldata = get_wave(cdata,samplesPsec)
c      ldata = get_sine(ss,maxlen,samplesPsec)
c      ldata = get_square(ss,maxlen,samplesPsec)
c      ldata = get_rndm(ss,maxlen,samplesPsec)
      if(ldata.le.0) goto 99
      call do_fft(ss,ldata,samplesPsec,average,over)
      goto 1
    2 write(6,*) ' Finished'
      stop 0
   99 write(6,*) ' Error ',ldata,' reading data'
      end


*=*=*=*= get_rndm.html =*=*=*=*

integer function get_rndm

      integer function get_rndm(ss,maxlen,samplesPsec)
      integer samplesPsec

      byte ss(maxlen)
      real val

      get_rndm = maxlen

      samplesPsec = maxlen

      do i=1,maxlen

         val = 127.*random(i)
         ss(i) = val

      end do
      end   

*=*=*=*= get_sine.html =*=*=*=*

integer function get_sine

      integer function get_sine(ss,maxlen,samplesPsec)
      integer samplesPsec

      byte ss(maxlen)
      real val,v1,v2,pi2

      save icall,pi2

      data icall /0/
c
      if(icall.eq.0) then
         icall = 1
         write(6,*) 'Enter frequency 1 in Hertz'
         read(5,*) hertz1

         write(6,*) 'Enter frequency 2 in Hertz'
         read(5,*) hertz2

         pi2 = 6.0*asin(sqrt(3.)*0.5)
      endif
      get_sine = maxlen

      samplesPsec = 7.1*max(hertz1,hertz2)
      v1 = pi2*real(Hertz1)/real(samplesPsec)
      v2 = pi2*real(hertz2)/real(samplesPsec)
      do i=1,maxlen

         val = 64.*(2.+ cos(v1*real(i-1)) +
     &                   cos(v2*real(i-1)))
         ss(i) = val

      end do
      end   
      

*=*=*=*= get_square.html =*=*=*=*

integer function get_square

      integer function get_square(ss,maxlen,samplesPsec)
      integer samplesPsec

      byte ss(maxlen)
      real val,v1,v2,pi2

      save icall,pi2

      data icall /0/
c
      if(icall.eq.0) then
         icall = 1
         write(6,*) 'Enter frequency 1 in Hertz'
         read(5,*) hertz1

         write(6,*) 'Enter frequency 2 in Hertz'
         read(5,*) hertz2

         pi2 = 6.0*asin(sqrt(3.)*0.5)
      endif
      get_square = maxlen

      samplesPsec = 7.1*max(hertz1,hertz2)
      v1 = pi2*real(Hertz1)/real(samplesPsec)
      v2 = pi2*real(hertz2)/real(samplesPsec)
      do i=1,maxlen

         a = -1.
         if(cos(v1*real(i-1)).gt.0.) a = 1.
         b = -1.
         if(cos(v2*real(i-1)).gt.0.) b = 1.
         if(hertz2.eq.0) b = 0.
         val = 64.*(2.+a+b)
         ss(i) = val

      end do
      end   
      
      
      

*=*=*=*= get_wave.html =*=*=*=*

integer function get_wave

      integer function get_wave(ss,samplesPsec)
      implicit integer(a-z)
      parameter (recl=512)
      integer*2 status

      character*127 cin

      character*(*) ss

      logical info

c
      structure /format_chunk/
       union
        map
         character*4 fmt

         integer*4   Length

         integer*2   FormatTag

         integer*2   Channels

         integer*4   SamplesSec

         integer*4   BytesSec

         integer*2   BlockSize

         integer*2   BitsSample

        end map
        map
         character*24 read

        end map
       end union
      end structure
      record /format_chunk/ fmt
c
      structure /header_chunk/
       union
        map
         character*4 riff

         integer*2   data1

         integer*2   data2

         character*4 filetype

        end map
        map
         character*12 read

        end map
       end union
      end structure
      record /header_chunk/ header
c
      structure /data_chunk/
       union
        map
         character*4  data

         integer*2    data1

         integer*2    data2

        end map
        map
         character*8  read

        end map
       end union
      end structure
      record /data_chunk/ datac
c
      save icall,info

c
      data icall /0/
c
      if(icall.eq.0) then
         icall = 1
         info = .false.
         call getarg(2,cin,status)
         if(status.gt.0.and.cin(:1).eq.'I') info = .true.
         call getarg(1,cin,status)
         if(status.le.0) 
     &     stop 'Usage : SPECTRUM <.wav file name> '
         if(cin(1:1).eq.'?') 
     &     stop 'Usage : SPECTRUM <.wav file name> '
c
         open(1,file=cin(:status),err=99,status='old',
     &        iostat=iostat,form='binary',recl=recl)
c
         lread = len(header.read)+len(fmt.read)+len(datac.read)
         len_data = recl-lread

         read(1) header.read,fmt.read,datac.read,ss(1:len_data)
         if(info) then
            write(6,*) 'RIFF      ',header.riff
            write(6,*) 'data1     ',header.data1
            write(6,*) 'data2     ',header.data2
            write(6,*) 'filetype  ',header.filetype
            write(6,*) 'fmt       ',fmt.fmt
            write(6,*) 'fmt length',fmt.length
            write(6,*) 'fmt tag   ',fmt.formattag
            write(6,*) 'Channels  ',fmt.channels
            write(6,*) 'Samples   ',fmt.samplessec
            write(6,*) 'Bytes     ',fmt.bytessec
            write(6,*) 'Blocksize ',fmt.blocksize
            write(6,*) 'Bits/smple',fmt.bitssample
            write(6,*) 'Data      ',datac.data
            write(6,*) 'Data1     ',datac.data1
         endif
         if(fmt.formattag.ne.1.or.
     &      header.filetype.ne.'WAVE') then
            stop 'Not a .wav file'
         endif
         get_wave = len_data

         samplesPsec = fmt.samplessec
         return
      endif
      iostat = 5000
      read(1,end=99,err=99) ss

      get_wave = recl

      return
   99 get_wave = -iostat

      end 



c
*=*=*=*= do_fft.html =*=*=*=*

subroutine do_fft

      subroutine do_fft(data,n,rate,average,count_over)
      implicit integer(a-z)
      parameter (segment=256,length=segment/2)
      complex in(segment)
      real d,sum(length),out(length),rdata(segment)
      byte data(n)
      integer rate,count,count_over

      logical average

      save count,sum

      data count /0/
   
      if(count.eq.0) then
         do i=1,length

            sum(i) = 0.0
         end do
      endif
c
      ipos = 0
    1 if(ipos+segment.gt.n) goto 2
      do i=1,segment

         d = real(iand(data(ipos+i),#FF))
         rdata(i) = d

         in(i) = cmplx( d , 0.0 )
      end do
      call fft(in,segment,-1)
c
c average
c
      count = count+1
      if(average) then
         do i=1,length

            sum(i) = sum(i)+abs(in(i))
            out(i) = sum(i)/real(count)
         end do
      else
         do i=1,length

            out(i) = abs(in(i))
         end do
         count_over = 0
      endif
c
      call display(rdata,out,segment,rate,count_over,count)
      ipos = ipos + segment

      if(average.and.count.ge.count_over) count = 0
      goto 1
    2 continue
c
      end

*=*=*=*= DISPLAY.html =*=*=*=*

SUBROUTINE DISPLAY

      SUBROUTINE DISPLAY(in,out,segment,rate,count_over,inum)
      include 'fgraph.fd'

c
      common /plot/ hz_low,hz_high,db_low,db_high,vlow,vhigh

      real hz_low,hz_high,db_low,db_high,vlow,vhigh

      integer rate,segment,length,count_over,inum

      integer*2 halfy,fullx,fully,xborder,yborder

      integer*2 status,handle1

      real in(segment),out(segment)
      character*5 ctext

      byte diagmask(8)
c
      record /videoconfig/ myscreen
      record /wxycoord/     xy
c
      save icall,halfy,fullx,fully,xborder,yborder

      data icall /0/
      data diagmask /8*#FF/
c
      length = segment/2
      if(icall.eq.0) then
         icall = 1
         hz_low = 0.0
         hz_high = real(rate)*0.5
         db_low = 0.0
         db_high = 10.0
         vlow = 0.
         vhigh = 255.
         status = setvideomode($MAXRESMODE)
         if(status.eq.0) stop 'Unable to start Graphics'
         handle1 = wggetactiveqq()
         status = wgcloseqq(handle1)
         handle1 = wgopenqq('Frequency Response'//char(0))
         status = wgsetactiveqq(handle1)
         status = setvideomode($VRES16COLOR)
c         status = clickqq(QWIN$TILE)
         call getvideoconfig(myscreen) 
         halfy = myscreen.numypixels/2
         fullx = myscreen.numxpixels-1
         fully = myscreen.numypixels-1
         xborder = fullx/5
         yborder = fully/8
         status = registerfonts("dummy")
         status = setfont('h8w8b')
         if(status.lt.0) then
c            write(6,*) ' Setting substitute font'
            status = setfont('f')
         endif
c
c Set axes and tick marks
c
         call setviewport(xborder,fully-yborder,fullx,
     &                    fully)
         status = setwindow(.true.,0.,0.,real(segment),1.)
         status = setcolor(6)
         call moveto_w(0.,1.,xy)
         status = lineto_w(real(segment),1.)
         ic = 0 
         do i=1,segment,10
            ic = ic+1
            call moveto_w(real(i),1.,xy)
            status = lineto_w(real(i),0.95)
            if(mod(ic,2).eq.1) then
               status=lineto_w(real(i),0.9) 
               write(ctext,'(f5.3)') real(i)/real(rate)
               call verttext(ctext)
            endif
         end do
         status = setcolor(8)
         call moveto_w(0.,0.10,xy)
         call outgtext(' Elapsed seconds from sample start')

         call setviewport(xborder,halfy-yborder+1,fullx,halfy)
         status = setwindow(.true.,hz_low,0.,hz_high,1.)
         status = setcolor(6)
         call moveto_w(hz_low,1.,xy)
         status = lineto_w(hz_high,1.)
         mhzh = nint(hz_high/1000.)*1000
         mhzl = ifix(hz_low/100.)*100
         ic = 0
         do i=mhzl,mhzh,100
            ic = ic + 1
            d = real(i)
            call moveto_w(d,1.,xy)
            status = lineto_w(d,0.95)
            if(mod(ic,2).eq.1) then
               status = lineto_w(d,0.9)
               write(ctext,'(i5)') i

               call verttext(ctext)
            end if
         end do
         call moveto_w(hz_low,0.10,xy)
         status = setcolor(8)         
         call outgtext(' Hertz')
         call setfillmask(diagmask)
      endif
c
c called every time for vertical axes and info
c
      call setviewport(0,0,xborder-1,fully)
      status = setwindow(.true.,0.,0.,1.,1.)
      status = setcolor(0)
      status = rectangle_w($GFILLINTERIOR,0.,0.,1.,1.)      
      status = setcolor(10)
      write(ctext,'(i5)') segment

      call moveto_w(0.,0.95,xy)
      call outgtext(ctext//' channels')
      call moveto_w(0.,0.90,xy)
      write(ctext,'(f5.1)') db_high

      call outgtext(' dB max '//ctext)
      call moveto_w(0.,0.85,xy)
      write(ctext,'(i5)') inum

      call outgtext('Sample '//ctext)
      call moveto_w(0.,0.80,xy)
      if(count_over.ne.0) then
         write(ctext,'(i5)') count_over

         call outgtext('out of '//ctext)
      endif



      call setviewport(xborder,halfy+yborder,fullx,
     &                 fully-yborder)
      status = setwindow(.true.,0.,vlow,real(segment),vhigh)
      status = setcolor(0)
      status = rectangle_w($GFILLINTERIOR,0.,vlow,real(segment),vhigh)
      call moveto_w(0.,vlow,xy)
      status = setcolor(3)
      do i=1,segment

         d = real(i)
         status = lineto_w(d,in(i))
      end do

      
      call setviewport(xborder,0,fullx,halfy-yborder)
      status = setwindow(.true.,hz_low,db_low,hz_high,db_high)
      channel_to_hz = real(rate)/real(segment)
      barw2 = channel_to_hz*0.5
      barhm = db_high

      call setfillmask(diagmask)
      amax = -9999.
      offset = 10000.
      if(out(1).gt.0.) offset = 10000./out(1)
      do i=2,length

         dB = out(i)
         if(dB.gt.amax) amax = dB

         d = real(i)*channel_to_hz-barw2*0.5
         status = setcolor(0)
         status = 
     &   rectangle_w($GFILLINTERIOR,d,db_low,d+barw2,db_high)
         status = setcolor(2)
         status = 
     &   rectangle_w($GFILLINTERIOR,d,db_low,d+barw2,dB)
      end do      
c
c Set scale to 1.5 times max amplitude for next sample
c
      db_high = 1.5*amax

      end


c
c  The following FFT code is from:
c
c  Arthur Wouk (wouk@brl-vgr)
c
c  *******************
c  Fast Fourier Transform
c  *******************
c
*=*=*=*= FFT.html =*=*=*=*

SUBROUTINE FFT

        SUBROUTINE FFT (X, N, K)
C       FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X
C       (A COMPLEX ARRAY OF DIMENSION N). SOURCE: Ferziger; Numerical
C       methods for engineering applications.
C
C       X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM.
C       N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<32769). C K="1" FOR FORWARD TRANSFORM. C K="-1" FOR INVERSE TRANSFORM. C IMPLICIT INTEGER (A-Z) INTEGER SBY2,S

        REAL GAIN, PI2, ANG, RE, IM

        COMPLEX X(N), XTEMP, T, U(16), V, W

        LOGICAL NEW

        DATA    PI2,GAIN,NO,KO /6.283185307, 1., 0, 0/
C
C       TEST FIRST CALL?
C
        NEW = ( NO .NE. N)
        IF ( .NOT. NEW) GO TO 2
C
C       IF FIRST CALL COMPUTE LOG2 (N).
C
        L2N = 0
        NO = 1
    1   L2N = L2N + 1
        NO = NO + NO

        IF (NO .LT. N) GO TO 1
        GAIN = 1./N

        ANG = PI2*GAIN

        RE = COS (ANG)
        IM = SIN (ANG)
C
C       COMPUTE COMPLEX EXPONENTIALS IF NOT FIRST CALL
C
    2   IF (.NOT. NEW .AND. K*KO .GE. 1) GO TO 4
        U(1) = CMPLX (RE, -SIGN(IM, FLOAT(K)))
        DO 3 I = 2,L2N

           U(I) = U(I-1)*U(I-1)
    3   CONTINUE
        K0 = K

C
C       MAIN LOOP
C
    4   SBY2 = N

        DO 7 STAGE = 1,L2N

           V = U(STAGE)
           W = (1., 0.)
           S = SBY2

           SBY2 = S/2
           DO 6 L = 1,SBY2

                DO 5 I = 1,N,S

                   P = I + L- 1
                   Q = P + SBY2

                   T = X(P) + X(Q)
                   X(Q) = ( X(P) - X(Q))*W

                   X(P) =T

    5           CONTINUE
                W = W*V

    6      CONTINUE
    7   CONTINUE
C
C       REORDER THE ELEMENTS BY BIT REVERSAL
C
        DO 9 I = 1,N

           INDEX = I-1
           JNDEX = 0
           DO 8 J = 1,L2N

                JNDEX = JNDEX+JNDEX

                ITEMP = INDEX/2
                IF (ITEMP+ITEMP .NE. INDEX) JNDEX = JNDEX + 1
                INDEX = ITEMP

    8      CONTINUE
           J = JNDEX + 1
           IF (J .LT. I) GO TO 9
           XTEMP = X(J)
           X(J) = X(I)
           X(I) = XTEMP

    9   CONTINUE
C
C       FORWARD TRANSFORM DONE
C
        IF (K .GT. 0) RETURN
C
C       INVERSE TRANSFORM
C
        DO 10 I = 1,N

           X(I) = X(I)*GAIN

   10   CONTINUE
        RETURN
        END
*=*=*=*= verttext.html =*=*=*=*

SUBROUTINE verttext

      SUBROUTINE verttext (text)
      
      INCLUDE 'FGRAPH.FD'


      CHARACTER*(*)  text

      CHARACTER*1    ch

      INTEGER*2      status, loop, ipos, jpos, indent

      RECORD         / xycoord / xy, xyorig
      RECORD         / fontinfo / fi

C     Get the beginning location for output
      CALL getcurrentposition(xyorig)
      ipos = xyorig.xcoord
      jpos = xyorig.ycoord
      
C     The font info is used to center the letters, and determine how
C       far each letter should be below the previous one
      status = getfontinfo(fi)
      
C     Output one char of the string at a time, in a descending column
      do loop = 1, len(text)
        ch = text(loop:loop)
        indent = ipos + (fi.pixwidth-getgtextextent(ch))/2
        CALL moveto(indent, jpos, xy)
        CALL outgtext(ch)
        jpos = jpos + fi.pixheight
      end do
      
C     Set new current position to the end of the vertical text column
      CALL moveto(ipos, jpos, xy)
      END