Subroutine F_PREPARE_MLS(mls,ierr)
!
! MLS sequence generator for up to 4 taps in a 24 bit register
!
! (c) Julian J. Bunn, 1999, julianb@altavista.net
      parameter (maxelements=24)
      integer bit(maxelements)
      integer tap(4,maxelements)
      integer ntap(maxelements)
      data NTAP / 0, 2, 2, 2, 2, 2, 2, 4, 2, 2, 2, 4, &
     &            4, 4, 2, 4, 2, 2, 4, 2, 2, 2, 2, 4/
      data TAP  / 0, 0, 0, 0, &
     &            1, 2, 0, 0, &
     &            1, 3, 0, 0, & !!! This used to be 2,3 !
     &            3, 4, 0, 0, &
     &            3, 5, 0, 0, &
     &            5, 6, 0, 0, &
     &            6, 7, 0, 0, &
     &            2, 3, 5, 8, &
     &            5, 9, 0, 0, &
     &            7,10, 0, 0, &
     &            9,11, 0, 0, &
     &            6, 8,11,12, &
     &            9,10,12,13, &
     &            4, 8,13,14, &
     &           14,15, 0, 0, &
     &            4,13,15,16, &
     &           14,17, 0, 0, &
     &           11,18, 0, 0, &
     &           14,17,18,19, &
     &           17,20, 0, 0, &
     &           19,21, 0, 0, &
     &           21,22, 0, 0, &
     &           18,23, 0, 0, &
     &           17,22,23,24/
      integer*1 mls(*)
! First calculate the sequence: Initialise the "nelement" bit
! shift register with all 1s, and then generate each bit of
! the sequence
      do I=1,NELEMENTS
         BIT(I) = 1
      end do
      do I=1,LSEQUENCE
         SEQ(I) = -1
         if(BIT(NELEMENTS).eq.0) SEQ(I) = 1
         IBIT = 0
         do J=1,NTAP(NELEMENTS)
            IBIT = IBIT + BIT( TAP(J,NELEMENTS) )
         end do
         do J=NELEMENTS,2,-1
            BIT(J) = BIT(J-1)
         end do
         BIT(1) = MOD(IBIT,2)
      end do
! Initialise the Playback buffer with the MLS samples
      do I=1,LPLAYBUFFER        !  2*(lsequence+1)
         INDEX = MOD(I,LSEQUENCE+1)
         MLS(I) = DACMAX
         if(SEQ(INDEX).eq.-1) MLS(I) = DACMIN
      end do
! Construct the Noise matrix: each row is the MLS right shifted one place
! Also contruct the P1tag for each column in the matrix: the tag for a
! column is simply calculated by treating the first nelement numbers as
! the bits set in the tag
      do ICOL=1,LSEQUENCE
         P1TAG(ICOL) = 0
         P2TAG(ICOL) = 0
      end do
      do IROW=1,NELEMENTS
         do ICOL=1,LSEQUENCE
            INP = ICOL - IROW + 1
            if(INP.le.0) INP = LSEQUENCE + INP
            if(SEQ(INP).eq.-1) then
               P1TAG(ICOL) = IBSET(P1TAG(ICOL),IROW-1)
            endif
         end do
      end do
! With the Tag, we can determine the inverse of the Permutation matrix P1.
! inv(P1) is defined as being the matrix which permutes the columns of
! the Noise matrix so that they are in ascending order of Tag value
! Calculate the P2Tag ... this is computed from the rows of the Noise
! matrix after permutation with P1
      do K=0,NELEMENTS-1
         INDEX = 2**K
         do I=1,LSEQUENCE
            ICOL = P1TAG(I)
            if(ICOL.eq.INDEX) then
               do IROW=1,LSEQUENCE
                  ! for this row (irow) and column (i) find if Seq bit is set
                  INP = I - IROW + 1
                  if(INP.le.0) INP = LSEQUENCE + INP
                  if(SEQ(INP).eq.-1) then
                       P2TAG(IROW) = IBSET(P2TAG(IROW),K)
                  endif
               end do
               goto 100
            endif
         end do
  100    continue
      end do
! The next section of code is not used: it calculates the Sylvester-type
! Hadamard matrix of order (lsequence+1)
! The base Hadamard matrix (of order 2) from which the total matrix will be evol
!      Hadamard(1,1) = 1
!      Hadamard(1,2) = 1
!      Hadamard(2,1) = 1
!      Hadamard(2,2) = -1
!      do i=2,nelements
!         last = 2**(i-1)
!         do irow=1,last
!            do icol=1,last
!               Hadamard(last+icol,irow) = Hadamard(icol,irow)
!               Hadamard(icol,last+irow) = Hadamard(icol,irow)
!               Hadamard(last+icol,last+irow) = -Hadamard(icol,irow)
!            end do
!         end do
!      end do
!      write(6,*) ' Hadamard matrix:'
!      do i=1,lsequence+1
!         write(6,'(10(i3,1x))') (Hadamard(i,j),j=1,lsequence+1)
!      end do
End
Subroutine FHT(ss)
! Takes the measured response (ss) to the MLS, manipulates it, and
! then performs the FHT to retrieve the impulse response of the system
! (c) Julian J. Bunn, 1999, julianb@altavista.net
      integer*2 ss(*)
! Take the sequence, and permute it with P1
      do I=1,LSEQUENCE
         SEQ(1+P1TAG(I)) = SS(I)
      end do
      SEQ(1) = 0
! Perform the Hadamard Transform
      IHALF = (LSEQUENCE+1)/2
      do ISHIFT=1,NELEMENTS
         ! calculate the current column values by computing the butterfly
         do IROW=1,LSEQUENCE+1
            if(IROW.le.IHALF) then
               I2 = 2*IROW
               ROW(IROW) = SEQ(I2-1) + SEQ(I2)
            else
               I2 = 2*(IROW-IHALF)
               ROW(IROW) = SEQ(I2-1) - SEQ(I2)
            endif
         end do
         do IROW=1,LSEQUENCE+1
            SEQ(IROW) = ROW(IROW)
         end do
      end do
! Remove the first element, of the permuted, transformed vector, and permute it
      do IROW=1,LSEQUENCE
         SEQ(IROW) = SEQ(IROW+1)
      end do
        OSEQ = 1./REAL(LSEQUENCE)
      do ICOL=1,LSEQUENCE
         HEARD(ICOL) = NINT( OSEQ * REAL(SEQ(P2TAG(ICOL))) )
      end do
      do I=LSEQUENCE+1,LRECORDBUFFER
         HEARD(I) = 0
      end do
End