function ad=sonify(S,D,varargin)
% sonify - Sonify sequences of durations and pitches
%
% sonify :: 
%    ([[N]] | seq(real))    ~'sequence of pitches in semitones',
%    (natural | seq(real))  ~'sequence of durations',
%    options {
%       env :: ([[N]]->[[N]])   ~'envelope shaping function',
%       transpose :: real /0    ~'add this to all pitches first';
%       bpm :: nonneg  /240     ~'number of time units per minute';
%
%       gen :: { 'sine'      ~'sine wave',
%                'square'    ~'frequency-quantised square wave',
%                'blsquare'  ~'bandlimited square wave',
%                'blsquare2' ~'alternative method',
%                'shepard'   ~'Shepard tones'
%              } / 'sine'    ~'determines waveform generated';
%       betap :: [[2]] /[2,8]  ~'Beta pdf parameters for Shepard tones'
%
%       fs  :: real    /11025   ~'intended sampling rate';
%       buffer :: noneg /0.125  ~'buffer size in seconds, [] means one per note'
%    }
% -> seq([[1,T]]).

	opts=prefs( ...
		'bpm',240,'transpose',0,    ...
		'gen','sine','betap',[2,8], ...
		'fs',11025, 'buffer',0.125, ...
		'env',[],'bl_cutoff',[],  ...
		varargin{:});

	env=opts.env;
	fm=440/opts.fs;
	dm=60*opts.fs/opts.bpm;
	tr=opts.transpose;
	num2hz=@(n)nan2x(0,fm*2.^((n+tr)/12));

	if isscalar(D), D=dm*D;
	else D=fndata(@(d)dm*d,dd(D),'compose',1); % D=dm*D with optimisation
	end
	C=fndata(@(n)(n+tr)/12,dd(S),'compose',1); % chroma values
	F=fndata(num2hz,dd(S),'compose',1);   % convert to normalised frequencies

	% cutoff for band-limited waveform generation
	if isempty(opts.bl_cutoff)
		cutoff=0.21+0.2*repeat(normalise01(cumsum(sample(stable(0,1.4,0),[1,1024]),2),2));
	else
		cutoff=opts.bl_cutoff;
	end

	switch opts.gen
		case 'square'
			% square wave with period quantized to sample period
			Y=blockdata(square,D,0,repeat(0.5),F);
		case 'blsquare'
			% cumsum within each buffer: doesn't drift but makes
			% audible clicks between buffers.
			Y=fndata(@(t)0.5*cumsum(t),blockdata(bpblit,D,0,cutoff,F));
		case 'blsquare2'
			% this version does cumsum across the sequence to remove 
			% discontinuities but requires high-pass filter to stop drift.
			[fb,fa]=butter(2,0.001,'high');
			Y=filter(0.5*fb,fa,cumsum(blockdata(bpblit,D,0,cutoff,F),2));
		case 'shepard'
			Y=0.04*shepardtone(D,C,opts.betap);
		case 'sine2'
			Y=blockdata(sine,D,0,F)./(F/0.01);
		otherwise
			Y=blockdata(sine,D,0,F);
	end

	% if supplied apply envelope to shape each note
	if nargin>2 && ~isempty(env),
		ad=fndata(env,Y,'sizecheck',1);
	else
		ad=Y;
	end

	% rebuffer to constant size
	if ~isempty(opts.buffer)
		ad=windowdata(ad,ceil(opts.buffer*opts.fs));
	end

