function [x,t]=ransiggen(ns,nt,T,frsd);



% Generates ns random signals with given power spectral density
% nt is the number of points (should be a power of 2)
% T is the total time for the signal
% frsd is a nb-by-2 matrix where nb is the number of points 
% defining the power spectral density.  The rows of frsd are
% the frequency and amplitude of each point defining the PSD
%


dt=T/nt;

t=linspace(0,T,nt);
df=1/(nt*dt);
fmax=1/(2*dt);

fr=frsd(:,1)';
sd=frsd(:,2)';
nb=length(fr);


if fr(nb) ~= fmax
  fr(nb+1) = fmax;
  sd(nb+1) = 0;
end

if fr(1) ~= 0
  fr=[0 fr];
  sd=[0 sd];
end

nt2=nt/2;
i=1:nt2;
mag=interp1(fr,sd,i*df);


%
%Interpolate (linear) the intermediate target psd values and convert to Fourier
%spectrum amplitudes
%


c=sqrt(df*mag);
for i=1:1:ns
  
  %Form Random phase angles
   ph=rand(1,nt2);
   ph=2*pi*ph-pi;
	
  %Form Fourier Transform of random signal
   ii=sqrt(-1);
   xft=cos(ph).*c+(sin(ph).*c)*ii;
   xft=[xft,0,conj(xft(nt2:-1:2))];
	
  %Inverse Fourier Transform to obtain time-domain signal
   xft=ifft(xft);

  %Calculate Error
   er=sum(imag(xft))/nt;
   %disp(['error in imaginary part = ',num2str(er)]);
   x(i,:)=real(xft);
   asd(i,:)=(fft(xft).*conj(fft(xft)))/df;
end
asd=asd(:,1:nt2);
freq=[0:df:((nt2-1)*df)];

