function out = trianglewave(t,f,shape)

%TRIANGLEWAVE generates triangle or sawtooth waves.
% trianglewave = TRIANGLEWAVE(T,FREQ,SHAPE)
%
% TRIANGLEWAVE generates a triangle wave of frequency FREQ for the input
% time vector T. If SHAPE is specified, it can also generate SAWTOOTH
% waves. The waves it generates has peaks of +1 and -1.
%
% TRIANGLEWAVE(T) generates a triangle wave with a default frequency of
% 40 Herz.
% TRIANGLEWAVE(T,FREQ) uses a frequency specified by FREQ. A FREQ value
% less than the Nyquist frequency is recommended.
% TRIANGLEWAVE(T,FREQ,SHAPE) generates a wave with the specific SHAPE.
% SHAPE could be either 'SAWTOOTH0' which is a inverse saw wave, or
% 'SAWTOOTH1' which is the conventional sharply drop saw wave. The
% default SHAPE is 'TRIANGLE'.
%
% The SAWTOOTH function from signal processing toolbox does the same job,
% using a different algorithm. SAWTOOTH is sometimes inaccurate.
%
% Example:
% t = 0:1/200:0.5;
% w = trianglewave(t,30);
% figure; plot(w);
% v = sawtooth(2*pi*30*t,0.5); % comparing with SAWTOOTH:
% figure; plot(v);
%
% See also SQUARE, SAWTOOTH.

% Code by Siyi Deng, sdeng@uci.edu
% HNL, COG SCI, UCI.
% 08-16-2006; Ver 1.0;
% 08-17-2006; Ver 1.01; efficency improvement.
% 08-21-2006; Ver 1.02; vector dimension check.

if nargin == 0, t = linspace(0,1/40,16000/40); end;
if nargin < 2, f = 40; end;
if nargin < 3, shape = 't'; end;

if ~isvector(t), error('ERROR: Time array dimension mismatch.'); end

ts = length(t);
fs = ts/(t(end) - t(1)); % sampling frequency;

if f>fs, error('ERROR: Input frequency exceeds sampling frequency.'); end

ds = fs/f; % samples per duration;

acme = floor(1:ds:ts);
intv = acme(2:end) - acme(1:end-1)-1;
intv = [intv intv(end)];
intv(intv<1) = 1;

lin = 1:1:ts;
ind = ceil(lin./ds);
ind(ind>length(acme)) = 1;

switch shape
case {'s0','S0','sawtooth0','SAWTOOTH0','Sawtooth0'}
out = (lin-acme(ind))./intv(ind); % height normalization by period;
out = out.*2-1;
case {'s1','S1','sawtooth1','SAWTOOTH1','Sawtooth1'}
out = (lin-acme(ind))./intv(ind) ;
out = (1-out).*2-1;
otherwise
apex = floor(ds/2:ds:ts);
out = (lin-apex(ind));

% positive part normalization;
pIntv = acme(2:end) - apex(1:end-1)-1;
pIntv = [pIntv pIntv(end)];
pIntv(pIntv<1) = 1;
pInd = ind;
pInd(out<0) = 1;
pos = out;
pos(pos<0) = 0;
pos = pos./pIntv(pInd);

% negtive part normalization;
nIntv = apex(1:end-1)-acme(1:end-1)+1;
nIntv = [nIntv nIntv(end)];
nIntv(nIntv<1) = 1;
nInd = ind;
nInd(out>0) = 1;
neg = out;
neg(neg>0) = 0;
neg = neg./nIntv(nInd);

% final output normalization;
out = 1-abs(pos + neg);
out = out.*2 -1;
out(out>1) =1;
out(out<-1) = -1;
end

if size(t,2) == 1, out = out'; end