% randsig() generates a random real signal based on an input power spectrum.
%
% The power spectrum will be interpreted as having frequencies increasing 
% monotonically from 0 to the Nyquist frequency.
%
% If the power spectrum has N elements, then the signal will have N elements.
% Note that the program actually generates a 2(N-1)-element signal and returns
% only half. This suppresses periodicity in the result.
%
function signal = randsig(powerspec)

N = numel(powerspec)
thetas = 2*pi*rand(1, N);
phases = exp(i*thetas);
phases(1) = 1;
phases(N) = 1;
signalf(1:N) = sqrt(powerspec).*phases;
signalf(2*N-2:-1:N+1) = conj(signalf(2:N-1));
   % The redundant part of the transform of a real signal
signalf;
signal = real(ifft(signalf));
signal = signal(1:N-1);