fmst.py

from numpy import *
def fmst(thein=None, p=None, angle=None):
#Function to implement fast m-sequence transform.
#
# out = fmst(in, p, angle)
#
# in      a column vector or a matrix whose columns consist of
#         samples of one period of a m-sequence based baseband
#         waveform.  The number of samples corresponding to
#         a sequence bit must be an integer value.
#
# p       a structure generated by the fmtsetup function.  This
#         is sequence law dependent and provides the fmt 
#         function information needed specific to the sequence.
#
# angle   modulation angle used in generating the input data
#         set.  This is in degrees.  Values greater than 180
#         result in the use of the "matched" angle, 
#         angle=atan(sqrt(L)) where L is the number of bits
#         in the sequence period.
#
#         The basic m-sequence transform technique is described
#         in the paper "On Fast M-Sequence Transforms" by
#         M.Cohn and A.Lempel, IEEE Transactions on Information
#         Theory, January 1977.
#
#         Hadamard processing based on T.G.Birdsall's pass invariant
#         fast Hardamard algorithm.  See also "Algorithms for
#         Discrete Fourier Transform and Convolution, 2nd Ed." 
#         by R.Tolimieri, M.An, and C.Lu, Springer-Verlag, 1997.
#         Chapter 2 discusses tensor product relationships.
#
#         Based on fmst.m by K.Metzger

# thein is supposed to be ALWAYS a vector... 
    out = []
    layers = p["layers"]     # layers in transform
    p_in   = p["in" ].astype(int)
    p_out  = p["out"].astype(int)
    LL = 2**layers           # sequence length plus 1
    L = LL - 1               # sequence length
    LL2 = int( LL/2 )

    rows = thein.size
    cols = 1

    if ( angle > 180 ): # if > 180 degrees use matched angle
       ta = sqrt(L)
    else:
       ta = tan(pi*angle/180) # otherwise use specified angle

# compute bias (aka pedestal) correction factor

    factor = ( ta*ta - L + 1j*ta*(L+1) )/( L*L + ta*ta )

#    out = zeros((rows,cols))  # create output using desired size
    rL = int( rows/L )
    out = zeros(rows)
    inresh = thein.reshape(L,rL)
  # extract a column from the input
  # place a zero at the column top
  # reshape to obtain one sample per bit and a number of
  # columns equal to the number of samples per bit
  # scramble values to allow use of Hadamard processing
    w0   = zeros((LL,rL)) + 1j*zeros((LL,rL))
    work = zeros((LL,rL)) + 1j*zeros((LL,rL))
    w0[1:,:] = inresh
    work[p_in-1,:] = w0 # O.K.
  # save('thework',work)
  
  # loop on the number of tensor products needed to form the
  # Hadamard matrix using the 2x2 version

    for lay_ctr in range(layers):
        tempa = copy( work[0:L :2,:] )
        tempb = copy( work[1:LL:2,:] ) #; print( tempb[:,0] )
        work[0:LL2 ,:] = tempa + tempb
        work[LL2:LL,:] = tempa - tempb
  # correct the bias (dc offset) for each column based on each
  # column's sum unscrambling the bit order at the same time
#    save('thework',work) # O.K.
#  work = work(p.out,:)-kron(ones(LL,1), work(1,:))*factor;
    work = work[p_out-1,:] - kron( ones((LL,1)), work[0,:] )*factor
  # save('thework',work) # O.K.
  # place the result into the output array doing a reinterleaving of samples
    out = ( work[1:LL+1,:] ).reshape((rows,1))

    return out