function [y, h] = randc(nSeq, Ryy)
%
% function [y h] = randc(nSeq, Ryy)
%
%   nSeq is the length of the sequence.
%   Ryy is the desired autocorrelation of the sequence (two-sided)
%   y is the random sequence
%   h is the FIR fitler (length(h) = length(Ryy))
%
% Produces a zero-mean, unit-variance correlated random sequence 
% with the given autocorrelation function Ryy.  
%
% Notes:
%   1. The length of Ryy must be odd.
%
% The sequence is generated by constructing
% a constant-phase FIR filter: H(f) = sqrt(Syy/Sxx) where
% Syy = fft(Ryy) and Sxx = fft(Rxx).  Rxx is the autocorrelation
% of white noise, ie all zeros except for 1 in the middle.
% H(f) is then converted to the time domain, h(t) = ifft(H(f))
% to yeild the coefficients of the FIR filter.  A zero-mean,
% unity-variance sequence x is generated and filtered by h
% to generate y.
%
% While magnitude of H(f) is determined, the phase is arbitrary.
% This makes the computation of H(f) ambiguous, and 
% many filters will yield the same Ryy.  Here, the real and imaginary 
% parts of H(f) are chosen to be equal. This results in a constant-
% phase filter.
%
% Author: Douglas N. Greve
% Date:   Aug 17, 1998
%
%
%


%
% fmri_randc.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:06 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

% Testing and usage: the code below shows how to use and test this
% function.  It generates a desired autocorrelation Ryy then uses
% this function to generate a random sequence.  The matlab xcorr
% function is used to generate the emperical autocorrelation function
% which is then plotted against the desired.  The results show
% that the difference between the desired and emperical goes to
% zero as the length of the sequence increases.
%   nSeq = 1024   
%   nEst = 8;     
%   Ryy = fir1(2*nEst,.25)';
%   Ryy = Ryy/max(Ryy);
%   y = randc(nSeq, Ryy);
%   nCoeff = length(Ryy);
%   Ryy2 = xcorr(y,(nCoeff-1)/2, 'biased');
%   Ryy2 = Ryy2/max(Ryy2);
%   nx = [1:nCoeff];
%   plot(nx,Ryy,'g',  nx,Ryy2,'r');

% Check that the length of Ryy is odd %
nCoeff = length(Ryy);
if( mod(nCoeff,2) == 0)
  error('Ryy must have an odd number of elements');
  return;
end

% Check that the sequence length is > 0 %
if( nSeq < 1)
  error('Sequence Length must be > 0');1
  return;
end

% normalilze and make sure dim is correct %%
if(size(Ryy,1) > size(Ryy,2))
  Ryy = Ryy/max(Ryy);
else 
  Ryy = Ryy'/max(Ryy);
end

% autocorrelation of uncorrelated noise %
Rxx = zeros(nCoeff,1);
Rxx((nCoeff-1)/2) = 1;

% power specta, use 2*nCoeff to avoid wrap-around
Sxx = fft(Rxx,2*nCoeff);
Syy = fft(Ryy,2*nCoeff);

% Compute the magnitude of H(f) %
mH = sqrt(abs(Syy ./ Sxx));
clear Sxx Syy Rxx;

% Compute Complex H(f) with equal Real and Imaginary Parts,
% Note:  |H|^2 = mH^2
H = mH - mH *1i; % constant phase = pi/4

% Convert filter to the time-domain 
htmp = fftshift(real(ifft(H)));
[hmax m] = max(htmp);
n = (nCoeff-1)/2;

% keep only nCoeff of the filter coefficients 
h  = zeros(nCoeff,1);
h = htmp([m-n:m+n]);
clear H htmp;

% generate a white noise sequence, use nSeq+2*nCoeff to avoid edge effects
x    = randn(nSeq+2*nCoeff,1);

% filter using the FIR
ytmp = filter(h,1,x);

% extract the nSeq samples from the middle
y = ytmp([nCoeff+1:nSeq+nCoeff]);
clear ytmp;

% Remove the mean
y = y - mean(y);

return;	

