Grid Enabled Analysis The GIOD Project General GIOD Partners Caltech's Centre for Advanced Computing Research Caltech's HEP department
|
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