changeset 34:c75bb62b90a9

Imported audio synthesis tools.
author samer
date Sun, 20 Jan 2013 19:05:05 +0000
parents 5b7d90b6393a
children f1ce7876346a
files dsp/synth/@blit/blit.m dsp/synth/@blit/block.m dsp/synth/@blit/block_sr.m dsp/synth/@bpblit/block.m dsp/synth/@bpblit/block_sr.m dsp/synth/@bpblit/bpblit.m dsp/synth/@hardsync/block.m dsp/synth/@hardsync/hardsync.m dsp/synth/@msine/block.m dsp/synth/@msine/block_sr.m dsp/synth/@msine/msine.m dsp/synth/@periodic/block.m dsp/synth/@periodic/block_sr.m dsp/synth/@periodic/periodic.m dsp/synth/@sawtooth/block.m dsp/synth/@sawtooth/block_sr.m dsp/synth/@sawtooth/sawtooth.m dsp/synth/@sine/block.m dsp/synth/@sine/block_sr.m dsp/synth/@sine/sine.m dsp/synth/@square/block.m dsp/synth/@square/block_sr.m dsp/synth/@square/square.m dsp/synth/README dsp/synth/addsynth.m dsp/synth/audio_sonifier.m dsp/synth/blockdata.m dsp/synth/dd.m dsp/synth/envadr.m dsp/synth/expchirp.m dsp/synth/expsin.m dsp/synth/fir2snd.m dsp/synth/formant_synth.m dsp/synth/harmsynth.m dsp/synth/linchirp.m dsp/synth/noise.m dsp/synth/noodle.m dsp/synth/osc.m dsp/synth/playarspec.m dsp/synth/playfirspec.m dsp/synth/poly2snd.m dsp/synth/shepardtone.m dsp/synth/sonify.m dsp/synth/sonify_disc.m dsp/synth/sonify_formant.m dsp/synth/sonify_voice.m dsp/synth/spec2snd.m dsp/synth/specsynth.m dsp/synth/srdata.m dsp/synth/unbuffer_nu.m
diffstat 50 files changed, 956 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@blit/blit.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=blit
+% blit - Band-limited impulse train generator
+
+o=class(struct,'blit');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@blit/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,25 @@
+function [y,phi]=block(o,phi,T,cutoff,f)
+% block - Generate block of blit signal data
+%
+% block :: 
+%    blit       ~'blit object'
+%    0--1       ~'initial phase (in cycles)'
+%    N:natural  ~'number of samples to compute',
+%    0--1       ~'normalised frequency (1=sampling freq)',
+%    0--0.5     ~'cutoff harmonics above this frequency [0.5]',
+% -> [[1,N]]    ~'band-limited impulse train',
+%    0--1       ~'initial phase for next block'.
+
+m = max(1,2*floor(cutoff/f)-1); % number of harmonics to keep
+y = m*f*diric(2*pi*(phi+f*(0:T-1)),m)-f;
+phi = mod(phi+f*T,1);
+
+% If diric is not available then use this
+% y=zeros(1,T);
+% phase = pi*(phi+f*(0:T-1));
+% modcheck=mod(phase,pi);
+% I=modcheck>0.001 & modcheck<(pi-0.001);
+% y(~I) = f*m;
+% y(I)  = f*sin(m*phase(I))./sin(phase(I));
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@blit/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,26 @@
+function [y,phi]=block(o,phi,cutoff,f)
+% block_sr - Generate block of blit signal data
+%
+% block_sr :: 
+%    blit       ~'blit object'
+%    0--1       ~'initial phase (in cycles)'
+%    N:natural  ~'number of samples to compute',
+%    0--1       ~'normalised frequency (1=sampling freq)',
+%    0--0.5     ~'cutoff harmonics above this frequency [0.5]',
+% -> [[1,N]]    ~'band-limited impulse train',
+%    0--1       ~'initial phase for next block'.
+
+m = max(1,2*floor(max(cutoff./f))-1); % number of harmonics to keep
+u = cumsum([phi,f],2);
+y = m*f.*diric(2*pi*u(1:end-1),m)-f;
+phi = mod(u(end),1);
+
+% If diric is not available then use this
+% y=zeros(1,T);
+% phase = pi*(phi+f*(0:T-1));
+% modcheck=mod(phase,pi);
+% I=modcheck>0.001 & modcheck<(pi-0.001);
+% y(~I) = f*m;
+% y(I)  = f*sin(m*phase(I))./sin(phase(I));
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@bpblit/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,16 @@
+function [y,phi]=block(o,phi,T,cutoff,f)
+% block - Generate block of bipolar blit signal data
+%
+% block :: 
+%    bpblit     ~'bpblit object'
+%    0--1       ~'initial phase (in cycles)'
+%    N:natural  ~'number of samples to compute',
+%    0--1       ~'normalised frequency (1=sampling freq)',
+%    0--0.5     ~'cutoff harmonics above this frequency [0.5]',
+% -> [[1,N]]    ~'band-limited impulse train',
+%    0--1       ~'initial phase for next block'.
+
+m = max(2,2*floor(cutoff/max(f,eps))); % number of harmonics to keep
+y = m*f*diric(2*pi*(phi+f*(0:T-1)),m);
+phi = mod(phi+f*T,2);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@bpblit/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,18 @@
+function [y,phi]=block_sr(o,phi,cutoff,f)
+% block_sr - Generate block of blit signal data
+%
+% block_sr :: 
+%    blit       ~'blit object'
+%    0--1       ~'initial phase (in cycles)'
+%    N:natural  ~'number of samples to compute',
+%    0--1       ~'normalised frequency (1=sampling freq)',
+%    0--0.5     ~'cutoff harmonics above this frequency [0.5]',
+% -> [[1,N]]    ~'band-limited impulse train',
+%    0--1       ~'initial phase for next block'.
+
+m = max(2,2*floor(max(cutoff./f))); % number of harmonics to keep
+u = cumsum([phi,f],2);
+y = m*f.*diric(2*pi*u(1:end-1),m);
+phi = mod(u(end),2);
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@bpblit/bpblit.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=blit
+% blit - Band-limited impulse train generator
+
+o=class(struct,'bpblit');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@hardsync/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,33 @@
+function [x,phi]=block(o,phi0,N,f0,f1)
+% block - Generate block of hardsync signal data
+%
+% block :: 
+%    hardsync,
+%    0--1       ~'initial phase of master'
+%    N:natural  ~'number of samples to compute',
+%    0--1       ~'normalised frequency of master (1=sampling freq)',
+%    0--1       ~'normalised frequency of slave (1=sampling freq)',
+% -> [[1,N]]    ~'signal data',
+%    0--1       ~'initial phase for next block'.
+
+	x=zeros(1,N); 
+	l0=round(1/f0);
+	
+	j=1; T=round(l0*(1-mod(phi0,1))); 
+	phi1=f1*mod(phi0,1);
+
+	while j<=N && T<=N
+		M=T-j+1;
+		x(j:T)=phi1+f1*(0:M-1);
+		j=T+1;
+		T=T+l0;
+		phi1=0;
+	end
+	phi0=0;
+
+	M=N-j+1;
+	x(j:N)=phi1+f1*(0:M-1);
+	phi1=phi1+f1*M;
+	phi0=f0*M;
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@hardsync/hardsync.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=blit
+% blit - Band-limited impulse train generator
+
+o=class(struct,'hardsync');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@msine/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,n,f)
+% block - generate sine wave
+
+t=2*pi*f*(0:n-1);
+y=sin(repmat(2*pi*phi,[1,n]) + t);
+phi=mod(phi+f*n,1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@msine/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,f)
+% block - generate sine wave
+
+u=cumsum([phi,f],2);
+y=sin(2*pi*u(:,1:end-1));
+phi=mod(u(:,end),1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@msine/msine.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=msine
+% msine - Sinusoidal signal generator
+
+o=class(struct,'msine');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@periodic/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,4 @@
+function [y,phi]=block(o,phi,n,f)
+
+y=o.fn(phi+f*(0:n-1));
+phi=mod(phi+f*n,1);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@periodic/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,9 @@
+function [y,phi]=block_sr(o,phi,f)
+% block_ar - Generate periodic signal with sample-rate control parameters
+%
+% block_sr :: periodic, real, [[N]] -> [[N]].
+
+PH=cumsum([phi,f],2);
+y=o.fn(PH(:,1:end-1));
+phi=mod(PH(:,end),1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@periodic/periodic.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function o=periodic(fn)
+% periodic - Quantised period square wave
+
+o.fn=fn;
+o=class(o,'periodic');  
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@sawtooth/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,n,f)
+% block - generate sawtooth wave
+	
+t=round(1/f);
+y=mod(phi+(0:n-1)/t,1);
+phi=mod(phi+n/t,1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@sawtooth/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,f)
+% block - generate sawtooth wave
+	
+u=cumsum([phi,1./round(1./f)],2);
+y=mod(u(1:end-1),1);
+phi=mod(u(end),1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@sawtooth/sawtooth.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=sawtooth
+% sawtooth - Quantised period sawtooth wave
+
+o=class(struct,'sawtooth');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@sine/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,n,f)
+% block - generate sine wave
+
+w=2*pi*f;
+y=sin(2*pi*phi + w*(0:n-1));
+phi=mod(phi+f*n,1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@sine/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,f)
+% block_sr - generate sine wave
+
+u=cumsum([phi,f],2);
+y=sin(2*pi*u(1:end-1));
+phi=mod(u(end),1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@sine/sine.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=sine
+% sine - Sinusoidal signal generator
+
+o=class(struct,'sine');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@square/block.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,n,duty,f)
+% block - generate square wave
+
+t=ceil(1/f);
+y=mod(phi+(0:n-1)/t,1)<duty;
+phi=mod(phi+n/t,1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@square/block_sr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function [y,phi]=block(o,phi,duty,f)
+% block - generate square wave
+
+u=cumsum([phi,1./ceil(1./f)]);
+y=mod(u(1:end-1),1)<duty;
+phi=mod(u(end),1);
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/@square/square.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,6 @@
+function o=square
+% square - Quantised period square wave
+
+o=class(struct,'square');  % no internal information
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/README	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,94 @@
+Notes on event scheduler
+
+
+rendersig allows collections of events to be rendered as signals
+using additive composition.
+
+event collections abstraction
+
+	The partition function should take a collection of events over a time
+	domain T and a time interval and return a cell array (therefore finite
+	collection) of events that are active on the given interval, and the
+	collection of events that are still active in later intervals.
+	
+		partition :: events(T), interval(T) -> {[K]->event(T)}, events(T).
+
+	What if head events could be an infinite list? This would allow us
+	to scale the load to the available CPU time. All we need in this case
+	is a corresponding implementation of a time-limited foldl.
+	
+	renderbuf uses something like
+
+		tstep :: point(T), duration(T) -> interval(T), point(T).
+		tstep(p,d)=(i,q) :-
+			endpoints(i,p,q),
+			hasDuration(i,d).
+	
+	To render one event to a given interval, 
+
+		accum :: signal(T), event(T) -> signal(T).
+
+	signal(T) represents a finite section of signal over a given domain, ie
+	signal(T) ::= ( (T->A), interval(T) ).
+
+	accum must redefine the signal function over the intersection of the
+	event interval and the signal's domain interval, ie
+
+		accum(x,e)=y :-
+			dom(y)=dom(x),
+			forall t. t in dom(x) =>
+				y(t) = if (t in dom(e)) x(t)+e(t) else x(t)
+
+			
+
+
+Generator classes
+
+	sine
+	msine
+	blit
+	bpblit
+	hardsync
+	square
+	sawtooth
+	periodic
+
+Array signal functions
+	expsin   - params, N:natural -> [[N]]
+	envadr   - params -> [[N]] -> [[N]]
+	expchirp - params -> dur -> [[N]]
+	linchirp - params -> dur -> [[N]]
+	
+Buffered data sequences
+
+	Block based
+		blockdata
+		specsynth   :: gen  -> N:natural   -> seq [[L]] -> seq [[1,K]]
+		shepardtone :: parms -> seq natural ->              seq real   -> seq [[1,_]]
+		harmsynth   :: [[K]] -> seq natural -> seq [[K]] -> seq real   -> seq [[1,_]]
+		addsynth    :: [[K]] -> seq natural -> seq [[K]] -> seq [[K]]  -> seq [[1,_]]
+
+	srdata
+	osc :: real -> seq [[N]] -> seq [[N]]
+	
+
+Event Based renderer
+	ping   	- decaying sinusoid event
+	wtevent	- wavetable event
+	sequence - arrange events one after the other
+	tshift	- shift event time
+	at     	- set event time
+
+	rendersig
+	esonify
+	esonify1
+
+Generator class
+
+	gen : new signal
+	transformer : transform one signal into another (pattern for audio rate control)
+
+
+Sequence based renderer
+	noodle
+	sonify
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/addsynth.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,13 @@
+function y=addsynth(N,A,F,phi)
+% addsynth - Additive synthesis
+%
+% addsynth ::
+%    N:natural    ~'size of buffers to produce',
+%    seq([[K]])   ~'sequence of amplitudes for each component',
+%    seq([[K]])   ~'sequence of frequencies of sine waves to generate',
+%    [[K]]        ~'initial phases'
+% -> seq([[1,N]]) ~'sum of components'.
+
+phi=repmat_to(phi,size(F));
+y=sum(zipdata(@(a,x)spdiag(a)*x,2,A,blockdata(msine,N,phi,F)),1);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/audio_sonifier.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,67 @@
+% audio_sonifier - Create sonifier object
+%
+% audio_sonifier ::
+%    nonneg ~ 'Duration of each note',
+%    options {
+%       transpose/0    :: int    ~'pitch transposition in semitones';
+%       gen    /'sine' :: {'sine','blsquare'};
+%       attack /5      :: nonneg ~'attack time in ms';
+%       decay  /200    :: nonneg ~'decay time in ms';
+%       release/5      :: nonneg ~'release time in ms';
+%       cutoff /0.4    :: nonneg ~'low pass cut-off as norm freq';
+%       fs     /11025  :: nonneg ~'sample rate'
+%    }
+% => sonifier.
+%
+% sonifier ::= 
+%   seq(natural) ~'thing to sonify', 
+%   (void=>void) ~'action on finish' 
+% => player.
+%
+% player ::= struct {
+%    start     :: void => void;    
+%    startq    :: (nonneg,nonneg) => void;    
+%    stop      :: void => void;
+%    isrunning :: void => bool;
+%    atend     :: void => bool;
+%    rewind    :: void => void;
+%    dispose   :: void => void
+% }.
+
+function f=audio_sonifier(device,dur,varargin)
+	audio_info = prefs( ...
+				'transpose',0,'gen','sine','fs',11025, 'cutoff',0.4, ...
+				'attack',5,'decay',200,'release',5,varargin{:});
+
+	if isempty(device)
+		audio_info.lineout = lineout(1,audio_info.fs);
+	else
+		audio_info.lineout = device;
+		audio_info.fs = rate(device);
+	end
+
+	f=@(seq,onfinish)sonify_audio(audio_info,dur,seq,onfinish);
+end
+
+function player=sonify_audio(audio,dur,seq,onfinish)
+	disp('preparing audio signal.');
+	sched  = playaudio_async( ...
+		sonify(seq,dur,audio_opts(audio)), audio.lineout, ...
+		'defer',1,'its',2^31,'onfinish',@(a)onfinish());
+
+	player=cpfields({'stop','isrunning','atend','rewind','dispose'},...
+		sched, struct(...
+			'start',  @()sched.startat(nows+1), ...
+			'startq', @(q,o)sched.startat(o+quant(q,nows))));
+%			'wait',   sig.wait));
+
+	disp('audio playback ready to go.');
+end
+
+function opts=audio_opts(A)
+	c=@(t)round(A.fs*t/1000);
+	opts=prefs('bpm',60, 'transpose',A.transpose-69, ...
+		'env',envadr(c(A.attack),c(A.decay),c(A.release)), ...
+		'gen',A.gen,'bl_cutoff',A.cutoff,'fs',A.fs,'buffer',0.125);
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/blockdata.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,41 @@
+function x=blockdata(gen,dur,s0,varargin)
+% blockdata - Signal as a sequence of blocks from signal generator object
+%
+% blockdata :: 
+%    gen(A,B{1:N}) ~'signal generator object with state A and N params B',
+%    seq(natural)  ~'number of samples to generate in each block',
+%    A             ~'initial state',
+%    seq(B{1})     ~'sequence of 1st param values',
+%    ...           
+%    seq(B{N})     ~'sequence of Nth param values',
+% -> seq([[1,_]])  ~'heterogenous sequence of arrays'.
+%
+% The main point of this construct is to maintain state
+% between calls to the block method of the signal generator, and
+% also to supply the generator with a sequence of parameter
+% values.
+
+	nargs=length(varargin)+1;
+	params=cellmap(@dd,varargin); % make sure all input sequences actually are
+	perm=[nargs+1,1,2:nargs];
+	x=zipaccum(@blockgen,s0,[{dd(dur)},params]); 
+
+	% !! NOTE: it would also be possible to NOT maintain state between calls,
+	% so that each block would start from the same initial state
+	% x=zipwith(@blockgen_nostate,dd(dur),params{:}); 
+
+	function [y,state]=blockgen(varargin)
+		% this is weird - for some reason Matlab crashes if I call block
+		% directly from here, but is ok if I go via the non-nested bblock
+		[y,state]=bblock(gen,varargin{perm});
+	end
+	
+%	function y=blockgen_nostate(theta)
+%		y=block(gen,s0,ceil(theta{1}),theta{2:end});
+%	end
+end
+
+function [y,s]=bblock(varargin) 
+	[y,s]=block(varargin{:}); 
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/dd.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,14 @@
+function x=dd(x)
+% dd - convert input sequence or array to sequence
+%
+% dd :: [[N,M]] -> seq([[N]]).
+% dd :: [[N]] -> seq(real).
+% dd :: seq([[N,M]]) -> seq([[N]]).
+
+if isseq(x)
+	if size(x,2)>1, x=concat(map(@slices,x)); end
+else
+	if isscalar(x),     x=repeat(x); 
+	else                x=slices(x);
+	end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/envadr.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,30 @@
+function ef=envadr(att,dec,rel)
+% envadr - Make an attack-decay-release envelope function
+%
+% envadr :: 
+%    natural ~'attack time in samples',
+%    real    ~'decay time constant in samples',
+%    natural ~'release time in seconds'
+% -> ([[N]]->[[N]]) ~'function to apply envelope to given signal'
+%
+% Attack and release are linear, decay is exponential.
+
+	e=expenv(600000,att,dec);
+	r=linspace(0,1,rel);
+	ef=@shape;
+
+	function x=shape(x)
+		L=length(x);
+		if L>0,
+			x=x.*e(1:L);
+			I=L - (0:length(r)-1);
+			x(I)=x(I).*r;
+		end
+	end
+end
+
+function x=expenv(t,att,decay)
+	x=exp(-(0:t-1)/decay);
+	x(1:att)=((0:att-1)/att).*x(1:att);
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/expchirp.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,11 @@
+function y=expchirp(fs,f0,f1,t)
+% expchirp - exponential chirp
+%
+% expchirp ::
+%    F:real   ~'sampling frequency',
+%    real     ~'initial frequency',
+%    real     ~'final frequency',
+%    D:nonneg ~'duration'
+% -> [[1,D*F]].
+
+y=chirp(0:(1/fs):t,f0,t,f1,'logarithmic',-90);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/expsin.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,22 @@
+function x=expsin(n,k,l,phi)
+% expsin - return real part of complex exponential.
+%
+% expsin :: 
+%    N:natural ~'length of signal to return',
+%    real      ~'normalised frequency 1=sample rate',
+%    real      ~'decay rate, 1/time constant'
+% -> [[N]].
+%
+% expsin :: 
+%    N:natural ~'length of signal to return',
+%    real      ~'normalised frequency 1=sample rate',
+%    real      ~'decay rate, 1/time constant'
+%    real      ~'initial phase'
+% -> [[N]].
+% 
+% If initial phase not given, defaults to -pi/2 (so initial
+% value is zero with +ve slope).
+
+if nargin<4, phi=pi/2; end
+t=0:n-1;
+x=real(exp((2*pi*k*i-l)*t-phi*i));
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/fir2snd.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,26 @@
+function X=fir2snd(A,u)
+% fir2snd - Generate sound by applying FIR filters to a noise grain
+%
+% fir2snd ::
+%    [[N,L]]              ~'array of M FIR filters of length N',
+%    ([[M]] | seq([[M]])) ~'noise grain or sequence of noise grains'
+% -> [[abs(M-N)+1,L]]       ~'columns of filtered noise'.
+
+n=size(A,1);
+l=size(A,2);
+m=size(u,1);
+
+% want to extract the valid portion of the convolution only
+if m>n, seg=n:m;
+else seg=m:n; end
+
+if ~isseq(u), u=repeat(u); end
+
+X=zeros(length(seg),l);
+for k=1:l
+	x=conv(A(:,k),head(u));
+	X(:,k)=x(seg);
+	u=next(u);
+end
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/formant_synth.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,13 @@
+% formant_synth - synthesise audio from format-based description
+function Y=formant_synth(block_size,q,amp,f0,fx)
+Y = amp*dynfilter( ...
+		map(@f2tf, fx), ...
+		blockdata(blit,block_size,0,0.45,f0) ...
+	);
+
+	function tf=f2tf(f)
+		zf=zeros(size(f));
+		[tf{1},tf{2}]=zp2tf([zf;zf],fq2poles(0.5*f,q),min(f));
+	end
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/harmsynth.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,13 @@
+function y=harmsynth(N,A,F0,phi)
+% harmsynth - Additive synthesis with harmonic frequencies
+%
+% harmsynth ::
+%    N:natural    ~'size of buffers to produce',
+%    seq([[K]])   ~'sequence of amplitudes for each component',
+%    seq(nonneg)  ~'sequence of fundamental frequencies',
+%    [[K]]        ~'initial phases'
+% -> seq([[1,N]]) ~'sum of components'.
+%
+% Note: the number of components must remain constant
+
+y=addsynth(N,A,F0*(1:size(A,1))',phi);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/linchirp.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,12 @@
+function y=linchirp(fs,f0,f1,t)
+% linchirp - Linear chirp
+%
+% linchirp ::
+%    F:real ~'sampling frequency',
+%    real ~'initial frequency',
+%    real ~'final frequency',
+%    D:nonneg ~'duration'
+% -> [[1,D*F]].
+
+y=chirp(0:(1/fs):t,f0,t,f1,'linear',-90);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/noise.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,10 @@
+% noise - white noise generator
+%
+% noise :: 
+%    D:natural ~'maximum duration to expect',
+% -> (   nonneg  ~'sampling frequency', 
+%        N:1..D  ~'duration' 
+%     -> [[1,N]] ~'some noise') ~'function to generate white noise'.
+function gen=noise(maxdur)
+	buf=randn(1,maxdur);
+	gen=@(fs,dur)buf(1:dur);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/noodle.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,35 @@
+function y=noodle(Pitch,Dur,N,varargin)
+% noodle - noodle with arbitrary note durations and meter
+%
+% noodle ::
+%    seq(natural) ~'sequence of top-line pitches in semitones',
+%    seq(nonneg)  ~'sequence of durations',
+%    natural      ~'bass note duration multiplier',
+%    options {
+%        decay1 :: real    /8000  ~'time constant for melody decay';
+%        decay2 :: real    /1600  ~'time constant for bass decay';
+%        period :: real    /1600  ~'samples per metrical unit';
+%        buffer :: nonneg  /0.4   ~'duration of buffers in seconds'
+%    }
+% -> seq([[1,N]]) ~'sequence of audio buffers'.
+
+opts=prefs('bpm',120,'buffer',0.4, ...
+	'decay1',8000,'decay2',1600, ...
+	'attack1', 16, 'attack2', 16, ...
+	'release1', 64, 'release2', 64, ...
+	'transpose',-5, ...
+	'fs',11025, ...
+	'betap1',[2.2,8], ...
+	'betap2',[1.1,12], ...
+	varargin{:});
+
+env = { ...
+	envadr(opts.attack1,opts.decay1,opts.release1), ...
+	envadr(opts.attack2,opts.decay2,opts.release1) };
+
+s1=sonify(Pitch,Dur,opts,'env',env{1},'betap',opts.betap1);
+s2=sonify(Pitch-12,N,opts,'env',env{2},'betap',opts.betap2);
+y=0.27*s1+0.50*s2;
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/osc.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,13 @@
+function [y,phi]=osc(f,phi0)
+% osc - sequence of oscillator output
+%
+% osc :: seq [[N]] ~'persample normalised frequencies' -> seq [[N]] ~'sinusoidal data'.
+%
+% osc :: 
+%    seq([[N]]) ~'persample normalised frequencies', 
+%    real       ~'initial phase' 
+% -> seq([[N]]) ~'sinusoidal data',
+%    seq([[N]]) ~'phase data'.
+
+if nargin<2, phi0=zeros(size(f,1),1); end
+y=srdata(sine,phi0,f);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/playarspec.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,10 @@
+function X=playarspec(S,h)
+% PLAYARSPEC - Return (or play) sound for a given power spectrum
+%
+% playarspec :: (S:spectrum,h:window)->x:signal
+% S can be an array of spectra
+% h is window or `grain shape' used to synthesis the sound
+% (What about overlap? See stereo reconstrunction functions)
+
+x=poly2snd(ac2poly(spec2ac(S)),h.*randn(size(h)));
+if nargout<1, sound(flatten(x),11025); else X=x; end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/playfirspec.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,7 @@
+function X=playfirspec(A,h)
+% PLAYFIRSPEC - Synthesis spectrum using FIR filters
+%
+% playfirspec :: [[N,M]], [[L]] -> [[T]].
+
+x=fir2snd(spec2fir(A),h.*randn(size(h)));
+if nargout<1, sound(flatten(x),11025); else X=x; end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/poly2snd.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,18 @@
+function X=poly2snd(A,u)
+% POLY2SND  - Convert array of filter polynomial coeffs to one long sound by filtering noise
+% Filter state is preserved as the function processes a number of filters
+%
+% poly2snd :: 
+% 	(A:array 1..N of Poly, signal~'noise to filter')
+%	->signal~'N concatentated sounds'
+%
+% Each ROW of A contains coefficients for one filter
+
+Z=[];
+for k=1:size(A,1)
+	[x,Z]=filter(1,A(k,:),u,Z);
+	if any(isnan(Z)), Z=[]; end;
+	X(:,k)=0.999*x/max(abs(x));
+end
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/shepardtone.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,23 @@
+function y=shepardtone(N,C,prof,H)
+% shepardtone - Shepard tone by additive synthesis
+%
+% shepardtone ::
+%    N:natural    ~'size of buffers to produce',
+%    seq(real)    ~'sequence of chroma values, 0--1 is one octave',
+% -> seq([[1,N]]) ~'sum of components'.
+%
+% Note: the number of components must remain constant
+
+if nargin<3, 
+	prof=@(z)betapdf(2*z,2,9); 
+elseif isvector(prof),
+	a=prof(1); b=prof(2);
+	prof=@(z)betapdf(2*z,a,b);
+end
+
+fc=0.25;
+if nargin<4, H=(-10:1)'; else H=H(:); end
+F=fndata(@(c)fc*2.^(H+mod(c,1)),C);  % frequencies
+A=fndata(prof,F); % component amplitudes
+y=addsynth(N,A,F,0);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/sonify.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,83 @@
+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
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/sonify_disc.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,39 @@
+function Y=sonify_disc(S,D,Gen,varargin)
+% sonify - Sonify sequences of durations and sound classes
+%
+% sonify :: 
+%    ([[N]] | seq(real))      ~'sequence of pitches in semitones',
+%    (natural | seq(real))   ~'sequence of durations',
+%    options {
+%       bpm :: nonneg  /240     ~'number of time units per minute';
+%       fs  :: real    /11025   ~'intended sampling rate';
+%       buffer :: noneg /0.125  ~'buffer size in seconds, [] means one per note'
+%    }
+% -> seq([[1,T]]).
+
+	opts=prefs('env',[],'bpm',240,'fs',11025,'buffer',0.125, varargin{:});
+
+	sample_rate=opts.fs;
+	dm=60*opts.fs/opts.bpm;
+
+	if isscalar(D), D=dd(dm*D);
+	else D=fndata(@(d)dm*d,dd(D),'compose',1); % D=dm*D with optimisation
+	end
+
+	Y=zipdata(@select_gen,2,D,dd(S));
+
+	% if supplied apply envelope to shape each note
+	if ~isempty(opts.env),
+		Y=fndata(opts.env,Y,'sizecheck',1);
+	end
+
+	% rebuffer to constant size
+	if ~isempty(opts.buffer)
+		Y=windowdata(Y,ceil(opts.buffer*opts.fs));
+	end
+
+	function y=select_gen(dur,k)
+		y=feval(Gen{k},sample_rate,floor(dur));
+	end
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/sonify_formant.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,23 @@
+% sonify_formant - sonify using formant synthesis
+function Y=sonify_formant(Controllers,S,varargin)
+
+	opts=prefs('bpm',240,'fs',11025,'buffer',0.125,'q',0.05,'cfilter',[],varargin{:});
+
+	CC = cellmap(@(f)f/opts.fs,Controllers);
+	f0 = window(map(@getf0,S));
+	fx = window(map(@getfx,S));
+
+	if ~isempty(opts.cfilter)
+		f0 = filter(opts.cfilter{1},opts.cfilter{2},f0);
+		fx = filter(opts.cfilter{1},opts.cfilter{2},fx,[],2);
+	end
+	Y  = formant_synth(256,opts.q,0.2,f0,fx);
+
+	% rebuffer to constant size
+	if ~isempty(opts.buffer)
+		Y=windowdata(Y,ceil(opts.buffer*opts.fs));
+	end
+
+	function z=getf0(s),z=CC{s}(:,1)'; end
+	function z=getfx(s),z=CC{s}(:,2:end)'; end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/sonify_voice.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,8 @@
+function Y=sonify_voice(S)
+
+Filt= decode_filter(S);
+Fun = decode_fundamental(S);
+
+Y = filter(Filt,source(Fun));
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/spec2snd.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,13 @@
+function Y=spec2snd(X,frame,hop,varargin)
+% spec2snd - Sound from sequence of power spectra
+%
+% spec2snd :: 
+% 	  seq([[M,L]])	~'power spectra in columns',
+%	  natural   	~'frame size for resynthesis',
+% 	  natural	   ~'hop size for unbuffer'
+% -> seq([[M]])	~'output sequence'.
+
+
+src=rndseq(gaussian,frame+size(X,1)-1);
+Y=unbuffer(fir2snd(spec2fir(X),src),hop);
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/specsynth.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,18 @@
+function Y=specsynth(S,gen,N,M)
+% specsynth - Spectral synthesis by FIR filtering white noise 
+%
+% specsynth ::
+%    seq([[L]])    ~'sequence of L-point power spectra',
+%    (N:natural=>[[1,N]]) ~'generator of noise buffers',
+%    N:natural     ~'size of buffers of noise signals to take'
+%    M:natural     ~'frame overlap'
+% -> seq([[1,M]])  ~'sequence of signal buffers'.
+
+% NOTE: this should really take gen to be seq natural -> seq [[1,_]]
+% eg
+% gen = @(n)rndzip(sampler(gaussian),n,getrndstate)
+% this should be  the basic type of a buffered signal generator
+
+if nargin<4, M=N; end
+U=rndmap(gen,repeat(N),getrndstate);
+Y=unbuffer(zipwith(@conv,map(@(s)spec2fir(sqrt(s)),S),U),M);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/srdata.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,29 @@
+function x=srdata(gen,s0,varargin)
+% srdata - Signal as a sequence of blocks from signal generator object
+%
+% srdata :: 
+%    gen(A,B{1:N})  ~'signal generator object',
+%    A              ~'initial state',
+%    seq(B{1}       ~'sequence of 1st param values',
+%    ...   
+%    seq(B{N}       ~'sequence of Nst param values',
+% -> seq([[1,_]]).
+%
+% The main point of this construct is to maintain state
+% between calls to the block method of the signal generator, and
+% also to supply the generator with a sequence of parameter
+% values.
+
+	x=mapaccum(@blockgen,s0,zip(varargin{:}),'sizecheck',1); 
+
+	function [y,state]=blockgen(theta,state)
+		% this is weird - for some reason Matlab crashes if I call block
+		% directly from here, but is ok if I go via the non-nested bblock
+		[y,state]=bblock(gen,state,theta);
+	end
+end
+
+function [y,s]=bblock(gen,s,theta)
+	[y,s]=block_sr(gen,s,theta{:});
+end
+	
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/synth/unbuffer_nu.m	Sun Jan 20 19:05:05 2013 +0000
@@ -0,0 +1,35 @@
+% unbuffer_nu - Non-uniform hop overlap and add
+%
+% unbuffer_nu :: 
+%    seq([[1,_]])  ~'sequence of overlapping frames',
+%    seq(natural)  ~'sequence of hop sizes'
+% -> seq([[1,_]])  ~'sequence of de-overlapped frames'.
+
+function Y=unbuffer_nu(X,hop)
+	Y=zipaccum(@olap,{hop,X},[]);
+
+	function [y,s1]=olap(hop,x,s)
+		ls=length(s);
+		lx=length(x);
+		if lx>=hop
+			if ls>=hop
+				% NB: this will fail if ls>lx, but this shouldn't happen since ls<=lx-hop
+				y=(x(1:hop)+s(1:hop))';
+				s1=[s(hop+1:ls)+x(hop+1:ls);x(ls+1:end)];
+			else
+				y=[x(1:ls)+s;x(ls+1:hop)]';
+				s1=x(hop+1:end);
+			end
+		else
+			if ls>=hop
+				y=[s(1:lx)+x;s(lx+1:hop)]';
+				s1=s(hop+1:end);
+			else
+				y=zeros(1,hop);
+				y(1:ls)=y(1:ls)+s';
+				y(1:lx)=y(1:lx)+x';
+			end
+		end
+	end
+end
+