sim_genSignal(t,sigInfo,numComponents)

function sig = sim_genSignal(t,sigInfo,numComponents)
%% Simple signal generator

%% Modified by John Harwell to allow each signal to be of a different type
%% Feb. 25, 2001 Consolidation

%% Copyright (C) by Charles. H. Anderson (All Rights Reserved)
%% Dept. Anatomy and Neurobiology
%% Washington Univ. School of Medicine
%% St. Louis, MO
%% [e-mail]



% Modified on 10/29/01 by John Harwell
% 
%   Removed debugging print statements.
%
% Modified on 10/31/01 by John Harwell
%   Renamed from "new_genSignal" to "sim_genSignal".

Nt = length(t);
T0 = t(1);
T = t(end);
dt = t(2)-t(1);
lor = Nt * dt;   % length of run
if(strcmp(sigInfo(1).SigType(1:4),'file'))
    if(exist([sigInfo(1).Filename,'.mat']))
         load(sigInfo(1).Filename);
         sig = interp1(ExtTime, ExtSignal',t, 'linear');
         sig = sig';
     else
         errstr = sprintf('Cannot find signal data file %s\n',sigInfo.Filename);
         error(errstr);
     end
else
    %%%%% Get all the components 
    for j=1:numComponents
		switch lower(sigInfo(j).SigType(1:4))
		case 'cons'
            sig(j,:) = sigInfo(j).Amplitude * ones(1,Nt);
        case 'nois'
            rms = 0.5;
			N = 2*floor(lor/(2*dt))+1; % An odd number
			domega = 2*pi/lor;
			omega = domega*[0:(N-1)/2];
            parms = [sigInfo(j).Frequency  2*pi*sigInfo(j).FrequencyHigh  -1]; %% Change the seed for identical reruns
			Amps = BLimitWhiteNoise(omega,parms);
			Amps = [Amps,fliplr(conj(Amps(2:end)))];
			S = real(ifft(Amps));
			rmsS = sqrt(sum(S.^2)/N);
			sig(j,:) =  sigInfo(j).Amplitude*S*(rms/rmsS);
	
		case 'ramp'
            ton  = sigInfo(j).TimeOn;  %%% * (T-T0);
            toff = sigInfo(j).TimeOff; %%% * (T-T0);
            iOn  = ceil(ton/dt)+2;
            iOff = ceil(toff/dt)+1;
            num = iOff - iOn + 1;
            sig(j,:) = zeros(1,Nt);
            sig(j,:) = ones(1,Nt) * sigInfo(j).dcOffset;
            amp = sigInfo(j).dcOffset;
            damp = sigInfo(j).Amplitude / num;
            for k=iOn:iOff
                sig(j,k) = amp;
                amp = amp + damp;
            end
            if (iOff < Nt)
               for k=iOff+1:Nt
                  sig(j,k) = amp;
               end
            end
		case 'step'
            sig(j,:) = zeros(1,Nt);
            sig(j,:) = ones(1,Nt) * sigInfo(j).dcOffset;
            ton  = sigInfo(j).TimeOn ;% * (T-T0); 
            toff = sigInfo(j).TimeOff;% * (T-T0);
            iOn  = ceil(ton/dt)+2;
            iOff = ceil(toff/dt)+1;
            num  = iOff-iOn+1;
            sig(j,iOn:iOff) = (sigInfo(j).Amplitude + sigInfo(j).dcOffset) * ones(1,num);
		%%%% Modulated sinewave %%%%%%%%
		case 'sine'
            sig(j,:) = (sigInfo(j).Amplitude) ...
                     * sin(2*pi*sigInfo(j).Frequency*t  + (pi/180)*sigInfo(j).Phase);
            sig(j,:) = sig(j,:) + (ones(1,length(sig(j,:)))*sigInfo(j).dcOffset);		%%%% Modulated sinewave %%%%%%%%
		case 'squa'
            squarewave = sin(2*pi*sigInfo(j).Frequency*t  + (pi/180)*sigInfo(j).Phase);
            sig(j,:) = (sigInfo(j).Amplitude)*((squarewave>0)-(squarewave<0));
            sig(j,:) = sig(j,:) + (ones(1,length(sig(j,:)))*sigInfo(j).dcOffset);
        %%%% External signal file %%%%%%%%%%%%    
       otherwise
            %msg = sprintf('Signal type choices "[Cons]tant","[Ramp]","Step","[Sine]wave"\n');
            msg = ['Signal type choices "[Band] Limited Noise", "[Cons]tant",' ...
                   '"Ramp","Step",''[Squa]re wave","File", ' ...
                   '"[Sine]wave".  Your choice was: "' sigInfo(j).SigType '"'];
            error(msg);
		end
    end
end


%----------------------------------------------------------------------------
function Amps = BLimitWhiteNoise(omega,parms)
%% Function that computes the amplitudes given the frequencies omega
%% rms level and bandwidth variance.
N = length(omega);
bandwidth1 = parms(1);
bandwidth2 = parms(2);
if(bandwidth1>bandwidth2)
   bandwidth1 = parms(2);
   bandwidth2 = parms(1);
end
if(length(parms)>2);
   RandomSeed=parms(3);  
   if(RandomSeed>0)      
      randn('state',RandomSeed);
   end
end   
Amps = zeros(1,N)+i*zeros(1,N);
index = find((abs(omega)<=bandwidth2)&(abs(omega)>=bandwidth1));
num = length(index);
Amps(index) = (randn(1,num)+i*randn(1,num));