function sig = sim_genSignal(t,sigInfo,numSignals,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
%% cha@shifter.wustl.edu

% 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
for i=1:numSignals
    for j=1:numComponents
		switch lower(sigInfo(i,j).SigType(1:4))
		case 'cons'
            sig(i,j,:) = sigInfo(i,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(i,j).Frequency  2*pi*sigInfo(i,j).FrequencyHigh  37]; %% Fixed seed !!! CHA
			Amps = BLimitWhiteNoise(omega,parms);
			Amps = [Amps,fliplr(conj(Amps(2:end)))];
			S = real(ifft(Amps));
			rmsS = sqrt(sum(S.^2)/N);
			sig(i,j,:) =  sigInfo(i,j).Amplitude*S*(rms/rmsS);
	
		case 'ramp'
            ton  = sigInfo(i,j).TimeOn;  %%% * (T-T0);
            toff = sigInfo(i,j).TimeOff; %%% * (T-T0);
            iOn  = ceil(ton/dt)+2;
            iOff = ceil(toff/dt)+1;
            num = iOff - iOn + 1;
            sig(i,j,:) = zeros(1,Nt);
            sig(i,j,:) = ones(1,Nt) * sigInfo(i,j).dcOffset;
            amp = sigInfo(i,j).dcOffset;
            damp = sigInfo(i,j).Amplitude / num;
            for k=iOn:iOff
                sig(i,j,k) = amp;
                amp = amp + damp;
            end
            if (iOff < Nt)
               for k=iOff+1:Nt
                  sig(i,j,k) = amp;
               end
            end
		case 'step'
            sig(i,j,:) = zeros(1,Nt);
            sig(i,j,:) = ones(1,Nt) * sigInfo(i,j).dcOffset;
            ton  = sigInfo(i,j).TimeOn ;% * (T-T0); 
            toff = sigInfo(i,j).TimeOff;% * (T-T0);
            iOn  = ceil(ton/dt)+2;
            iOff = ceil(toff/dt)+1;
            num  = iOff-iOn+1;
            sig(i,j,iOn:iOff) = (sigInfo(i,j).Amplitude + sigInfo(i,j).dcOffset) * ones(1,num);
		%%%% Modulated sinewave %%%%%%%%
		case 'sine'
            sig(i,j,:) = (sigInfo(i,j).Amplitude) ...
                     * sin(2*pi*sigInfo(i,j).Frequency*t  + (pi/180)*sigInfo(i,j).Phase);
            sig(i,j,:) = sig(i,j,:) + (ones(1,1,length(sig(i,j,:)))*sigInfo(i,j).dcOffset);
        %%%% External signal file %%%%%%%%%%%%    
        case 'file'
            load(sigInfo(i,j).Filename);
            sig(i,j,:) = interp1(ExtTime, ExtSig,t, 'linear');        
		otherwise
            %msg = sprintf('Signal type choices "[Cons]tant","[Ramp]","Step","[Sine]wave"\n');
            msg = ['Signal type choices "[Band] Limited Noise", "[Cons]tant",' ...
                   '"[Ramp]","Step",' ...
                   '"[Sine]wave".  Your choice was: "' sigInfo(i,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);
seed=parms(3);
bandwidth1 = parms(1);
bandwidth2 = parms(2);
if(bandwidth1>bandwidth2)
   bandwidth1 = parms(2);
   bandwidth2 = parms(1);
end
if(length(parms)>2);
%   RandomSeed=parms(3);  %% Modified 10/13/03 by CHA
%   if(RandomSeed>0)          %% Bandlimited signals were all the same!!
%      randn('state',RandomSeed);
%   end
end   
Amps = zeros(1,N)+i*zeros(1,N);
index = find((abs(omega)<=bandwidth2)&(abs(omega)>=bandwidth1));
num = length(index);
if seed>0
    randn('state',seed);
end
Amps(index) = (randn(1,num)+i*randn(1,num));