function [x,t]=sum_sine(f,Tinterv,fs,phi,amp)
%% SUM_SINE  Function to produce samples of a signal that consists of a
%%           sum of sinewaves
%%
%%   Usage: 
%%     [x,t]=sum_sine(f,Tinterv,fs,phi)
%%
%%   Variables:
%%         f = vector containing sinusoidal frequencies (Hz)
%%        fs = sampling frequency (Hz)  
%%   Tinterv = length of time interval (seconds)
%%       phi = phase offset (in radians) for each sinewave [OPTIONAL]
%%            (f and phi MUST be the same length)
%%
%%         x = samples of the sinusoidal signal
%%         t = vector of sample times (for plotting)

f=f(:);				% Make sure f is a column vector

if nargin<4 			% Phase is optional
  phi=zeros(size(f));		% Set to zero if unspecified
  amp=ones(size(f));
elseif nargin<5 		% Amplitude is optional
  phi=phi(:);
  amp=ones(size(f));
else
  phi=phi(:);			% Otherwise make it a column vector
  amp=amp(:);
end

if isempty(phi); phi=zeros(size(f)); end;
if isempty(amp); amp=ones(size(f)); end;

t=0:1/fs:Tinterv;		% Vector of times
s=sin(2*pi*f*t+phi*ones(size(t)));  % Compute sinewaves
x=amp'*s;			% Sum them

return