function out = fmst(in, p, angle)
%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.
%
%         22April2000 initial production version...K.Metzger

layers = p.layers;           % layers in transform
LL = 2^layers;               % sequence length plus 1
L = LL-1;                    % sequence length
LL2 = LL/2;
[rows, cols] = size(in);     % determine dimensions of input

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

% compute bias (aka pedestal) correction factor

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

out =zeros(rows,cols);       % create output using desired size
for c=1:cols                 % loop on the number of input columns
  
  % 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
  
  work(p.in,:) = [zeros(1,rows/L); reshape(in(:,c), rows/L,L).'];
  %tic; 
  
  % loop on the number of tensor products needed to form the
  %   Hadamard matrix using the 2x2 version
  
  for lay_ctr = 1:layers
    
    % the algorithm is layer invariant and does all
    %   sample per bit interleaves simultaneously
    %
    % work = [work(1:2:L,:)+work(2:2:LL,:); work(1:2:L,:)-work(2:2:LL,:)];
    % but runs faster if done in stages
    
    tempa = work(1:2:L,:);
    tempb = work(2:2:LL,:); 
    work(1:LL2,:)    = tempa+tempb;
    work(LL2+1:LL,:) = tempa-tempb;
    
    % processing can also be implemented as
    %
    % work([1:2:L 2:2:LL],:) = ...
    % [work(1:LL/2,:)+work(LL/2+1:LL,:); work(1:LL/2,:)-work(LL/2+1:LL,:)];
    % but runs faster if done in stages
    
    %tempa = work(1:LL2,:);
    %tempb = work(LL2+1:LL,:);
    %work(1:2:L,:)  = tempa+tempb;
    %work(2:2:LL,:) = tempa-tempb;

  end
  %toc
  
  % correct the bias (dc offset) for each column based on each
  %   column's sum unscrambling the bit order at the same time
  
  work = work(p.out,:)-kron(ones(LL,1), work(1,:))*factor;
  
  % place the result into the output array doing a reinterleaving
  %   of samples
  
  out(:,c) = reshape(work(2:LL,:).', rows,1);
  
end

return;