initial commit

dev/externals/reverb/CMakeLists.txt

@@ -0,0 +1,31 @@
+if(APPLE)
+  set(CMAKE_C_FLAGS "-O2 -Wall -mmacosx-version-min=10.8  -Dunix -fomit-frame-pointer -funroll-loops")
+  set(CMAKE_CXX_FLAGS "-O2 -Wall -mmacosx-version-min=10.8  -Dunix -fomit-frame-pointer -funroll-loops -std=c++11 -stdlib=libc++")
+else()
+  if(MINGW)
+    set(CMAKE_C_FLAGS "-O3 -DWIN32 -D_WIN32 -fomit-frame-pointer  -funroll-loops")
+    set(CMAKE_CXX_FLAGS "-O3 -DWIN32 -D_WIN32 -fomit-frame-pointer  -funroll-loops -static-libgcc -static-libstdc++")
+  else()
+    set(CMAKE_C_FLAGS "-O3 -Wall -Dlinux -Dunix -fomit-frame-pointer -funroll-loops")
+  endif()
+endif()
+
+link_directories(../lib)
+
+include_directories(../../../include ../include)
+  
+add_executable(rmverb rmverb.cpp reverberator.cpp wavetable.cpp)
+target_link_libraries(rmverb portsf sfsys ${EXTRA_LIBRARIES})
+my_install(rmverb)
+
+add_executable(reverb reverb.cpp reverberator.cpp wavetable.cpp)
+target_link_libraries(reverb portsf sfsys ${EXTRA_LIBRARIES})
+my_install(reverb)
+
+add_executable(rmresp roomresp.cpp)
+target_link_libraries(rmresp portsf sfsys ${EXTRA_LIBRARIES})
+my_install(rmresp)
+
+add_executable(tapdelay tdelaymain.cpp reverberator.cpp wavetable.cpp)
+target_link_libraries(tapdelay portsf sfsys ${EXTRA_LIBRARIES})
+my_install(tapdelay)

dev/externals/reverb/LARGERM.TXT

@@ -0,0 +1,7 @@
+1.0               # Input gain for data file (float)
+0                 # Open Path = 0 Closed Path = 1 (int)
+43    31    19       # Room Dimensions (X,Y,Z) (float)
+0.99              # Reflectivity of walls (float)
+5                 # Number of rooms per dimension, 1 or greater
+20    15    2      # Listener coord (X,Y,Z) (float)
+7    15    2   0 # Source coord and time at coord (X,Y,Z,t) (float)

dev/externals/reverb/MEDIUMRM.TXT

@@ -0,0 +1,7 @@
+1.0               # Input gain for data file (float)
+0                 # Open Path = 0 Closed Path = 1 (int)
+21   16    5       # Room Dimensions (X,Y,Z) (float)
+0.9              # Reflectivity of walls (float)
+5                 # Number of rooms per dimension, 1 or greater
+12    8    2      # Listener coord (X,Y,Z) (float)
+5    8    2   0 # Source coord and time at coord (X,Y,Z,t) (float)

dev/externals/reverb/POS.TXT

@@ -0,0 +1,7 @@
+1.0               # Input gain for data file (float)
+0                 # Open Path = 0 Closed Path = 1 (int)
+20    10    6       # Room Dimensions (X,Y,Z) (float)
+0.5              # Reflectivity of walls (float)
+5                 # Number of rooms per dimension, 1 or greater
+15    5    2      # Listener coord (X,Y,Z) (float)
+5    4    2   0 # Source coord and time at coord (X,Y,Z,t) (float)

dev/externals/reverb/allpass.c

@@ -0,0 +1,252 @@
+/*
+ * Copyright (c) 1983-2013  Composers Desktop Project Ltd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+/*	allpass.c
+ *	allpass filter
+ *	A. Bentley, Composers' Desktop Project, Version Nov 1987
+ */
+
+
+/*
+ *	$Log: allpass.c%v $
+ * Revision 3.3  1994/10/31  22:08:15  martin
+ * first rcs version
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>	
+#include <math.h>
+#include <sfsys.h>
+#include <cdplib.h>
+static char *sfsys_h_rcsid = SFSYS_H_RCSID;
+static char *cdplib_h_rcsid = CDPLIB_H_RCSID;
+
+void usage();
+void allpass_long(long,long,long,long,const short*,short*,long*,short *);
+
+#define BUFLEN	(32768)
+#define FP_ONE  (65536.0)
+#define LOG_PS  (16)
+
+PROGRAM_NUMBER(0xf93c2700);
+
+void
+main(argc,argv)
+int argc;
+char *argv[];
+{
+	int ifd, ofd;
+	int sampsize;
+	int quickflag = 1;	/* default is fast fixed point calc */
+	register int rv;
+	float	ip, op;
+#ifdef ORIGINAL_VERSION
+	long	lip, lop;
+	long	i;
+#endif
+	float 	*delbuf1;
+	float   *delbuf2;
+	long	rdsamps;
+	long	*ldelbuf1;
+	short	*sdelbuf2;
+	short	*sinbuf;
+	short	*soutbuf;
+	double  delay; 
+	float	gain, prescale = 1.0f;
+	long	lgain, lprescale;
+	long 	delsamps;
+	long	dbp1 = 0, dbp2 = 0;
+	long	chans;
+	long	isr;
+	double	sr;
+
+	if(sflinit("allpass")){
+		sfperror("Allpass: initialisation");
+		exit(1);
+	}
+    	while(argv[1] != 0 && argv[1][0] == '-') {
+	      switch(argv[1][1]) {
+	      case 'f':	
+		 quickflag--;
+		 argv++;
+		 argc--;
+		 break;
+	      default:
+	         usage();
+	         break;
+	      }
+	}
+	if(argc < 5 || argc > 6) usage();
+
+	if(argv[1] != 0) {
+	      	if((ifd = sndopen(argv[1])) < 0) {
+		   fprintf(stderr,"Allpass: cannot open input file %s\n", argv[1]);
+		   exit(1);
+		}
+	}
+	if(sndgetprop(ifd,"sample rate",(char *) &isr, sizeof(long)) < 0){
+		fprintf(stderr,
+		  "\nAllpass: failed to read sample rate for file %s",argv[1]);
+		exit(1);
+	}
+	sr = (double) (isr);
+	if(sndgetprop(ifd,"channels",(char *) &chans, sizeof(long)) < 0){
+		fprintf(stderr,
+		  "\nAllpass: failed to read channel data for file %s", argv[1]);
+		exit(1);
+	}
+	if(chans != 1){
+		fprintf(stderr,"\nAllpass works only on mono files");
+		exit(1);
+	}
+	if(sndgetprop(ifd,"sample type",(char *) &sampsize, sizeof(long)) < 0){
+		fprintf(stderr,
+		  "\nAllpass: failed to read sample type for file %s", argv[1]);
+		exit(1);
+	}
+	/* get command line args */
+	delay = abexpr(argv[3],sr);
+	if(delay < 0.0){
+		fprintf(stderr,"\nAllpass: invalid delay parameter");
+		exit(1);
+	}
+	gain = (float)abexpr(argv[4],sr);
+	if((gain < -1.0) || (gain > 1.0)){
+		fprintf(stderr,"\nAllpass: gain out of range");
+		exit(1);
+	}
+
+	if((ofd = sndcreat(argv[2],-1,sampsize)) < 0) {
+		  fprintf(stderr,"Cannot open output file %s\n",argv[2]);
+		  exit(1);
+	}
+
+	lgain = (long) (gain * FP_ONE);
+	if(argc > 5) {
+		prescale = (float)abexpr(argv[5],sr);
+		if((prescale < -1.0) || (prescale > 1.0))
+		   fprintf(stderr,
+			"Allpass: warning - prescale exceeds +/-1\n");
+	}
+	lprescale = (long) (prescale * FP_ONE);		
+	/* allocate buffer for delay line */
+	delsamps = (long) (delay * sr);
+	if(quickflag) {
+	   if(( sinbuf = (short *) malloc(BUFLEN * sizeof(short))) == NULL) {
+		fprintf(stderr,"\nNot enough memory for input buffer");
+		exit(1);
+	   }
+	   if(( soutbuf = (short *) malloc(BUFLEN * sizeof(short))) == NULL) {
+		fprintf(stderr,"\nNot enough memory for output buffer");
+		exit(1);
+	   }
+	   if(( ldelbuf1 = (long *) calloc(delsamps, sizeof(long))) == NULL) {
+		fprintf(stderr,"\nNot enough memory for delay buffer 1");
+		exit(1);
+	   }
+	   if(( sdelbuf2 = (short *) calloc(delsamps, sizeof(short))) == NULL) {
+		fprintf(stderr,"\nNot enough memory for delay buffer 2");
+		exit(1);
+	   }
+	}else{
+	   if(( delbuf1 = (float *) calloc(delsamps, sizeof(float))) == NULL) {
+		fprintf(stderr,"\nNot enough memory for delay buffer 1");
+		exit(1);
+	   }
+	   if(( delbuf2 = (float *) calloc(delsamps, sizeof(float))) == NULL) {
+		fprintf(stderr,"\nNot enough memory for delay buffer 2");
+		exit(1);
+	   }
+	}
+	stopwatch(1);
+    
+	for(;;){
+		if((rv = fgetfloat(&ip,ifd)) < 0){
+			fprintf(stderr,"\nError in reading file"); 
+			exit(1);
+		}
+		if(!rv) break;
+			ip *= prescale;
+		op = (-gain) * ip;
+		op += delbuf1[dbp1];
+		op += gain * delbuf2[dbp2];
+			delbuf1[dbp1++] = ip;
+		delbuf2[dbp2++] = op;
+
+		if(dbp1 >= delsamps) dbp1 = 0;
+		if(dbp2 >= delsamps) dbp2 = 0;
+
+   		if(fputfloat(&op,ofd) < 1){
+			fprintf(stderr,"\nError writing to output file"); 
+			exit(1);
+		}	
+    }
+	stopwatch(0);
+
+	sndclose(ifd);
+	if(sndputprop(ofd,"sample rate", (char *) &isr, sizeof(long)) < 0){
+		fprintf(stderr,"\nAllpass: failed to write sample rate");
+	}
+	if(sndputprop(ofd,"channels",(char *) &chans, sizeof(long)) < 0){
+		fprintf(stderr,"\nAllpass: failed to write channel data");
+	}
+	sndclose(ofd);
+}
+
+void
+usage()
+{
+	fprintf(stderr,"\nAllpass - A Groucho Program: $Revision: 3.3 $");
+	fprintf(stderr,"\n%s\n",
+	"usage: allpass [-f] infile outfile delay gain [prescale]");
+	exit(0);
+}
+
+/*RWD: my allpass functions */
+//TODO: optimize using ptr arithmetic
+void allpass_long(long prescale, 
+				   long buflen, 
+				   long l_gain,
+				   long l_delsamps,
+				   const short *s_inbuf,
+				   short *s_outbuf,
+				   long *l_delbuf1,
+				   short *s_delbuf2)
+{
+	int i; 
+	static long	db_p1 = 0;
+	static long db_p2 = 0;
+	//do allpass on block of samples
+	for(i = 0; i < buflen; i++){
+		   long lip,lop;
+		   lip = s_inbuf[i] * prescale;
+		   lop = (-l_gain) * (lip >> LOG_PS);
+		   lop += l_delbuf1[db_p1];
+		   lop += l_gain * s_delbuf2[db_p2];
+
+		   l_delbuf1[db_p1++] = lip;
+		   s_outbuf[i] = s_delbuf2[db_p2++] = (short) (lop >> LOG_PS);
+
+		   if(db_p1 >= l_delsamps) db_p1 = 0;
+		   if(db_p2 >= l_delsamps) db_p2 = 0;
+    }
+	   //block done
+	//return &soutbuf[0];			//if we return anything?
+}

dev/externals/reverb/delay.c

@@ -0,0 +1,546 @@
+/*
+ * Copyright (c) 1983-2013  Composers Desktop Project Ltd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+/* delay.c  
+ * generalised delay line featuring
+ * feedthrough, feedforward and feedback
+ * A. Bentley, Composers' Desktop Project, Nov 1987.
+ * Revised and extended to stereo and command-line mode 
+ * by Richard Dobson Oct 1993.
+ */ 
+
+static char *rcsid = "$Id: delay.c%v 3.3 1994/10/31 22:30:12 martin Exp $";
+/*
+ *	$Log: delay.c%v $
+ * Revision 3.3  1994/10/31  22:30:12  martin
+ * first rcs version
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "sfsys.h"
+#include "cdplib.h"
+static char *sfsys_h_rcsid = SFSYS_H_RCSID;
+static char *cdplib_h_rcsid = CDPLIB_H_RCSID;
+
+#define remainder	my_remainder
+
+#define BUFLEN	  (4096)
+#define	DELBUFSIZ (176400)	/* allows 2sec delay in stereo */
+#define ONE_FLT   (65536.0)
+
+void domono();
+void dostereo();
+void usage();
+void delinit(int cmndline,char *argv[]);
+void initfbuf();
+void initlbuf();
+
+int 	Funcflag = 0;
+int 	Quickflag = 1;
+long 	Nchans;
+long	Delaysamps;
+long    bufcnt;
+long	remainder;
+float 	*Fdelbuf;
+long	*Ldelbuf;
+double 	Sr;
+double	Delay;
+float	Feedforgain;
+float	Feedbckgain;
+float	Delaygain;
+float	Prescale = 1.0F;
+double	Invert =1.0;
+double  Rescale;
+long	Lfeedforgain;
+long	Lfeedbckgain;
+long	Ldelaygain;
+long	Lprescale;
+double  Trailtime = 0.;
+long	Trailsamps = 0;
+int	ifd;	 	/* input file descriptor */
+int	ofd;		/* output file descriptor */
+long 	isr;
+long	sampsize;
+short  	*inbuf;
+short	*outbuf;
+
+static char version[] =
+	"~CDPVER~GROUCHO: DELAY Portable version 3.0, REVISED BY R. DOBSON OCT 1993.";
+ 
+PROGRAM_NUMBER(0xf93c2705);
+
+void
+main(argc,argv)
+int argc;
+char *argv[];
+{
+int cmndline = 0;
+char *ptr;
+	
+	if(sflinit("delay")) {
+		sfperror("Delay: initialisation");
+		exit(1);
+	}
+    	while(argv[1] != 0 && argv[1][0] == '-') {
+	      switch(argv[1][1]) {
+	          case 'f':	
+		 	Quickflag--;
+		 	break;
+		  case 'p':
+			ptr = (char *) &(argv[1][2]);
+			if(*ptr==0)
+			{
+				usage();
+			}
+			Prescale = (float)atof(ptr);
+			if(Prescale==0.0)
+			{
+			    printf("\ndelay: -p error, prescale = 0.0!\n");
+			    usage();
+			}
+			break;
+		   case 'i':
+			Invert = -1.0;
+			break;
+	          default:
+	         	usage();
+	         	break;
+	      }
+	argv++; argc--;
+	}
+	if(argc < 3) usage();
+	if(argc>4)
+	{
+	    if(argc<7)
+		{
+		fprintf(stderr,"\ndelay: incorrect no of parameters");
+		usage();
+		}
+	    else
+		{
+		cmndline=1;
+		}
+	}
+	if((inbuf = 
+		(short *) malloc((BUFLEN) * sizeof(short))) == NULL) {
+		fprintf(stderr,"\ndelay: not enough memory for input buffer");
+		exit(1);
+	}
+	if((outbuf = 
+		(short *) malloc(BUFLEN * sizeof(short))) == NULL) {
+		fprintf(stderr,"\ndelay: not enough memory for output buffer");
+		exit(1);
+	}
+      	if((ifd = sndopen(argv[1])) < 0) {
+	  fprintf(stderr,"\ndelay: unable to open input file\n");
+ 	  exit(1);
+	} 
+	if(sndgetprop(ifd,"sample rate",(char *) &isr, sizeof(long)) < 0) {
+		fprintf(stderr,"\ndelay: failed to get sample rate");
+		exit(1);
+	}
+	if(sndgetprop(ifd,"channels",(char *) &Nchans, sizeof(long)) < 0) {
+		fprintf(stderr,"\ndelay: failed to get channel data");
+		exit(1);
+	}
+	if(sndgetprop(ifd,"sample type",(char *) &sampsize, sizeof(long)) < 0){
+		fprintf(stderr,"\ndelay: failed to get sample type");
+		exit(1);
+	}
+	if(Nchans > 2) {
+	  fprintf(stderr,"\ndelay: too many channels! Mono or stereo only.\n");
+	  exit(1);
+	}
+
+	Sr = (double) isr;
+
+	delinit(cmndline,argv);
+
+	if((ofd = sndcreat(argv[2],-1,sampsize)) < 0) {
+	  fprintf(stderr,"\ndelay: unable to open output file\n");
+	  exit(1);
+	}
+
+	if(!Quickflag) 
+		initfbuf();
+	else
+		initlbuf();
+	stopwatch(1);
+	if(Nchans==1) {domono();} else {dostereo();}
+	stopwatch(0);
+	sndclose(ifd);
+   	if(sndputprop(ofd,"sample rate",(char *) &isr, sizeof(long)) < 0){
+	fprintf(stderr,"\ndelay: failed to write sample rate");
+   	}
+   	if(sndputprop(ofd,"channels", (char *) &Nchans, sizeof(long)) < 0){
+	fprintf(stderr,"\ndelay: failed to write channel data");
+   	}
+	free(inbuf);
+	free(outbuf);
+	if(!Quickflag) {free(Fdelbuf);} else {free(Ldelbuf);}
+   	sndclose(ofd);
+	sffinish();
+}
+
+
+void
+domono()
+{
+	int 	i;	
+	long	rdsamps;
+	float	input,output;
+	long	linput,loutput;
+	long	ltemp1,ltemp2;
+	long	ipptr=0,opptr=0;
+
+	do		/* while rdsamps */
+	{
+	   	if((rdsamps = fgetsbuf(inbuf, BUFLEN, ifd)) < 0) {
+        	fprintf(stderr,"\ndelay: failure to read input file\n");
+          	exit(1);
+	    }
+	    
+	    for(i=0;i<rdsamps;i++)
+		{
+			
+				/* delay line */
+			input = inbuf[i] * Prescale;
+			output = (input * Feedforgain) +		  //dry signal
+				(Fdelbuf[opptr] * Delaygain);
+			Fdelbuf[ipptr] = Fdelbuf[opptr++] * Feedbckgain; 
+			Fdelbuf[ipptr++] += input ;
+			outbuf[i] = (short) output;				  //dry + wet
+
+			if(ipptr >= Delaysamps) ipptr -= Delaysamps;
+			if(opptr >= Delaysamps) opptr -= Delaysamps;
+			 
+			
+			if(ipptr < 0 || opptr < 0) {			
+				 printf(
+			  "\ninternal error, ipptr=%d,opptr=%d\n",ipptr,opptr);
+			}	
+		}	
+        if(fputsbuf(outbuf,rdsamps,ofd) < rdsamps) {
+		     fprintf(stderr,
+				"\ndelay: failure in writing to output file\n");
+				exit(1);
+		}
+		inform(rdsamps,Sr);
+	} while(rdsamps > 0);
+
+	    	/* now do trailer	*/
+    rdsamps=BUFLEN;
+	if (Trailsamps>0) {
+		do	{	/* while bufcnt */
+		    
+		  	if(!bufcnt) 
+				rdsamps = remainder;
+					    	
+			for(i=0;i<rdsamps;i++) {
+		       	output=(Fdelbuf[opptr] * Delaygain);
+		       	Fdelbuf[ipptr++]  = 
+				Fdelbuf[opptr++] *Feedbckgain;
+		       	outbuf[i] = (short) output;
+		       	if(ipptr >= Delaysamps) ipptr -= Delaysamps;
+		       	if(opptr >= Delaysamps) opptr -= Delaysamps;
+		    }
+					  		
+        	if(fputsbuf(outbuf,rdsamps,ofd) < rdsamps){
+				fprintf(stderr,
+					"\ndelay: failure in writing to output file\n");
+				exit(1);
+		  	}
+    		inform(BUFLEN,Sr);
+			bufcnt -=1;
+    	} while(bufcnt>=0);
+    }				/* end if Trailsamps */
+}				/* end domono	*/
+
+void
+dostereo()			/* do stereo delay */
+{
+int 	i;	
+long	rdsamps;
+float	Linput,Rinput;
+float 	Loutput,Routput;
+long	Llinput,Rlinput;
+long	Lloutput,Rloutput;
+long	Ltemp1,Ltemp2;
+long	Rtemp1,Rtemp2;
+long	Lipptr = 0,Ripptr = 1;
+long	Lopptr = 0,Ropptr = 1;
+
+    do			/* while rdsamps */
+    {
+	if((rdsamps = fgetsbuf(inbuf,BUFLEN,ifd)) < 0)
+	{
+		fprintf(stderr,"\ndelay: failure to read input file\n");
+		exit(1);
+	}
+	    for(i=0;i<(rdsamps-1);i+=2)
+	    {
+		if(!Quickflag)
+		{
+			Linput=inbuf[i]*Prescale;
+			Rinput=inbuf[i+1]*Prescale;
+			Loutput=(Linput*Feedforgain) + 
+				(Fdelbuf[Lopptr]*Delaygain);
+			Routput=(Rinput*Feedforgain) + 
+				(Fdelbuf[Ropptr]*Delaygain);
+			Fdelbuf[Lipptr] = Fdelbuf[Lopptr]*Feedbckgain;
+			Lopptr+=2;
+			Fdelbuf[Ripptr] = Fdelbuf[Ropptr]*Feedbckgain;
+			Ropptr+=2;
+			Fdelbuf[Lipptr]+=Linput;
+			Fdelbuf[Ripptr]+=Rinput;
+			Lipptr+=2;
+			Ripptr+=2;
+			outbuf[i] = (short) Loutput;
+			outbuf[i+1] = (short) Routput;
+
+			if(Lipptr >= Delaysamps) Lipptr -= Delaysamps;
+			if(Ripptr >= Delaysamps) Ripptr -= Delaysamps;
+			if(Lopptr >= Delaysamps) Lopptr -= Delaysamps;
+			if(Ropptr >= Delaysamps) Ropptr -= Delaysamps;
+			}
+		else
+		{
+			Llinput = inbuf[i] * Lprescale;
+			Rlinput = inbuf[i+1] * Lprescale;
+			Ltemp1 = Llinput >> 16;
+			Rtemp1 = Rlinput >> 16;
+			Ltemp2 = Ldelbuf[Lopptr] >> 16;
+			Rtemp2 = Ldelbuf[Ropptr] >> 16;
+				Lopptr += 2;
+				Ropptr += 2;
+			Lloutput = (Ltemp1 * Lfeedforgain) + 
+				(Ltemp2 * Ldelaygain);
+			Rloutput = (Rtemp1 * Lfeedforgain) +
+				(Rtemp2 * Ldelaygain);
+			Ldelbuf[Lipptr] = Ltemp2 * Lfeedbckgain;
+			Ldelbuf[Ripptr] = Rtemp2 * Lfeedbckgain;
+			Ldelbuf[Lipptr] += Llinput;
+			Ldelbuf[Ripptr] += Rlinput;
+				Lipptr += 2;
+				Ripptr += 2;
+			outbuf[i] = (short)(Lloutput>>16);
+			outbuf[i+1] = (short)(Rloutput>>16);
+
+			if(Lipptr >= Delaysamps) Lipptr -= Delaysamps;    
+			if(Ripptr >= Delaysamps) Ripptr -= Delaysamps;
+			if(Lopptr >= Delaysamps) Lopptr -= Delaysamps;
+			if(Ropptr >= Delaysamps) Ropptr -= Delaysamps;
+		}
+	}
+	if(fputsbuf(outbuf,rdsamps,ofd) < rdsamps) {
+		fprintf(stderr,"\ndelay: failure in writing to output file\n");
+		exit(1);
+	}
+	inform(rdsamps/2,Sr);
+   }while(rdsamps > 0);		/* end do */
+
+/* now do trailer */
+
+    rdsamps=BUFLEN;
+    if(Trailsamps>0)
+    {
+	do
+	{
+		if(!bufcnt) rdsamps = remainder;
+		if(!Quickflag)
+		    {
+		    for(i=0;i<(rdsamps-1);i+=2)
+			{
+			Loutput = Fdelbuf[Lopptr]*Delaygain;
+			Routput = Fdelbuf[Ropptr]*Delaygain;
+			Fdelbuf[Lipptr] = Fdelbuf[Lopptr]*Feedbckgain;
+			Fdelbuf[Ripptr] = Fdelbuf[Ropptr]*Feedbckgain;
+			Lopptr+=2;
+			Ropptr+=2;
+			Lipptr+=2;
+			Ripptr+=2;
+			outbuf[i] = (short) Loutput;
+			outbuf[i+1] = (short) Routput;
+			if(Lipptr >= Delaysamps) Lipptr -= Delaysamps;
+			if(Ripptr >= Delaysamps) Ripptr -= Delaysamps;
+			if(Lopptr >= Delaysamps) Lopptr -= Delaysamps;
+			if(Ropptr >= Delaysamps) Ropptr -= Delaysamps;
+			}
+		    }
+		else
+		    {
+		    for(i=0;i<(rdsamps-1);i+=2)
+		    	{
+			Ltemp2 = Ldelbuf[Lopptr] >> 16;
+				Lopptr+=2;
+			Rtemp2 = Ldelbuf[Ropptr] >> 16;
+				Ropptr+=2;
+			Lloutput = Ltemp2 * Ldelaygain;
+			Rloutput = Rtemp2 * Ldelaygain;
+			Ldelbuf[Lipptr] = Ltemp2*Lfeedbckgain;
+			Ldelbuf[Ripptr] = Rtemp2*Lfeedbckgain;
+			Lipptr +=2;
+			Ripptr +=2;			
+			outbuf[i] = (short)(Lloutput>>16);
+			outbuf[i+1] = (short)(Rloutput>>16);
+			if(Lipptr >= Delaysamps) Lipptr -= Delaysamps;
+			if(Ripptr >= Delaysamps) Ripptr -= Delaysamps;
+			if(Lopptr >= Delaysamps) Lopptr -= Delaysamps;
+			if(Ropptr >= Delaysamps) Ropptr -= Delaysamps;
+			}
+		    }
+		if(fputsbuf(outbuf,rdsamps,ofd) < rdsamps) {
+		fprintf(stderr,"\ndelay: failure in writing to output file\n");
+		exit(1);
+		}
+		inform(rdsamps/2,Sr);
+		bufcnt -=1;
+	    } while (bufcnt>=0);
+	}				/* end if */
+   }					/* end dostereo */		
+ 
+void
+usage()
+{
+	fprintf(stderr,"delay - A Groucho Program: $Revision: 3.3 $\n");
+	fprintf(stderr,
+"usage: delay [-f][-pN][-i] infile outfile [4 cmndline parameters]\n");
+	fprintf(stderr,
+"	where -f sets output in floatsams, -p sets prescale = N for infile,\n");
+	fprintf(stderr,
+"	and -i inverts dry signal (for allpass filters).\n");
+	fprintf(stderr,"\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n",
+"Cmndline option: after [flags] filenames, enter all parameters in this order:",
+"	delaytime mix feedbackgain trailertime",
+"	 0.0  	< delaytime (msecs) 	<= maxdelay",
+"	 0.0 	<= mix		 	<= 1.0",
+"	-1.0 	<= feedbackgain    	<= 1.0",
+"	 0.0 	<= trailertime (secs)	<= 30.0",
+"where maxdelay =  8000ms @22050 mono",
+"               =  4000ms @22050 stereo, and 44100 mono",
+"               =  2000ms @44100 stereo"); 
+	exit(0);
+}
+
+void
+delinit(int cmndlin, char *argv[])
+{
+	double temp,maxdelay,mix;
+	char msg[80];
+	/*	Some jiggery pokery is needed in setting the
+	 *	fixed point variables to avoid overloading
+	 *	the input to the delay line. The input is
+	 *	divided by two in conjunction with prescaling
+	 *	and the feedthrough and delay outputs are
+	 *	raised by two to compensate.
+	 */
+	maxdelay = ((double) (DELBUFSIZ) / Sr) * 1000.0; /* milliseconds */
+	if(Nchans==2) {maxdelay /= 2;}
+	if(cmndlin) 
+	{
+		temp = atof(argv[3]);
+		if (temp > maxdelay)
+		{
+			printf("delay: delay time too long\n");
+			exit(1);
+		}
+		Delaysamps = (long) ((Sr *temp) / 1000.0);
+		mix = atof(argv[4]);
+		if((mix < 0.0) || (mix > 1.0))
+		{
+			printf("delay: mix value out of range\n");
+			exit(1);
+		}
+		Feedbckgain = (float)atof(argv[5]);
+		if((Feedbckgain < -1.0) || (Feedbckgain > 1.0))
+		{
+			printf("delay: feedbackgain out of range\n");
+			exit(1);
+		}
+		Trailtime = atof(argv[6]);
+		if((Trailtime < 0.0) || (Trailtime > 1000.0))
+		{
+			printf("delay: trailtime out of range\n");
+			exit(1);
+		}
+				
+	}
+	else		/* ask for params from terminal */
+	{
+		printf("max delay in msec = %f\n",maxdelay);
+		sprintf(msg,
+"Give delay time in milliseconds			:");
+		temp  = accexpr(msg,0.0,maxdelay,0.0,0,Sr);
+		if(temp==0.0)
+		{
+	printf("\nzero delay time!!! No point in running program. Bye!\n");
+	exit(1);
+		}
+		Delaysamps = (long) ((Sr * temp) / 1000.0);
+		printf("\ndelay length in samples = %ld \n",Delaysamps);
+		sprintf(msg,
+"Give mix proportion (0.0 <= level <= 1.0)		:");
+		mix = accexpr(msg,0.0,1.0,0.0,0,Sr);
+		sprintf(msg,
+"Give feedback gain (-1.0 <= level <= 1.0)		:");
+		Feedbckgain = (float)accexpr(msg,-1.0,1.0,0.0,0,Sr);
+		sprintf(msg,
+"Give trailer time (default = 0 <= trailtime <= 30.0)	:");
+		Trailtime = accexpr(msg, 0.0,30.0,0.0,1,Sr);
+	}
+/* now massage parameters 			*/
+/* nb delaysamps is per channel here 	*/
+ 
+	Feedforgain = (float)((1.0 - mix)*Invert);	
+	Lfeedbckgain = (long) (Feedbckgain * ONE_FLT);
+	Delaygain = (float)mix;
+	Ldelaygain = (long) (Delaygain * ONE_FLT * 2.0);	
+	Rescale = (1.0 / (Feedbckgain + 1.0)); /* i/p compensation */
+	Prescale *= (float)Rescale;
+	Lprescale= (long) (Prescale * ONE_FLT * .5);
+	Feedforgain /= (float)Rescale;
+	Lfeedforgain = (long) (Feedforgain * ONE_FLT * 2.0);	
+	Trailsamps = (long) (Trailtime * Sr);
+	if(Nchans==2) Trailsamps *= 2;
+	bufcnt = Trailsamps/BUFLEN;
+	remainder = Trailsamps % BUFLEN;
+}
+
+void
+initfbuf()
+{
+	Delaysamps*=Nchans;
+	if((Fdelbuf = (float *) calloc(Delaysamps+2, sizeof(float))) == NULL){
+	   fprintf(stderr,"\ndelay: not enough memory for main delay buffer");
+	   exit(1);
+	}
+}
+
+void
+initlbuf()
+{
+	Delaysamps*=Nchans;
+	if((Ldelbuf = (long *) calloc(Delaysamps+4, sizeof(long))) == NULL){
+	    fprintf(stderr,"\ndelay: not enough memory for main delay buffer");
+	    exit(1);
+	}
+}

dev/externals/reverb/diffuse.h

@@ -0,0 +1,62 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+//diffuse.h
+//definition of classes for small, meduim large dense reverberators
+// based on gardner, using nested allpasses
+#ifndef __DIFFUSE_INCLUDED__
+#define __DIFFUSE_INCLUDED__
+
+#include "reverberator.h"
+#include "tapdelay.h"
+
+typedef struct nested_allpass_data{
+	unsigned int time1;
+	unsigned int time2;
+	unsigned int time3;
+	double gain1;
+	double gain2;
+	double gain3;
+}
+NESTDATA;
+
+
+class small_diffuser {
+public:
+	small_diffuser(unsigned int pre_delay,
+					const NESTDATA *apdata1,
+					const NESTDATA *apdata2,double lp_gain,double lp_coeff);
+	virtual ~small_diffuser();
+	bool create(void);
+	float tick(float insam);
+	//void clear(void);
+
+private:
+	NESTDATA ap1_data,ap2_data;
+	nested_allpass *ap1,*ap2;
+	delay *predelay;
+	lowpass1 *lp1;
+	unsigned int predelay_time;
+	double lpgain,lpcoeff;
+	float out1,out2;
+};
+
+
+#endif

dev/externals/reverb/lpcomb.cpp

@@ -0,0 +1,169 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+//lpcomb.cpp
+extern "C" {
+#include <sfsys.h>
+#include <cdplib.h>	   //NB requires stdio.h etc - time to change this?
+}
+
+#include "reverberator.h"
+
+//straight from Csound: nreverb (vdelay.c)
+static bool prime( int val )
+{
+    int i, last;
+	
+    last = (int)sqrt( (double)val );
+    for ( i = 3; i <= last; i+=2 ) {
+      if ( (val % i) == 0 ) return false;
+    }
+    return true;
+}
+
+static int to_prime(double dur,double sr)
+{
+		int time = (int) (dur * sr);
+		if(time % 2== 0) time += 1;
+		while(!prime(time)) 
+			time +=2;
+		return time;
+}
+
+void
+usage()
+{
+
+	printf("\nusage: lpcomb  infile outfile looptime rvbtime filtgain");
+	exit(0);
+}
+
+void
+main(int argc,char *argv[])
+
+{
+	int		ifd, ofd;	
+	//double  delaytime,gain;
+	double	rvbtime,gain2,filtgain,dloop; 	
+	long	isr,chans,sampsize,ilooptime;
+	double	sr;
+	lpcomb	*comb1 = 0;
+
+
+	if(sflinit("lpcomb")){
+		sfperror("lpcomb: initialisation");
+		exit(1);
+	}
+
+
+	if(argc < 6)
+		usage();
+	dloop = atof(argv[3]);
+	if(dloop <= 0.0){
+		fprintf(stderr,"\nlpcomb: invalid looptime parameter");
+		exit(1);
+	}
+
+	rvbtime = atof(argv[4]);
+	if(rvbtime <= 0.0){
+		fprintf(stderr,"\nlpcomb: invalid delay parameter");
+		exit(1);
+	}
+	filtgain = atof(argv[5]);
+	if((filtgain < 0.0) || (filtgain > 1.0)){
+		fprintf(stderr,"\nlpcomb: filter gain %lf out of range",filtgain);
+		exit(1);
+	}
+
+	if(argv[1] != 0) {
+	    if((ifd = sndopen(argv[1])) < 0) {
+		   fprintf(stderr,"Allpass: cannot open input file %s\n", argv[1]);
+		   exit(1);
+		}
+	}
+	if(sndgetprop(ifd,"sample rate",(char *) &isr, sizeof(long)) < 0){
+		fprintf(stderr,
+		  "\nAllpass: failed to read sample rate for file %s",argv[1]);
+		exit(1);
+	}
+	sr = (double) (isr);
+	if(sndgetprop(ifd,"channels",(char *) &chans, sizeof(long)) < 0){
+		fprintf(stderr,
+		  "\nAllpass: failed to read channel data for file %s", argv[1]);
+		exit(1);
+	}
+	if(chans != 1){
+		fprintf(stderr,"\nAllpass works only on mono files");
+		exit(1);
+	}
+	if(sndgetprop(ifd,"sample type",(char *) &sampsize, sizeof(long)) < 0){
+		fprintf(stderr,
+		  "\nAllpass: failed to read sample type for file %s", argv[1]);
+		exit(1);
+	}
+
+
+	if((ofd = sndcreat(argv[2],-1,sampsize)) < 0) {
+		  fprintf(stderr,"Cannot open output file %s\n",argv[2]);
+		  exit(1);
+	}
+	ilooptime  = to_prime(dloop, sr);	
+	gain2 = exp( log(0.001)* (dloop/rvbtime)) * (1. - filtgain) ;
+	printf("\ngain2 = %lf\tfiltgain = %lf",gain2,filtgain);
+	comb1 = new lpcomb(ilooptime,gain2,filtgain,1.0);
+	if(!comb1->create()){
+		printf("\nFailed to create comb filter");
+		exit(1);
+	}
+
+	stopwatch(1);
+
+	for(;;){
+		int rv;
+		float ip,op;
+
+		if((rv = fgetfloat(&ip,ifd)) < 0){
+			fprintf(stderr,"\nError in reading file"); 
+			exit(1);
+		}
+		if(!rv)
+			break;
+		op = comb1->tick(ip);
+		if(fputfloat(&op,ofd) < 1){
+			fprintf(stderr,"\nError writing to output file"); 
+			exit(1);
+		}
+	}
+	stopwatch(0);
+	sndclose(ifd);
+	if(sndputprop(ofd,"sample rate", (char *) &isr, sizeof(long)) < 0){
+		fprintf(stderr,"\nAllpass: failed to write sample rate");
+	}
+	if(sndputprop(ofd,"channels",(char *) &chans, sizeof(long)) < 0){
+		fprintf(stderr,"\nAllpass: failed to write channel data");
+	}
+	sndclose(ofd);
+	delete comb1;
+	sffinish();
+}

dev/externals/reverb/readme.txt

@@ -0,0 +1,121 @@
+readme file for CDP reverb programs
+
+ROOMRESP
+
+simple utility program to generate early reflections for reverbs, using "mirror room" model.
+
+The format for early reflection data is a series of lines containing
+<time> <value> pairs. The first time is non-zero (since the direct sound
+is presumed to be there):
+
+0.0082 0.500000
+0.0121 0.204198
+0.0134 0.169094
+0.0145 0.143638
+0.0170 0.093391
+0.0180 0.084478
+0.0200 0.067530
+.. etc
+
+Data obtained from roomresp is used internally for the early reflections in both "rmverb" and "reverb". 
+Roomresp usually generates many more taps than are necessary or efficient; typically the first
+19 taps (say) will be used. The two reverb programs offer the facility to import reflection data from
+roomresp directly into the reverb (see the -e flag option).  As a FIR filter is used to process
+these, the CPU load of the program is directly affected by the number of reflections used.
+
+the choice of early reflections has an influence on the character and tonality of the reverb 
+- experimentation is recommended here. Hand-written reflection data is especially good for "creative" reverbs
+with heavy ringing and colour!
+
+RMVERB
+
+The program "rmverb" uses nested allpass filters based on Bill Gardner's published models; hence three 
+room models are offered, "small", "medium" and "large".
+ 
+The "large room" model is modified with smaller delays (so really "large-ish"), to reduce the 
+obviousness of the slap-back echoes generated by the model. 
+
+Note that this relies on the fact that each reverb element is implemented in a self-contained manner, 
+unit-generator style (e.g. including plain delays between elements). If it is required to optimise 
+the architecture by having all elements share one delay line,  this change will not work - the outer 
+delay will need to be smaller than the  inner ones.
+
+Two allpass filters per channel are aded to the output, to obtain
+a stereo effect from the internal mono reverb. According to Gardner, his concept was to use one "diffuser" 
+(configuration of nested allpass filters) per speaker; this suggests that any thinness in the output can 
+be overcome (at the cost of processing time) by paralleling several diffusers per channel, with slightly 
+different parameters.
+
+Additionally (NEW from CDP!) each nested allpass filter has been extended with the addition of a first-order
+filter in the outer feedback loop. This acts in addition to the global feedback loop defined by Gardner.
+This is indicated in "rmverb" by the -d flag ("double damping"). This seems to reduce the metallic 
+ringing at the very end of the reverb decay, at the cost of shortening the overall decay time. It would
+easily be possible, with this change, to eliminate the global feedback filter and rely entirely on the
+filter inside the nested allpasses, to recover some of the decay time. This feature is experimental.
+
+As Gardner notes, pushing the feedback amount to obtain long decay times is not recommended, as the timbre
+deteriorates. The program is implemented in a relatively "raw" state; with a raw feedback level parameter 
+rather than a "reverb time"; some mapping will need to be computed for this!
+
+The programs also allow the delay between direct sound and first reflection (as defined in the reflection data)
+to be modified - see the -p flag option.	  Standard filters are provided to filter the input signal itself,
+and the signal entering the reverb. The latter is usual, the former is rarely necessary. 
+Some filters are second-order, offering a 12dB slope (currently fixed in the program). It would be easy to 
+make the slope variable as another user parameter.	 The sources include a library of low and high pass filters,
+first order and second order (shelving).
+
+
+ REVERB
+
+The program "reverb" uses a standard moorer/schroeder comb+allpass structure, with	6 comb filters 
+(with lp filters in the fedback loop) feeding 4 allpass filters. A conventional,"reverb time" is provided, but
+no "room size". This wouild have to be obtained via more or less extensive early reflections. 
+
+
+Not shown in the programs usage message is a possble extra argument which gives the name of a text data file
+containing delay times (in msecs) for the comb and allpass filters.	 This will be a plain text file containing
+10 positive floating-point values. The internal preset times corresponds to this data file:
+
+50 56 61 68 72 78 20 14 9 6
+
+ - which gives the canonical Moorer values.
+
+
+NOTE one importnat difference in the "absorb" argument to the programs.
+Reverb:  range 0 - 1, large values = more absorption.
+Rmverb:  value in Hzf for feedback filter (See Gardner for details). He suggests:
+	small room: 4200 Hz
+	medium:  2500 Hz
+	large:  2600 Hz
+
+Rmverb offers direct setting of this value. This then applies both the the global feedback loop
+and to the internal nested-allpass loop as described above.
+
+In both reverb programs, ~all~ delay times are adjusted to give prime number integral lengths in samples. This means
+that performance will not be ~identical~ at all sample rates. Naturally, if the sample rate is to be fixed, in a
+dsp implementation, the possibility exists to fine-tune delays for that sample rate.
+
+
+Both programs generate multi-channel (surround) signals by using a separate allpass filter on each channel for
+de-correlation. The direct sound is fed only to the front stereo pair.
+
+
+Note that these are command-line programs, and the various dsp elements are not fully re-entrant in terms of 
+running re-initializations.	 This would not be hard to do, howewver, especially given the goal is a 
+re-implementation on a DSP architecture.
+
+CDP plans to extend these algorithms with variable delay lines, and also to increase density by addition of 
+elements.  To say nothing of real-time versions in plugin form!
+
+ Richard Dobson
+ CDP
+ September 14th 2004
+
+
+
+
+
+
+
+
+

dev/externals/reverb/reflect.cpp

@@ -0,0 +1,189 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+//NB it is most important to avoid periodicities in these
+// using rmresp, best to place source and listener near opposite corners, or both near one corner 
+//preset early reflections - could get these from a brkfile, of course
+//except small room ; seems best to have src and listener same end of room, to get even earlier reflections
+//11 8 3.3
+//0.9
+//5
+//9 3 1.75
+//9 5.7 1.29
+//normed to 0.5
+
+deltap smalltaps[] = {
+	{0.0082,0.500000},
+	{0.0121,0.204198},
+	{0.0134,0.169094},
+	{0.0145,0.143638},
+	{0.0170,0.093391},
+	{0.0180,0.084478},
+	{0.0200,0.067530},
+	{0.0219,0.063096},
+	{0.0226,0.053176},
+	{0.0233,0.044833},
+	{0.0237,0.048586},
+	{0.0243,0.046058},
+	{0.0250,0.043713},
+	{0.0256,0.037388},
+	{0.0261,0.044475},
+	{0.0265,0.034818},
+	{0.0271,0.033360},
+	{0.0276,0.035772},
+	{0.0281,0.034382},
+	{0.0285,0.030051},
+	{0.0287,0.033059},
+	{0.0301,0.027091},
+	{0.0304,0.026513},
+	{0.0306,0.026200},
+	{0.0309,0.023001},
+	{0.0319,0.024114},
+	{0.0326,0.020657},
+	{0.0335,0.021782},
+	{0.0340,0.019019},
+	{0.0356,0.017387},
+	{0.0398,0.017155},
+	{0.0408,0.014690},
+	{0.0412,0.014424},
+	{0.0416,0.014160},
+	{0.0425,0.012175},
+	{0.0429,0.011971},
+	{0.0439,0.011468},
+	{0.0451,0.010854}
+	};
+
+//dim: 21 16 5
+// refl: 0.95
+//nrooms 5
+//listenerpos
+//17 8 2
+//srcpos:
+//6 9.35 2
+//RWD NB does sound worse with with only 32 taps, or refl-order 3
+// this is the 'no-compromise' reverb!
+deltap mediumtaps[] = 
+{
+	{0.0332,0.993303},
+	{0.0353,0.834879},
+	{0.0377,0.729745},
+	{0.0447,0.494138},
+	{0.0548,0.345328},
+	{0.0561,0.313134},
+	{0.0570,0.319440},
+	{0.0577,0.296282},
+	{0.0583,0.290651},
+	{0.0598,0.276075},
+	{0.0615,0.274630},
+	{0.0625,0.240116},
+	{0.0627,0.251368},
+	{0.0641,0.240392},
+	{0.0644,0.226004},
+	{0.0684,0.200374},
+	{0.0690,0.218340},
+	{0.0700,0.201354},
+	{0.0713,0.194250},
+	{0.0718,0.191280},
+	{0.0728,0.176802},
+	{0.0740,0.171020},
+	{0.0752,0.165815},
+	{0.0770,0.166316},
+	{0.0778,0.147079},
+	{0.0780,0.154272},
+	{0.0791,0.149852},
+	{0.0816,0.148066},
+	{0.0825,0.137700},
+	{0.0826,0.130403},
+	{0.0836,0.134167},
+	{0.0862,0.132653} ,
+	{0.0870,0.117789},
+	{0.0871,0.123637},
+	{0.0881,0.120782},
+	{0.0913,0.106847},
+	{0.0929,0.114356},
+	{0.0937,0.106863},
+	{0.0946,0.104723},
+    {0.0975,0.103831}
+	
+};
+
+//normed to 0.5:
+//43 31 19
+//0.95
+//7
+//7.6 15 2
+//30 23.7 6
+
+deltap largetaps[] = 
+{
+	{0.0729,0.500000},
+	{0.0758,0.439456},
+	{0.0975,0.265770},
+	{0.0997,0.241550},
+	{0.1151,0.190794},
+	{0.1162,0.187246},
+	{0.1180,0.172387},
+	{0.1247,0.154486},
+	{0.1320,0.137708},
+	{0.1330,0.135754},
+	{0.1344,0.139864},
+	{0.1346,0.125902},
+	{0.1360,0.129780},
+	{0.1404,0.115604},
+	{0.1449,0.114371},
+	{0.1464,0.112074},
+	{0.1477,0.115801},
+	{0.1492,0.107884},
+	{0.1540,0.096167}
+};
+
+//Moorer
+deltap longtaps[18] = {
+					{0.0043,0.841},
+					{0.0215,0.504},
+					{0.0225,0.491},
+					{0.0268,0.379},
+					{0.0270,0.380},
+					{0.0298,0.346},
+					{0.0458,0.289},
+					{0.0485,0.272},
+					{0.0572,0.192},
+					{0.0587,0.193},
+					{0.0595,0.217},
+					{0.0612,0.181},
+					{0.0707,0.180},
+					{0.0708,0.181},
+					{0.0726,0.176},
+					{0.0741,0.142},
+					{0.0753,0.167},
+					{0.0797,0.134}
+};
+
+
+deltap shorttaps[] = 
+{
+					{0.0099,1.02},			  //my times!
+					{0.0154,0.818},
+					{0.0179,0.635},
+					{0.0214,0.719},
+					{0.0309,0.267},
+					{0.0596,0.242}
+					};
+

dev/externals/reverb/reflect.h

@@ -0,0 +1,33 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+//reflect.h
+#ifndef __REFLECT_INCLUDED
+#define __REFLECT_INCLUDED
+#include "reverberator.h"
+
+
+extern deltap smalltaps[];
+extern deltap mediumtaps[];
+extern deltap largetaps[];
+extern deltap longtaps[];
+extern deltap shorttaps[];
+#endif
+

dev/externals/reverb/reverb.cpp

@@ -0,0 +1,745 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+/* MODULE FOR REVERB.EXE */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+#include <algorithm>
+#include <vector>
+using namespace std;
+/***** FORMULA FOR LOOP GAIN *****
+
+  g = pow(10, - ( ((60*ilooptime)/reverbtime) / 20))
+
+  derived from Moore:
+
+	reverbtime = 60/(-20.log10(g)) * ilooptime
+
+  then we set igain = g(1 - filtgain), officially...
+
+  */
+extern "C" {
+#include <sfsys.h>
+#include <osbind.h>
+#include <cdplib.h>	   //NB requires stdio.h etc - time to change this?
+}
+#include "reverberator.h"
+#include <string.h>
+void usage(void);
+long readtaps(FILE *fp, deltap **taps,double sr);
+
+//my library of early reflections...
+//by doing it this way I can use the sizeof() operator on the array names
+#include "reflect.cpp"
+
+enum cmdargs {
+	PROGNAME,
+	INFILE,
+	OUTFILE,
+	RGAIN,
+	MIX,
+	REVERBTIME,	
+	DAMPING,
+	LPFREQ,
+	TRTIME,
+	CONFIG
+};
+
+const char* cdp_version = "5.0.0";
+
+///////////// Moorer reverb
+int
+main(int argc,char *argv[])
+
+{
+	int		ifd=-1, ofd=-1,i;	
+	//double  delaytime,gain;
+	double	rvbtime, damping, trailertime = 0.0; 	
+	long	chans,outchans = 2;
+	double	sr;
+	//all this stuff from Brandt  Csound implementation of posh Moorer
+	MOORERDATA mrdata = {0.050,0.056,0.061,0.068,0.072,0.078,0.02,0.014,0.009,0.006};	
+	double predelay = 0.0;
+	moorer *moorerverb = 0;
+	tapdelay *tdelay = 0;
+	deltap *ptaps = 0;	
+	allpassfilter **chan_ap = 0;		//ptr to array of allpassfilter ptrs
+	tonecontrol *tc_l_lowcut = 0, *tc_r_lowcut = 0;
+	// could be replaced with a tonecontrol for deeper slope
+	onepole *l_highcut = 0, *r_highcut = 0, *lp = 0;
+	FILE *fp_earlies = 0; 
+	long ntaps;
+	bool want_mono = false;
+	bool want_floats = false;
+	bool usertaps = false;
+	double lowcut_freq = 0.0,highcut_freq = 0.0,lp_freq = 0.0;
+	float dry_gain,wet_gain,diff_gain;
+	CHPEAK *peaks;
+	SFPROPS props;
+	// test vcomb with variable delay taps 
+	//vmtcomb4 vcomb;
+	//fastlfo lfo1,lfo2,lfo3,lfo4;
+	if((argc==2) && strcmp(argv[1],"--version")==0) {
+		fprintf(stdout,"%s\n",cdp_version);
+		fflush(stdout);
+		return 0;
+	}
+
+	if(sflinit("reverb")){
+		sfperror("reverb: initialisation");
+		exit(1);
+	}
+	if(argc < CONFIG) usage();
+
+	while((argc > 1) && argv[1][0]=='-'){
+		char *arg = argv[1];
+		switch(*(++arg)){
+		case('p'):
+			if(*(++arg) =='\0'){
+				fprintf(stderr,"\npredelay requires a parameter");
+				usage();
+			}			
+			predelay = atof(arg);
+			if(predelay <0.0){
+				fprintf(stderr,"\npredelay must be >= 0.0");
+				usage();
+			}
+			predelay *= 0.001;			//arg is in msecs
+			break;
+		case('c'):			
+			if(*(++arg) == '\0') {
+				fprintf(stderr,"\n:reverb: -c flag requires a value");
+				usage();				
+			}
+			outchans = atoi(arg);
+			if(outchans < 1 || outchans > 16){
+				fprintf(stderr,"\nreverb: impossible channel value requested");
+				usage();
+			}
+			if(outchans==1)
+				want_mono = true;
+			break;
+		case('e'):
+			if(*(++arg)=='\0'){
+				fprintf(stderr,"\nreverb: -e flag needs a filename");
+				usage();
+			}
+			if((fp_earlies = fopen(arg,"r")) ==0){
+				fprintf(stderr,"\nreverb: unable to open breakpoint file %s",arg);
+				exit(1);
+			}
+			usertaps = true;
+			break;
+		case('f'):
+			want_floats = true;
+			if(*(++arg) != '\0')
+				fprintf(stderr,"\nreverb: WARNING: -f flag does not take a parameter");
+			break;
+		case('L'):
+			if(*(++arg) == '\0'){
+				fprintf(stderr,"\nreverb: Lowcut flag needs frequency argument");
+				usage();
+			}			
+			lowcut_freq = atof(arg);
+			if(lowcut_freq <= 0.0){
+				fprintf(stderr,"\nreverb: Lowcut freq must be greater than zero");
+				usage();
+			}
+			break;
+		case('H'):
+			if(*(++arg) == '\0'){
+				fprintf(stderr,"\nreverb: Highcut flag needs frequency argument");
+				usage();
+			}			
+			highcut_freq = atof(arg);
+			if(highcut_freq <= 0.0){
+				fprintf(stderr,"\nreverb: Highcut freq must be greater than zero");
+				usage();
+			}
+			break;
+		default:
+			fprintf(stderr,"\nreverb: illegal flag option %s",argv[1]);
+			usage();
+			break;
+		}
+		argc--;
+		argv++;
+	}
+	peaks = (CHPEAK *) malloc(sizeof(CHPEAK) * outchans);
+	if(peaks==NULL){
+		fputs("\nNo memory for PEAK data\n",stderr);
+		exit(1);
+	}
+	chan_ap = new allpassfilter*[outchans];
+	if(chan_ap==0){
+		fputs("\nreverb: no memory for multi-channel diffusion",stderr);
+		exit(1);
+	}
+
+
+	if(argc < CONFIG)
+		usage();
+
+	if(argv[INFILE] != 0) {
+	      	if((ifd = sndopenEx(argv[INFILE],0,CDP_OPEN_RDONLY)) < 0) {
+		   fprintf(stderr,"\nreverb: cannot open input file %s\n", argv[INFILE]);
+		   exit(1);
+		}
+	}
+	if(!snd_headread(ifd,&props)){
+		fprintf(stderr,"\nerror reading infile header: %s\n",rsferrstr);
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	if(props.type != wt_wave){
+		fprintf(stderr,"\ninfile is not a waveform file");
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	
+	sr = (double) props.srate;
+	chans = props.chans;
+
+	if(chans > 2){
+		fprintf(stderr,"\nreverb works only on mono or stereo files");
+		exit(1);
+	}
+	
+	
+	diff_gain = atof(argv[RGAIN]);
+	if(diff_gain < 0.0 || diff_gain > 1.0 ){
+		printf("\nreverb: rgain must be >= 0.0 and  <= 1.0");
+		usage();
+	}
+	dry_gain = (float) atof(argv[MIX]);			  //global output gain from diffuser
+	if(dry_gain < 0.0 || dry_gain >= 1.0 ){
+		printf("\nreverb: mix must be  between 0.0 and 1.0");
+		usage();
+	}		
+	wet_gain = 1.0f - dry_gain;
+	//probably not very scientific, but it works intuitively...
+	//wet_gain *= 2;		  
+	
+
+	rvbtime = atof(argv[REVERBTIME]);	  
+	if(rvbtime <= 0.0){
+	   fprintf(stderr,
+		"\nreverb: rvbtime must be > 0\n");
+		exit(1);
+	}
+
+	damping = atof(argv[DAMPING]);
+	if(damping < 0.0 || damping > 1.0){
+		fprintf(stderr,"\nreverb: absorb must be in the range 0.0 - 1.0\n");
+		exit(1);
+	}
+	double filt_limit = sr / 3.0;
+	lp_freq = atof(argv[LPFREQ]);
+	if(lp_freq < 0.0 || lp_freq > filt_limit){
+		printf("\napass: lp_freq must be within 0.0 to %.4lfHz",filt_limit);
+		usage();
+	}
+
+	trailertime = atof(argv[TRTIME]);
+	if(trailertime < 0.0)
+		trailertime = 0.0;
+	
+	if(argc==CONFIG+1){
+		int got = 0;
+		double dtime1,dtime2,dtime3,dtime4,dtime5,dtime6,
+			adtime1,adtime2,adtime3,adtime4;
+		FILE *fp = fopen(argv[CONFIG],"r");
+		if(!fp)
+			printf("\nreverb: can't open datafile %s: using presets",argv[CONFIG]);
+		else{
+			got = fscanf(fp,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
+				&dtime1,
+				&dtime2,
+				&dtime3,
+				&dtime4,
+				&dtime5,
+				&dtime6,
+				&adtime1,			 //allpasses
+				&adtime2,
+				&adtime3,
+				&adtime4				
+				);
+			fclose(fp);
+			if(got != 10)
+				printf("\nreverb: error reading looptime values");
+			else {
+				//create prime delay times from the msec vals
+				//get this in a loop later on...				
+				//RWD TODO: check times are not too huge!
+				mrdata.ctime1 = dtime1 * 0.001;
+				mrdata.ctime2 = dtime2 * 0.001;
+				mrdata.ctime3 = dtime3 * 0.001;
+				mrdata.ctime4 = dtime4 * 0.001;
+				mrdata.ctime5 = dtime5 * 0.001;
+				mrdata.ctime6 = dtime6 * 0.001;
+				mrdata.atime1 = adtime1 * 0.001;
+				mrdata.atime2 = adtime2 * 0.001;
+				mrdata.atime3 = adtime3 * 0.001;
+				mrdata.atime4 = adtime4 * 0.001;
+			}
+		}
+	}
+
+	//create input low/high filters, if wanted
+	if(lowcut_freq > 0.0){
+		if(highcut_freq > 0.0 && highcut_freq <= lowcut_freq){
+			fprintf(stderr,"\nreverb: Highcut frequency	 must be higher than Lowcut frequency");
+			usage();
+		}
+		//lowcut is based on low shelving eq filter; here fixed at 12dB
+		tc_l_lowcut = new tonecontrol(lowcut_freq,-12.0,LOW_SHELVE,sr);
+		if(!tc_l_lowcut->create()){
+			fprintf(stderr,"\nreverb: unable to create Lowcut filter");
+			exit(1);
+		}
+		if(chans==2){
+			tc_r_lowcut = new tonecontrol(lowcut_freq,-12.0,LOW_SHELVE,sr);
+			if(!tc_r_lowcut->create()){
+				fprintf(stderr,"\nreverb: unable to create Lowcut filter");
+				exit(1);
+			}
+		}
+	}
+
+	if(highcut_freq > 0.0){		
+		l_highcut = new onepole(highcut_freq,sr,LOW_PASS);
+		if(chans==2)
+			r_highcut = new onepole(highcut_freq,sr,LOW_PASS);
+	}
+
+	if(lp_freq > 0.0)
+		lp		= new onepole(lp_freq,sr,LOW_PASS);
+
+
+	moorerverb =  new moorer(&mrdata,rvbtime,damping,sr);
+	//further external allpasses will be used to create multi-channel diffusion
+	
+	for(i = 0; i < outchans; i++){
+		chan_ap[i] = new allpassfilter(sr,to_prime(0.005 + 0.00025 * (double)i,sr),0.4 ,0);
+		if(chan_ap[i] == 0) {
+			fprintf(stderr,"\nreverb: no memory for output diffusers");
+			exit(1);
+		}
+		if(!chan_ap[i]->create()){
+			fprintf(stderr,"\nreverb: no memory for output diffusers");
+			exit(1);
+		}
+	}
+	if(usertaps){
+#ifdef _DEBUG
+		assert(fp_earlies);
+#endif
+		ntaps = readtaps(fp_earlies,&ptaps,sr);
+		if(ntaps==0){
+			fprintf(stderr,"\nreverb: error reading tapfile");
+			exit(1);
+		}
+		printf("\nloaded %ld early reflections from tapfile",ntaps);
+		fclose(fp_earlies);
+		fp_earlies = 0;
+	}
+	else {
+		if(rvbtime <= 0.6){
+			ntaps = sizeof(smalltaps) / sizeof(deltap);
+			ptaps = smalltaps;
+			printf("\nsmall-room: using %ld early reflections",ntaps);
+		}
+		else if(rvbtime > 0.6 && rvbtime <= 1.3){
+			ntaps = sizeof(mediumtaps) / sizeof(deltap);
+			ptaps = mediumtaps;
+			printf("\nmedium-room: using %ld early reflections",ntaps);
+		}
+		else{
+			ntaps = sizeof(largetaps) / sizeof(deltap);
+			ptaps = largetaps;
+			printf("\nlarge-room: using %ld early reflections",ntaps);
+		}
+	}
+	tdelay = new tapdelay(ptaps,ntaps,sr);
+
+	//if user prescribes a predelay, adjust all taptimes
+	if(predelay > 0.0){
+		double current_delay = ptaps[0].pos;
+		double adjust = predelay - current_delay;
+		for(int i = 0; i < ntaps; i++)
+			ptaps[i].pos += adjust;
+	}
+	
+	if(!moorerverb->create()){
+		fprintf(stderr,"\nreverb: unable to create reverberator");
+		exit(1);
+	}
+	if(!tdelay->create()){
+		fprintf(stderr,"\nreverb: unable to create early reflections");
+		exit(1);
+	}
+	
+	//double maxgain = tdelay->get_maxgain();
+	//if(maxgain > 1.0)
+	//	tdelay->set_gain(maxgain * early_gain);
+	printf("\npredelay = %.4lf secs",ptaps[0].pos);
+	//wet_gain *= floor(maxgain);
+	if(want_floats)
+		props.samptype =FLOAT32;
+	props.chans = outchans;
+
+	if((ofd = sndcreat_ex(argv[OUTFILE],-1,&props,SFILE_CDP,CDP_CREATE_NORMAL)) < 0) {
+		  fprintf(stderr,"\nreverb: cannot open output file %s\n",argv[OUTFILE]);
+		  exit(1);
+	}
+
+
+	printf("\n");
+	long outsamps = 0;
+	long step = (long) (0.25 * sr);
+#ifdef NOTDEF
+	// vcomb test
+	long seed = lfo1.init((double)inprops.srate,0.0,seed);
+	lfo1.set_WaveType(LFO_RANDOM);
+	lfo1.set_freq(0.5);
+	lfo1.set_mod(1.0); /* no change - we control modf range externally */
+	lfo2.init((double)inprops.srate,0.0,seed);
+	lfo2.set_WaveType(LFO_RANDOM);
+	lfo2.set_freq(0.55);
+	lfo2.set_mod(1.0);
+	lfo3.init((double)inprops.srate,0.0,seed);
+	lfo3.set_WaveType(LFO_RANDOM);
+	lfo3.set_freq(0.61);
+	lfo3.set_mod(1.0);
+	lfo4.init((double)inprops.srate,0.0,seed);
+	lfo4.set_WaveType(LFO_RANDOM);
+	lfo4.set_freq(0.66);
+	lfo4.set_mod(1.0);
+	vcomb.init(double)srate,0.8);
+	double delay1,delay2,delay3,delay4;
+	double delaymod = 0.05;
+	delay1 = 0.05;
+	delay2 = 0.056;
+	delay3 = 0,061;
+	delay4 = 0.068;
+	vcomb.init(double)srate,delay4 + (delay4 * delaymod));
+	vcomb.setfgains(
+#endif
+	//stopwatch(1);
+	//run main single-sample loop
+	float l_ip,r_ip,l_out,r_out,l_direct,r_direct,mono_direct;
+	for(;;){
+		int rv;
+		float out,earlies;
+		//read mon/left channnel
+		if((rv = fgetfloatEx(&l_ip,ifd,0)) < 0){
+			fprintf(stderr,"\nreverb: error reading infile"); 
+			exit(1);
+		}
+		if(!rv) break;	   // end of inputfile	- handle any trailertime
+		//apply any conditioning to direct signal
+		if(tc_l_lowcut)
+			l_ip = (float) tc_l_lowcut->tick((double)l_ip);
+		if(l_highcut)
+			l_ip = (float) l_highcut->tick((double)l_ip);
+		l_direct = l_ip;
+		mono_direct = l_direct;
+		r_direct = l_direct;
+		//handle R channnel if active
+		if(chans==2){
+			if((rv = fgetfloatEx(&r_ip,ifd,0)) < 0){
+				fprintf(stderr,"\nError in reading file"); 
+				exit(1);			
+			}		
+			if(!rv) break;					
+			//apply any conditioning to direct signal
+			if(tc_r_lowcut)
+				r_ip = (float) tc_r_lowcut->tick((double)r_ip);
+			if(r_highcut)
+				r_ip = (float) r_highcut->tick((double)r_ip);
+			r_direct = r_ip;
+			if(want_mono)
+				mono_direct = (l_direct + r_direct) * 0.5f;			
+		}
+		//mono_direct = (mixed) input to reverberator
+		//possibly also filter it...
+		if(lp)
+			mono_direct = (float) lp->tick(mono_direct);					//input lowpass
+		earlies = tdelay->tick(mono_direct);				// ---> early reflections 	
+		//send (filtered) mono output from reverberator
+		out = earlies + moorerverb->tick(earlies + mono_direct) * diff_gain;
+		//and send thru multi_channel diffusers
+		//NB must handle first 2 chans as special cases, then loop thru the rest
+		//l_out = want_mono ? mono_direct : l_direct;
+		//l_out *= dry_gain;
+		l_out = chan_ap[0]->tick(out) * wet_gain;
+		if(want_mono)			
+			l_out += (mono_direct * dry_gain);
+		else 
+			l_out += l_direct * dry_gain;
+		
+		
+		if((float)fabs((double)l_out) > peaks[0].value) {
+			peaks[0].value = (float)fabs((double)l_out);
+			peaks[0].position = outsamps;
+		}
+		if(fputfloatEx(&l_out,ofd) < 1){
+			fprintf(stderr,"\nreverb: error writing output file"); 
+			exit(1);
+		}
+ 
+		if(outchans>=2){
+			r_out = r_direct * dry_gain;
+			r_out += chan_ap[1]->tick(out) * wet_gain;			
+			if((float)fabs((double)r_out) > peaks[1].value) {
+				peaks[1].value = (float)fabs((double)r_out);
+				peaks[1].position = outsamps;
+			}
+			if(fputfloatEx(&r_out,ofd) < 1){
+				fprintf(stderr,"\nreverb: error writing output file"); 
+				exit(1);
+			}
+		}
+
+		//now any further channels
+		for(i=2;i < outchans; i++){						
+			l_out = chan_ap[i]->tick(out) * wet_gain;		   			
+			if((float)fabs((double)l_out) > peaks[i].value) {
+				peaks[i].value = (float)fabs((double)l_out);
+				peaks[i].position = outsamps;
+			}
+			if(fputfloatEx(&out,ofd) < 1){
+				fprintf(stderr,"\nreverb: error writing output file"); 
+				exit(1);
+			}
+		}
+		outsamps++;
+		if((outsamps % step) == 0)
+			//inform(step,sr);
+            fprintf(stdout,"%.2f\r",(double)outsamps/(double)sr);
+
+	}
+	//end of infile; play out reverb tail as much as we're allowed
+	int trtime  = (int)( trailertime * sr);	
+	for(i = 0; i < 	trtime; i++){
+		float tr_l_ip,tr_r_ip = 0.0f, tr_out,tr_earlies,tr_l_direct,tr_mono_direct,tr_r_direct,
+			tr_l_out,tr_r_out;
+		tr_l_ip = 0.0;
+		if(tc_l_lowcut)
+			tr_l_ip = (float) tc_l_lowcut->tick(tr_l_ip);
+		if(l_highcut)
+			tr_l_ip = (float) l_highcut->tick(tr_l_ip);
+		tr_l_direct = tr_l_ip;
+		tr_mono_direct = tr_l_direct;
+		tr_r_direct = tr_l_direct;
+		//handle R channnel if active
+		if(chans==2){
+			tr_r_ip = 0.0;	   // inout is zero, except if using input filters
+			// get any samples from input conditioning filters
+			if(tc_r_lowcut)
+				tr_r_ip = (float) tc_r_lowcut->tick(tr_r_ip);						   
+			if(r_highcut)
+				tr_r_ip = (float) r_highcut->tick(tr_r_ip);
+			tr_r_direct = tr_r_ip;
+			if(want_mono)
+				tr_mono_direct = (tr_l_direct + tr_r_direct) *0.5f;			
+		}
+		if(lp)
+			tr_mono_direct = (float) lp->tick(tr_mono_direct);
+		tr_earlies = tdelay->tick(tr_mono_direct);
+		//send (filtered) mono output from reverberator
+		tr_out = tr_earlies +  moorerverb->tick(tr_earlies * tr_mono_direct) * diff_gain;
+		tr_l_out = chan_ap[0]->tick(tr_out) * wet_gain;
+		if(want_mono)			
+			tr_l_out += (tr_mono_direct * dry_gain);
+		else 
+			tr_l_out += tr_l_direct * dry_gain;
+		if((float)fabs((double)tr_l_out) > peaks[0].value) {
+			peaks[0].value = (float)fabs((double)tr_l_out);
+			peaks[0].position = outsamps;
+		}
+		if(fputfloatEx(&tr_l_out,ofd) < 1){
+			fprintf(stderr,"\nreverb: error writing to output file"); 
+			exit(1);
+		}
+		if(outchans>=2){
+			tr_r_out = tr_r_direct * dry_gain;
+			tr_r_out += chan_ap[1]->tick(tr_out) * wet_gain;	
+			if((float)fabs((double)tr_r_out) > peaks[1].value) {
+				peaks[1].value = (float)fabs((double)tr_r_out);
+				peaks[1].position = outsamps;
+			}
+			if(fputfloatEx(&tr_r_out,ofd) < 1){
+				fprintf(stderr,"\nreverb: error writing to output file"); 
+				exit(1);
+			}
+		}
+
+		for(int j = 2;j < outchans;j++){
+			float ch_out;
+			ch_out = chan_ap[j]->tick(tr_out)  * wet_gain;
+			if((float)fabs((double)ch_out) > peaks[j].value) {
+				peaks[j].value = (float)fabs((double)ch_out);
+				peaks[j].position = outsamps;
+			}
+			if(fputfloatEx(&ch_out,ofd) < 1){
+				fprintf(stderr,"\nreverb: error writing to output file"); 
+				exit(1);
+			}
+		}
+		outsamps++;
+		if((outsamps % step) ==0)
+			//inform(step,sr);
+            fprintf(stdout,"%.2f\r",(double)outsamps/(double)sr);
+	}
+	//stopwatch(0);
+	printf("\nPEAK data:\n");
+	for(i=0;i < outchans;i++)
+		printf("Channel %d: %.4f at frame %d: %.4lf secs\n",i+1,
+			peaks[i].value,peaks[i].position, (double) peaks[i].position / (double)props.srate);
+	sndputpeaks(ofd,outchans,peaks);
+	sndcloseEx(ifd);
+	
+	if(sndcloseEx(ofd) < 0)
+		fprintf(stderr,"\nwarning: error closing outfile");
+	delete [] peaks;
+	if(moorerverb)
+		delete moorerverb;
+	for(i=0;i < outchans;i++)
+		delete chan_ap[i];
+	delete [] chan_ap;
+	delete tdelay;
+	if(usertaps)
+		delete [] ptaps;
+	return 0;
+}
+
+
+
+void
+usage(void)
+{
+	fprintf(stderr,"\nreverb: Multi-channel reverberator\n\tVersion 1.1 Release 4 CDP 1998,1999\n");
+	fprintf(stderr,
+	"usage: reverb [flags] infile outfile rgain mix rvbtime absorb\n"
+	"                                 lpfreq trailertime [times.txt]\n"
+	"flags    : any or all of the following:\n"
+	"   -cN   : create outfile with N channels (1 <= N <= 16 : default =  2)\n"
+	"   -eFILE: read early reflections from breakpoint textfile FILE\n"
+	"   -f    : force floating-point output (default: infile format)\n"
+	"   -HN   : apply Highcut filter with cutoff freq N Hz to infile (6db/oct)\n"
+	"   -LN   : apply Lowcut filter with cutoff freq N Hz to infile	(12dB/oct)\n"
+	"   -pN   : force reverb predelay to N msecs (shifts early reflections)\n"
+	"rgain    : set level of dense reverb (0.0 <= rgain <= 1.0)\n"    
+	"mix      : dry/wet balance (source and reverb): 1.0<-- mix -->= 0.0\n"
+	"rvbtime  : reverb decay time (to -60dB) in seconds\n"
+	"absorb   : degree of hf damping to suggest air absorption: 0.0 <= absorb <= 1.0\n"
+	"lpfreq   : lowpass filter cutoff freq (Hz) applied at input to reverb\n" 
+	"           to disable either filter (absorb, lpfreq), use the value 0\n" 
+	"trtime   : trailer time added to outfile for reverb tail (secs)\n"
+	"times.txt: list of delay times (msecs) for 6 comb and 4 allpass filters\n"
+	);
+	exit(0);
+}
+
+long readtaps(FILE *fp, deltap **taps,double sr)
+{
+	vector<deltap> vtaps;
+	deltap thistap;
+	char line[256];
+	int ntaps = 0;
+	int linecount = 1;
+	bool error = false;
+	double time= 0.0, current_time = 0.0, val = 0.0;
+
+	double sample_duration = 1.0 / sr;
+	if(fp==0 || sr <= 0.0){
+		*taps = 0;
+		return 0;
+	}
+	while(fgets(line,256,fp))	{
+		int got;
+		if(line[0] == '\n' || line[0]== ';'){	//allow comment lines; this does not check for leading white-space...
+			linecount++;
+			continue;
+		}
+
+		if((got = sscanf(line,"%lf%lf",&time,&val))<2){			
+			fprintf(stderr,"\nerror in reflections file: line %d",linecount);
+			error =true;
+			break;
+		}
+		if(time < 0.0){
+			fprintf(stderr,"\nerror in tapfile: line %d: bad time value",linecount);
+			error = true;
+			break;
+		}
+		//if (first) val = 0.0, or VERY close
+		if(time < sample_duration)	{ //non-zero time must be at least one sample on!			
+			fprintf(stderr,"\nWARNING: very small taptime %.4lf treated as zero",time);
+			time = 0.0;			
+		}
+
+		if(current_time != 0.0 && time==current_time){
+			fprintf(stderr,"\nWARNING: duplicate time in line %d: ignored",linecount);
+			linecount++;
+			continue;
+		}
+		if(time < current_time){
+			fprintf(stderr,"\nerror in tapfile: time out of sequence: line %d",linecount);
+			error = true;
+			break;
+		}
+		current_time = time;
+		thistap.pos = time;
+		if(val <=0.0 || val > 1.0){
+			fprintf(stderr,"\nerror: bad amplitude in tapfile: line %d",linecount);
+			error = true;
+			break;
+		}
+		thistap.val = val;
+		vtaps.push_back(thistap);
+		linecount++;
+	}
+	
+	ntaps = vtaps.size();
+	if(ntaps==0){
+		fprintf(stderr,"\ntapfile contains no data!");
+		error = true;
+	}
+	if(error){
+		*taps = 0;
+		return 0;
+	}
+	*taps = new deltap[ntaps];
+	if(taps==0)
+		return 0;
+	vector<deltap>::iterator I =  vtaps.begin(); 
+	for(int i = 0; I != vtaps.end();i++, I++){
+		(*taps)[i].pos = I->pos;
+		(*taps)[i].val = I->val;
+	}
+
+	return ntaps;
+}

dev/externals/reverb/reverberator.cpp

@@ -0,0 +1,1488 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+// reverberator.cpp: implementation of the reverberator class.
+// see also tapdelay.cpp
+//TODO: simple delay for use with tapped delayline
+// NB blame VC6 for horrible code, new returns null...
+//////////////////////////////////////////////////////////////////////
+#include <stdio.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <math.h>
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+#include "reverberator.h"
+#ifndef PI
+# define PI 3.141592653589793238462643
+#endif
+#ifndef TWOPI
+#define TWOPI (2.0 * PI)
+#endif
+
+#define ROOT2O2 (0.7071067811965475244)
+#define SPN 1.65436e-24			   //RWD: use a proper MSVC version...?
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+//for example: 
+//const deltap taps[4] = {{0.1,0.9},{0.25,.75},{0.5,.5},{0.618,.25}};
+//NB this accepts a zero taptime
+
+tapdelay::tapdelay(const deltap *taps,unsigned int ctaps,double sr) : dtaps(taps)
+{
+	tapbuf = 0;
+	ntaps = ctaps;
+	buflen = 0;
+	reader = 0;
+	tapptr = 0;
+	tapcoeff = 0;
+	srate = sr;
+	output_gain = 1.0;
+	max_gain = 1.0;
+}
+#define SAFETY (2)
+bool tapdelay::create(void)
+{
+#ifdef _DEBUG
+	assert(tapbuf==0);
+	assert(ntaps > 0);
+#endif
+	if(tapbuf)
+		return false;				//or some informative errorval...
+	if(ntaps <=0)
+		return false;
+	buflen = (unsigned int)((dtaps[ntaps-1].pos * srate) + SAFETY);
+	tapbuf = new double[buflen];
+	if(!tapbuf)
+		return false;
+	memset(tapbuf,0,buflen*sizeof(double));
+	
+	tapptr = new unsigned int[ntaps];
+	if(!tapptr){
+		delete [] tapbuf; tapbuf = 0;
+		return false;
+	}
+	tapcoeff = new double[ntaps];	
+	if(!tapcoeff) {
+		delete [] tapbuf; tapbuf = 0;
+		delete [] tapptr; tapptr = 0;
+		return false;
+	}
+	//copy taps data, and get scalefac from total amps
+	double sum = 0.0;
+	for(unsigned int i = 0;i < ntaps;i++){
+		//RWD.10.98 allow zero taptime, for input mix control...
+		if(dtaps[i].pos == 0.0)
+			tapptr[i] = 0;
+		else
+			tapptr[i] = buflen - (int)(dtaps[i].pos * srate);	
+			
+		tapcoeff[i] = dtaps[i].val;
+		sum += fabs(tapcoeff[i]);
+#ifdef _DEBUG
+		assert(tapptr[i] < buflen);
+#endif
+	}
+	//some idiot might give us empty taps! */
+	if(sum == 0.0){
+		delete [] tapbuf; tapbuf = 0;
+		delete [] tapptr; tapptr = 0;
+		delete [] tapcoeff; tapcoeff = 0;
+		return false;
+	}
+	max_gain = (double)ntaps / sum;
+	one_over_ntaps = 1.0 / (double)ntaps;
+	return true;
+}
+
+double tapdelay::get_maxgain(void)
+{
+	return max_gain;
+}
+
+void tapdelay::set_gain(double gain)
+{
+	output_gain = gain;
+}
+
+double tapdelay::tick(double insam)
+{
+	double val = 0.0;
+	tapbuf[reader++] = insam;
+	if(reader == buflen)
+		reader = 0;
+	for(unsigned int i = 0;i < ntaps;i++) {
+		val += 	((tapbuf[tapptr[i]++]) * (tapcoeff[i]));
+		if(tapptr[i] == buflen)
+			tapptr[i] = 0;
+	}
+	val *= one_over_ntaps;
+	return val * output_gain;
+
+}
+
+tapdelay::~tapdelay()
+{
+	delete [] tapbuf;
+	delete [] tapptr;
+	delete [] tapcoeff;
+}
+
+
+delay::delay(unsigned int len, double gain)
+{
+	delbuf= 0;
+	dgain = gain;
+	ptr = 0;
+	buflen = len;
+}
+
+delay::~delay()
+{
+	delete delbuf;
+}
+
+bool delay::create(void)
+{
+	 if(buflen <= 0)
+		 return false;
+	 if(delbuf)
+		 return false;
+
+	 delbuf = new double[buflen];
+	 if(!delbuf) {		
+		 return false;
+	 }
+	 memset(delbuf,0,buflen*sizeof(double));
+	 return true;
+}
+
+double delay::tick(double insam)
+{
+	double val;
+	val = delbuf[ptr];
+	delbuf[ptr++] = insam;
+	if(ptr==buflen)
+		ptr = 0;
+	return  val*dgain;
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+// allpassfilter 
+//////////////////////////////////////////////////////////////////////
+//global helper functions: really belong in a new cdplib...except it's nice to have proper bool retval
+//straight from Csound: nreverb (vdelay.c)
+static bool prime( int val )
+{
+    int i, last;
+	
+    last = (int)sqrt( (double)val );
+    for ( i = 3; i <= last; i+=2 ) {
+      if ( (val % i) == 0 ) return false;
+    }
+    return true;
+}
+
+unsigned int to_prime(double dur,double sr)
+{
+		unsigned int time = (unsigned int) (dur * sr);
+		if(time % 2== 0) time += 1;
+		while(!prime(time)) 
+			time +=2;
+		return time;
+}
+#define VMODMAX (0.5)
+//standard allpass, with optional pre_delay
+allpassfilter::allpassfilter(long sr,unsigned int buflen, double ingain,unsigned int predelay = 0)
+{
+    vmod                = 0.1;
+    vfreq               = 0.5; 
+	rvbuf1				= 0;
+	pre_dline			= 0;
+	writer1 = reader1	= 0;
+    sinelfo             = 0;
+    noiselfo            = 0;
+    srate               = sr;
+	rvblen				=  buflen + (unsigned int)(buflen * VMODMAX);
+	gain				= ingain;
+	pre_delay			= predelay;
+    
+}
+
+bool allpassfilter::create(void)
+{
+	if((gain < -1.0) || (gain > 1.0))		//do I want negative gain vals?
+		return false;
+	rvbuf1 = new double[rvblen];	
+	if(!rvbuf1 ){
+#ifdef _DEBUG
+		printf("\nallpass: failed to allocate buffer size %d",rvblen);
+#endif		
+		return false;
+	}
+				
+	if(pre_delay > 0){
+		pre_dline = new delay(pre_delay,1.0);
+		if(!pre_dline->create()){
+#ifdef _DEBUG
+		    printf("\nallpass: can't create pre_dline" );
+#endif		
+		    delete []rvbuf1; rvbuf1 = 0; 
+		    return false;
+	    }
+	}
+	memset(rvbuf1,0,rvblen*sizeof(double));
+    sinelfo = new fastlfo();
+    sinelfo->init((double) srate,0.0,0,1);
+    sinelfo->set_WaveType(LFO_SINE);
+    sinelfo->set_freq(vfreq);
+    sinelfo->set_mod(1.0);
+    noiselfo = new fastlfo();
+    noiselfo->init((double) srate,0.0,0,1);
+    noiselfo->set_WaveType(LFO_RAND_GAUSS);
+    noiselfo->set_freq(vfreq * 2.5);
+    noiselfo->set_mod(1.0);
+
+	return true;
+}
+
+void allpassfilter::set_vmod(double amount)
+{
+    double numod = VMODMAX;
+    if(amount < numod)
+        numod = amount;
+    vmod = numod;
+}
+
+void allpassfilter::set_vfreq(double freq)
+{
+    vfreq = freq;
+    sinelfo->set_freq(vfreq);
+    noiselfo->set_freq(vfreq * 2.5);
+}
+
+/* TODO: optimize this:  separate  funcs for pre-dline and plain */
+double allpassfilter::tick(double insamp)
+{
+	double output, input; 
+	input = insamp;
+	output = (-gain) * input;
+	output += rvbuf1[reader1++];
+	input += gain * output;
+
+	if(pre_dline)
+		rvbuf1[writer1++] = pre_dline->tick(input);
+	else
+		rvbuf1[writer1++] = input;
+	if(reader1 == rvblen)
+		reader1 = 0;
+	if(writer1 == rvblen)
+		writer1 = 0;
+
+	return output;
+}
+/* TODO: implement this! */
+double allpassfilter::vtick(double insamp)
+{
+	double output, input; 
+	input = insamp;
+	output = (-gain) * input;
+	output += rvbuf1[reader1++];
+	input += gain * output;
+
+	if(pre_dline)
+		rvbuf1[writer1++] = pre_dline->tick(input);
+	else
+		rvbuf1[writer1++] = input;
+	if(reader1 == rvblen)
+		reader1 = 0;
+	if(writer1 == rvblen)
+		writer1 = 0;
+
+	return output;
+}
+
+
+
+allpassfilter::~allpassfilter()
+{
+	delete [] rvbuf1;
+	delete pre_dline;
+}
+
+//////////////////////////////////////////////////////////////////////
+//// nested_allpass
+//overall delay = ap1.pre_delay + ap2.pre_deay + ap1.length + ap2.length + post_delay
+nested_allpass::nested_allpass(double srate,double lpfreq, unsigned int outertime, unsigned int time1, unsigned int time2,
+							   double gain, double gain1, double gain2, 
+							   unsigned int delay1 = 0, unsigned int delay2 = 0)
+							   
+{
+	outer_gain		= gain;
+	outer_time		= outertime;
+	ap1_gain		= gain1;
+	ap2_gain		= gain2;	
+	ap1_length		= time1;
+	ap2_length		= time2;
+	ap1_delay		= delay1;
+	ap2_delay		= delay2;
+	buf				= 0;
+	if(ap2_gain == 0.0)
+		ap2_length = 0;
+	ap1				= ap2	= 0;		
+	writer = reader = 0;	
+	lp = 0;
+	lpfreq_ = lpfreq;
+	sr_ = srate;
+
+}
+
+nested_allpass::~nested_allpass()
+{	
+	delete [] buf;
+	delete ap1;
+	delete ap2;
+	//RWD
+	delete lp;
+}
+
+bool nested_allpass::create(void)
+{
+	if(outer_gain < -1.0 || outer_gain > 1.0)
+		return false;	
+	if(outer_time <= 0)
+		return false;
+	buf = new double[outer_time];
+	if(!buf)
+		return false;
+	memset(buf,0,outer_time * sizeof(double));
+	if(ap1_length > 0 && ap1_gain > 0.0){
+		ap1 = new allpassfilter((long)sr_,ap1_length,ap1_gain,ap1_delay);
+		if(!ap1->create()){
+			delete [] buf; buf = 0;
+			return false;
+		}
+	}
+	
+	if(ap2_length >0 && ap2_gain > 0.0){
+		ap2 = new allpassfilter((long)sr_,ap2_length,ap2_gain,ap2_delay);
+		if(!ap2->create()){
+			delete ap1;	 ap1 = 0;
+			delete [] buf;  buf = 0;
+			return false;
+		}
+	}
+	// RWD experiment: internal absorption damping
+	lp = new onepole(lpfreq_,sr_,LOW_PASS);
+	if(!lp){
+		delete ap1;	 ap1 = 0;
+		delete ap2;	 ap2 = 0;
+		delete [] buf;  buf = 0;
+		return false;
+	}
+	damping_ = false;
+	return true;
+}
+
+//messy casts - sort this out later...
+double nested_allpass::tick(double insam)
+{
+
+	double output, input,nest_out; 
+	input = insam;
+	output = (-outer_gain) * input;
+	output += buf[reader++];
+	input += outer_gain * output;
+	//input is to two allpasses, and the post_delay, which acts as the overall allpass
+	if(!ap2) {	//no second allpass, may be one or no first allpass
+		if(ap1)			
+			nest_out =  ap1->tick(input);
+		else
+			//behave as normal allpass
+			nest_out =  input;
+	}
+	else		
+		nest_out = 	 ap2->tick(ap1->tick(input));
+
+	//RWD experiment!
+	if(damping_)
+	   buf[writer++] = lp->tick(nest_out);
+	else
+		buf[writer++] = nest_out;
+	//wrap indices
+	if(writer == outer_time)
+		writer = 0;
+	if(reader == outer_time)
+		reader = 0;
+			
+	return output;
+}
+
+
+//////////////////////////////////////////////////////////////////////
+// lpcomb 
+//////////////////////////////////////////////////////////////////////
+
+lpcomb::lpcomb(unsigned int buflen, double fbgain, double lpcoeff, double pre = 1.0)
+{
+	combbuf = 0;
+	lp = 0;
+	gain = fbgain;
+	lpgain = lpcoeff;
+	rvblen = buflen;
+	prescale  = pre;
+	reader1 = writer1 = 0;
+}
+
+bool lpcomb::create(void)
+{
+	if(prescale <= 0.0) 
+		return false;
+	if((gain < -1.0 ) || (gain > 1.0))		//do I want negative gain vals?
+		return false;
+
+	combbuf = new double[rvblen];
+	if(!combbuf)
+		return false;
+	lp = new lowpass1(lpgain);
+	if(!lp){
+		delete [] combbuf; combbuf = 0;
+		return false;
+	}
+	memset(combbuf,0,rvblen*sizeof(double));
+	return true;
+}
+
+ //wet output only, unscaled
+double lpcomb::tick(double insam)
+{
+	double input,output,lpsig;
+	input = insam * prescale;
+	
+	output = combbuf[reader1++];		 //do gain scaling outside
+	lpsig = lp->tick(output);
+
+	combbuf[writer1++] = input + (gain * lpsig) ;	//feedback scaled output
+	
+
+	if(writer1 == rvblen)
+		writer1 = 0;
+	if(reader1 == rvblen)
+		reader1 = 0;
+	return  output;
+}
+
+lpcomb::~lpcomb()
+{
+	delete [] combbuf;
+	delete lp;
+}
+ 
+ //////////////////////////////////////////////////////////////////////
+// lowpass 	: perhaps a bit ott as a class, but it allows fairly independent fiddling
+//////////////////////////////////////////////////////////////////////
+lowpass1::lowpass1(double filtgain)
+{
+	temp = output = 0.0;
+	gain = filtgain;
+}
+
+lowpass1::~lowpass1()
+{
+}
+ //IIR lowpass
+inline double lowpass1::tick(double insam)
+{	 
+	output  = temp;
+	temp = insam + (output * gain);
+	return output;
+
+}
+// NB allows zero frequency: straight pass-through
+onepole::onepole(double freq,double srate, int low_high)
+{
+	if(freq == 0.0){
+		a1 = 1.0;
+		a2 = 0.0;
+	}
+	else {
+		b = 2.0 - cos(TWOPI * (freq/srate));
+		a2 = sqrt((b*b) - 1.0) - b;
+		a1 = 1.0 + a2;
+		if(low_high==HIGH_PASS)
+			a2 = -a2;
+	}
+	lastout = 0.0;
+}
+
+onepole::~onepole()
+{
+	 //nothing to do
+}
+
+double onepole::tick(double input)
+{
+	double output;
+	output  =  a1*input - a2*lastout;
+	lastout =  output;
+	return output;
+}
+
+//uses code from CDP eq.c
+
+tonecontrol::tonecontrol(double frq,double dbfac,int type,double srate)
+{
+	freq = frq;
+	boost = dbfac;
+	sr = srate;
+	tc_type = type;
+	x1 = x2 = y1 = y2 = 0.0;
+	a0 = a1 = a2 = b1 = b2 = 0.0;
+}
+
+tonecontrol::~tonecontrol()
+{
+   //nothing to do;
+}
+
+bool tonecontrol::create(void)
+{
+	if(sr <= 0.0)
+		return false;
+	if(freq <= 0.0)
+		return false;
+	//normalize freq
+	freq /= sr;
+	switch(tc_type){
+	case(LOW_SHELVE):
+		lshelve();
+		break;
+	case(HIGH_SHELVE):
+		hshelve();
+		break;
+	default:
+		return false;
+		break;
+	}
+	return true;
+}
+
+double tonecontrol::tick(double input)
+{
+	double out;
+
+	if(boost== 0.0)
+		return input;
+
+	out = ((a0 * input) + (a1 * x1) + (a2 * x2) - (b1 * y1) - (b2 * y2));
+	x2 = x1;
+	x1 = input;
+	y2 = y1;
+	y1 = out;
+	return out;
+}
+
+void tonecontrol::lshelve(void)
+{
+	double a, A, F, tmp, b0, recipb0, asq, F2, gamma2, siggam2, gam2p1;
+	double gamman, gammad, ta0, ta1, ta2, tb0, tb1, tb2, aa1, ab1;
+
+	a = tan(PI * (freq - 0.25));	/* Warp factor */ 
+	asq = a * a;
+	A = pow(10.0, boost/20.0);	/* Cvt dB to factor */
+	if((boost < 6.0) && (boost > -6.0)) F = sqrt(A);
+	else if (A > 1.0) F = A/sqrt(2.0);
+	else F = A * sqrt(2.0);
+	  /* If |boost/cut| < 6dB, then doesn't make sense to use 3dB point.
+	     Use of root makes bandedge at half the boost/cut amount
+	  */
+
+	F2 = F * F;
+	tmp = A * A - F2;
+	if(fabs(tmp) <= SPN) gammad = 1;
+	else gammad = pow( (F2 - 1)/ tmp, 0.25);	/* Fourth root */
+	gamman = sqrt(A) * gammad;
+
+   /* Once for the numerator */
+
+	gamma2 = gamman * gamman;
+	gam2p1 = 1 + gamma2;
+	siggam2 = 2 * ROOT2O2 * gamman;
+	
+	ta0 = gam2p1 + siggam2;
+	ta1 = -2 * (1 - gamma2);
+	ta2 = gam2p1 - siggam2;
+
+   /* And again for the denominator */
+
+	gamma2 = gammad * gammad;
+	gam2p1 = 1 + gamma2;
+	siggam2 = 2 * ROOT2O2 * gammad;
+	
+	tb0 = gam2p1 + siggam2;
+	tb1 = -2 * (1 - gamma2);
+	tb2 = gam2p1 - siggam2;
+
+   /* Now bilinear transform to proper centre frequency */
+
+	aa1 = a * ta1;
+	a0 = ta0 + aa1 + asq * ta2;
+	a1 = 2 * a * (ta0 + ta2) + (1 + asq) * ta1;
+	a2 = asq * ta0 + aa1 + ta2;
+
+	ab1 = a * tb1;
+	b0 = tb0 + ab1 + asq * tb2;
+	b1 = 2 * a * (tb0 + tb2) + (1 + asq) * tb1;
+	b2 = asq * tb0 + ab1 + tb2;
+	
+   /* Normalise b0 to 1.0 for realisability */
+
+	recipb0 = 1 / b0;
+	a0 *= recipb0;
+	a1 *= recipb0;
+	a2 *= recipb0;
+	b1 *= recipb0;
+	b2 *= recipb0;
+
+}
+
+void tonecontrol::hshelve(void)
+{
+	double a, A, F, tmp, b0, recipb0, asq, F2, gamma2, siggam2, gam2p1;
+	double gamman, gammad, ta0, ta1, ta2, tb0, tb1, tb2, aa1, ab1;
+
+	a = tan(PI * (freq - 0.25));	/* Warp factor */ 
+	asq = a * a;
+	A = pow(10.0, boost/20.0);	/* Cvt dB to factor */
+	if(boost < 6.0 && boost > -6.0) F = sqrt(A);
+	else if (A > 1.0) F = A/sqrt(2.0);
+	else F = A * sqrt(2.0);
+	  /* If |boost/cut| < 6dB, then doesn't make sense to use 3dB point.
+	     Use of root makes bandedge at half the boost/cut amount
+	  */
+
+	F2 = F * F;
+	tmp = A * A - F2;
+	if(fabs(tmp) <= SPN) gammad = 1;
+	else gammad = pow( (F2 - 1)/ tmp, 0.25);	/* Fourth root */
+	gamman = sqrt(A) * gammad;
+
+   /* Once for the numerator */
+
+	gamma2 = gamman * gamman;
+	gam2p1 = 1 + gamma2;
+	siggam2 = 2 * ROOT2O2 * gamman;
+	
+	ta0 = gam2p1 + siggam2;
+	ta1 = -2 * (1 - gamma2);
+	ta2 = gam2p1 - siggam2;
+
+	ta1 = -ta1;
+	
+   /* And again for the denominator */
+
+	gamma2 = gammad * gammad;
+	gam2p1 = 1 + gamma2;
+	siggam2 = 2 * ROOT2O2 * gammad;
+	
+	tb0 = gam2p1 + siggam2;
+	tb1 = -2 * (1 - gamma2);
+	tb2 = gam2p1 - siggam2;
+
+	tb1 = -tb1;
+
+   /* Now bilinear transform to proper centre frequency */
+
+	aa1 = a * ta1;
+	a0 = ta0 + aa1 + asq * ta2;
+	a1 = 2 * a * (ta0 + ta2) + (1 + asq) * ta1;
+	a2 = asq * ta0 + aa1 + ta2;
+
+	ab1 = a * tb1;
+	b0 = tb0 + ab1 + asq * tb2;
+	b1 = 2 * a * (tb0 + tb2) + (1 + asq) * tb1;
+	b2 = asq * tb0 + ab1 + tb2;
+	
+   /* Normalise b0 to 1.0 for realisability */
+
+	recipb0 = 1 / b0;
+	a0 *= recipb0;
+	a1 *= recipb0;
+	a2 *= recipb0;
+	b1 *= recipb0;
+	b2 *= recipb0;
+
+}
+
+
+/***************** GARDNER DIFFUSERS ***************/
+
+// NOTE: the published Gardner designs are here modified by adding a onepole lp filter
+// in the internal outer feedback loop of each nested allpass filter. 
+// When combined with the overall feedback filter, this will of course reduce the reverb decay time.
+// The option exists to bypass this overall feedback filter, and rely solely on the internal ones.
+
+
+// NOTE: to get an infinite reverb, all feedback filters have to be bypassed by setting lp_freq to 0.0
+// otherwise, the output will still decay eventually !
+
+small_diffuser::small_diffuser(unsigned int pre_delay, 
+							   const NESTDATA *apdata1,
+							   const NESTDATA *apdata2,double fb_gain,double lp_freq,double srate){
+
+	  ap1_data = *apdata1;
+	  ap2_data = *apdata2;
+	  predelay_time = pre_delay;
+	  lpgain = fb_gain;
+	  lpfreq = lp_freq;
+	  ap1 = ap2 = 0;
+	  lp1 = 0;
+	  predelay = 0;
+	  out1 = out2 = 0.0;
+	  sr = srate;
+	  damping_ = false;
+}
+
+small_diffuser::~small_diffuser()
+{
+	delete ap1;
+	delete ap2;
+	delete predelay;	
+	delete lp1;
+}
+
+bool small_diffuser::create(void)
+{
+	if(ap1_data.time1 == 0 || ap1_data.gain1 <= 0.0){
+#ifdef _DEBUG
+		printf("\ndiffuser: bad parameters(1)");
+#endif		
+		return false;
+	}
+	if(ap1_data.time2 ==0 || ap1_data.gain2 <= 0.0){
+#ifdef _DEBUG
+		printf("\ndiffuser: bad parameters(2)");
+#endif		
+		return false;
+	}
+	if(ap2_data.time1 == 0 || ap2_data.gain1 < 0.0){
+#ifdef _DEBUG
+		printf("\ndiffuser: bad parameters(3)");
+#endif		
+		return false;
+	}
+	if(ap2_data.time2 ==0 || ap2_data.gain2 < 0.0){
+#ifdef _DEBUG
+		printf("\ndiffuser: bad parameters(4)");
+#endif		
+		return false;
+	}
+
+	if(sr <=0.0){
+#ifdef _DEBUG
+		printf("\ndiffuser: bad srate parameter)");
+#endif		
+		return false;
+	}
+	if(lpfreq <0.0){
+#ifdef _DEBUG
+		printf("\ndiffuser: bad freq parameter)");
+#endif		
+		return false;
+	}
+
+	ap1 = new nested_allpass(sr,lpfreq,ap1_data.time1,ap1_data.time2,ap1_data.time3,
+							ap1_data.gain1,ap1_data.gain2,ap1_data.gain3);
+	if(!ap1->create()){
+#ifdef _DEBUG
+		printf("\ndiffuser: can't create first allpass");
+#endif		
+		return false;
+	}
+	if(ap2_data.gain1 != 0.0){	 //allow ap to eliminate second block
+		ap2 = new nested_allpass(sr,lpfreq,ap2_data.time1,ap2_data.time2,ap2_data.time3,
+							ap2_data.gain1,ap2_data.gain2,ap2_data.gain3);
+		if(!ap2->create()){
+#ifdef _DEBUG
+			printf("\ndiffuser: can't create second allpass");
+#endif		
+			delete ap1; ap1 = 0;
+			return false;
+		}
+	}
+	if(predelay_time > 0){
+		predelay = new delay(predelay_time,1.0);
+		if(!predelay->create()){
+#ifdef _DEBUG
+			printf("\ndiffuser: can't create predelay");
+#endif		
+			delete ap1;	ap1 = 0;
+			delete ap2;	ap2 = 0;
+			return false;
+		}
+	}
+	// NB if lpfreq == 0, onepole is no-op - returns input
+	lp1 = new onepole(lpfreq,sr,LOW_PASS);
+	if(lp1 == 0){
+		delete [] predelay; predelay = 0;
+		delete ap1;	ap1 = 0;
+		delete ap2;	ap2 = 0;
+		return false;
+	}
+	ap1->set_damping(damping_);
+	ap2->set_damping(damping_);
+	return true;
+}
+
+
+double small_diffuser::tick(double insam) 
+{
+	double filter_out;	
+	double lp_in;	
+	double output;
+	double ip;
+	// may only contain one nested allpass filter:
+	// input to filter is either out1 or out2
+	if(ap2)
+		lp_in = out2;	   //= both in series
+	else
+		lp_in = out1;
+	// RWD NOTE: if the damping_ mechanism is preferred, this filter can be omitted and ...
+	filter_out =  lp1->tick(lp_in);
+	// ... just use this				 
+	//	filter_out = lp_in ;
+	
+	ip = insam + lpgain * filter_out;	
+	if(predelay)
+		ip = predelay->tick(ip);
+	out1 = ap1->tick( ip);
+	if(ap2){
+		out2 = ap2->tick(out1);
+		output =  (out1 + out2)  * 0.5;
+	}
+	else
+		output =  out1;
+	return  output;
+}
+
+//post-delay almost certainly doe not need to be prime...
+medium_diffuser::medium_diffuser(double post_delay,const NESTDATA *apdata1,
+					const NESTDATA *apdata2,
+					double gain,
+					double lp_freq,
+					double srate)
+{
+	ap1_data = *apdata1;
+	ap2_data = *apdata2;
+	postdelay_time = to_prime(post_delay,srate);
+	md_gain = gain;
+	lpfreq = lp_freq;
+	ap1 = ap2 = 0;
+	ap3 = 0;
+	lp1 = 0;
+	delay1 = delay2 = postdelay = 0;
+	sr = srate;
+	out1 = out2 =out3 = 0.0;
+	damping_  = false;
+}
+
+medium_diffuser::~medium_diffuser()
+{	
+		delete ap1;	
+		delete ap2;	
+		delete ap3;	
+		delete lp1;	
+		delete delay1;	
+		delete delay2;
+}
+
+
+bool medium_diffuser::create(void)
+{
+	if(ap1_data.time1 == 0 || ap1_data.gain1 <= 0.0){
+#ifdef _DEBUG
+		printf("\nmedium_diffuser: bad parameters(1)");
+#endif		
+		return false;
+	}
+	if(ap1_data.time2 ==0 || ap1_data.gain2 <= 0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad parameters(2)");
+#endif		
+		return false;
+	}
+	if(ap2_data.time1 == 0 || ap2_data.gain1 < 0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad parameters(3)");
+#endif		
+		return false;
+	}
+	if(ap2_data.time2 ==0 || ap2_data.gain2 < 0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad parameters(4)");
+#endif		
+		return false;
+	}
+	if(sr <=0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad srate parameter)");
+#endif		
+		return false;
+	}
+	if(lpfreq <0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad freq parameter)");
+#endif		
+		return false;
+	}
+
+	ap1 = new nested_allpass(sr,lpfreq,ap1_data.time1,ap1_data.time2,ap1_data.time3,
+							ap1_data.gain1,ap1_data.gain2,ap1_data.gain3);
+	if(!ap1->create()){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: can't create first diffuser");
+#endif		
+		return false;
+	}
+	ap2 = new nested_allpass(sr,lpfreq,ap2_data.time1,ap2_data.time2,ap2_data.time3,
+							ap2_data.gain1,ap2_data.gain2,ap2_data.gain3);
+	if(!ap2->create()){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: can't create second diffuser");
+#endif	
+		delete ap1; ap1 = 0;
+		return false;
+	}
+
+	ap3 = new allpassfilter((long)sr,to_prime(0.03,sr),0.5,to_prime(0.005,sr));
+	if(!ap3->create()){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: can't create internal allpass filter");
+#endif
+		delete ap1; ap1 = 0;
+		delete ap2; ap2 = 0;
+		return false;
+	}
+
+
+	if(postdelay_time > 0){
+		postdelay = new delay(postdelay_time,1.0);
+		if(!postdelay->create()){
+#ifdef _DEBUG
+			printf("\nmedium diffuser: can't create postdelay");
+#endif		
+
+			delete ap1; ap1 = 0;
+			delete ap2; ap2 = 0;
+			delete ap3; ap3 = 0;
+			return false;
+		}
+	}
+
+	//internal fixed delays
+	delay1 = new delay(to_prime(/*0.007*/0.067,sr),1.0);
+	if(!delay1->create()){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: can't create internal delay1");
+#endif
+		delete ap1; ap1 = 0;
+		delete ap2; ap2 = 0;
+		delete ap3; ap3 = 0;
+		delete [] postdelay; postdelay = 0;
+		return false;
+	}
+		delay2 = new delay(to_prime(0.015,sr),1.0);
+	if(!delay2->create()){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: can't create internal delay2");
+#endif
+		delete ap1; ap1 = 0;
+		delete ap2; ap2 = 0;
+		delete ap3; ap3 = 0;
+		delete [] postdelay; postdelay = 0;
+		delete delay1; delay1 = 0;
+		return false;
+	}
+	
+	lp1 = new onepole(lpfreq,sr,LOW_PASS);
+	if(lp1 == 0){
+		delete ap1; ap1 = 0;
+		delete ap2; ap2 = 0;
+		delete ap3; ap3 = 0;
+		delete [] postdelay; postdelay = 0;
+		delete delay1; delay1 = 0;
+		delete delay2; delay2 = 0;
+		return false;
+	}
+	if(lpfreq == 0.0)
+		damping_ = false;
+	ap1->set_damping(damping_);
+	ap2->set_damping(damping_);
+	return true;
+}
+
+double medium_diffuser::tick(double insam)
+{
+	double filter_out;
+	double output;
+
+	//feedback from postdelay, thru lp filter
+	// RWD NOTE: if the damping_ mechanism is preferred, this filter can be omitted and...
+	filter_out = md_gain * lp1->tick(postdelay->tick(out3));
+	// ... just use this
+	//filter_out = (float)(md_gain * postdelay->tick(out3));
+	
+	//first nested-allpass, takes input plus feedback
+	out1 = ap1->tick(insam + filter_out);
+	//inner allpass with predelay, followed by plain delay
+	out2 = delay1->tick(ap3->tick(out1));
+	//second nested-allpass, takes direct input
+	out3 = ap2->tick(insam + md_gain * delay2->tick(out2));
+	output = 0.5 * (out1 + out2 + out3);
+	return output;
+}
+
+
+large_diffuser::large_diffuser(const NESTDATA *apdata1,
+					const NESTDATA *apdata2,
+					double gain,
+					double lp_freq,
+					double srate)
+{
+	ap1_data = *apdata1;
+	ap2_data = *apdata2;
+	ld_gain = gain;
+	lpfreq = lp_freq;
+	ap1 = ap2 = 0;
+	ap3 = ap4 = 0;
+	lp1 = 0;
+	delay1 = delay2 = delay3 = delay4 = 0;
+	sr = srate;
+	out1 = out2 =out3 = 0.0;
+	damping_ = false;
+}
+
+large_diffuser::~large_diffuser()
+{	
+	delete ap1;	
+	delete ap2;
+	delete ap3;
+	delete ap4;
+ 	delete lp1;
+	delete delay1;
+	delete delay2;
+	delete delay3;
+	delete delay4;
+}
+
+bool large_diffuser::create(void)
+{
+	if(ap1_data.time1 == 0 || ap1_data.gain1 <= 0.0){
+#ifdef _DEBUG
+		printf("\nmedium_diffuser: bad parameters(1)");
+#endif		
+		return false;
+	}
+	if(ap1_data.time2 ==0 || ap1_data.gain2 <= 0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad parameters(2)");
+#endif		
+		return false;
+	}
+	if(ap2_data.time1 == 0 || ap2_data.gain1 < 0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad parameters(3)");
+#endif		
+		return false;
+	}
+	if(ap2_data.time2 ==0 || ap2_data.gain2 < 0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad parameters(4)");
+#endif		
+		return false;
+	}
+	if(sr <=0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad srate parameter)");
+#endif		
+		return false;
+	}
+	if(lpfreq <0.0){
+#ifdef _DEBUG
+		printf("\nmedium diffuser: bad freq parameter)");
+#endif		
+		return false;
+	}
+	ap1 = new nested_allpass(sr,lpfreq,ap1_data.time1,ap1_data.time2,ap1_data.time3,
+							ap1_data.gain1,ap1_data.gain2,ap1_data.gain3);
+	if(!ap1->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create first nested_allpass");
+#endif		
+		return false;
+	}
+	ap2 = new nested_allpass(sr,lpfreq,ap2_data.time1,ap2_data.time2,ap2_data.time3,
+							ap2_data.gain1,ap2_data.gain2,ap2_data.gain3);
+	if(!ap2->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create second netsed_allpass");
+#endif	
+		cleanup();
+		return false;
+	}
+
+	ap3 = new allpassfilter((long) sr,to_prime(0.008,sr),0.3);
+	if(!ap3->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create first internal allpass filter");
+#endif
+		cleanup();
+		return false;
+	}
+
+	ap4 = new allpassfilter((long)sr,to_prime(0.012,sr),0.3);
+	if(!ap4->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create second internal allpass filter");
+#endif
+		cleanup();
+		return false;
+	}
+	//internal fixed delays
+	delay1 = new delay((unsigned int)(0.004 * sr),1.0);
+	if(!delay1->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create internal delay1");
+#endif
+		cleanup();
+		return false;
+	}
+		delay2 = new delay((unsigned int)(0.017 * sr),1.0);
+	if(!delay2->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create internal delay2");
+#endif
+		cleanup();
+		return false;
+	}
+	delay3 = new delay((unsigned int)(0.031 * sr),1.0);
+	if(!delay3->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create internal delay3");
+#endif
+		cleanup();
+		return false;
+	}
+	delay4 = new delay((unsigned int)(0.003 * sr),1.0);
+	if(!delay4->create()){
+#ifdef _DEBUG
+		printf("\nlarge diffuser: can't create internal delay4");
+#endif
+		cleanup();
+		return false;
+	}
+	
+	lp1 = new onepole(lpfreq,sr,LOW_PASS);
+	if(lp1 == 0){
+		cleanup();
+		return false;
+	}
+	if(lpfreq == 0.0)
+		damping_ = false;
+	
+	ap1->set_damping(damping_);
+	ap2->set_damping(damping_);
+	return true;
+}
+void  large_diffuser::cleanup()
+{
+	delete lp1; lp1 = 0;
+	delete delay4; delay4 = 0;
+	delete delay3; delay3 = 0;
+	delete delay2; delay2 = 0;
+	delete delay1; delay1 = 0;
+	delete ap4; ap4 = 0;
+	delete ap3; ap3 = 0;
+	delete ap2; ap2 = 0;
+	delete ap1; ap1 = 0;
+}
+
+// NOTE: the three output taps from the diffuser can lead to discrete slap-back echoes in the output.
+// possibly authentic, but may not suit all purposes. See rmverb.cpp: some delay lengths reduced from Gardner.
+// This is possible because delays are discrete: they do not share single delay line buffer.
+// true also of medium diffuser, but not so apparent as delays are smaller
+double large_diffuser::tick(double insam)
+{
+	double filter_out;
+	double output;
+	//feedback from lpfilter
+	// RWD NOTE: if the damping_ mechanism is preferred, this filter can be omitted and ...
+	filter_out = ld_gain * lp1->tick(out3);		
+	// ... just use this
+	//filter_out = (float)(ld_gain * out3);			  
+
+	//first tap, from two plain allpasses, and trailing delay
+	out1 = delay1->tick(ap4->tick(ap3->tick(insam + filter_out)));
+	//second tap, from first nested allpass and delays at entrance and exit
+	out2 = delay3->tick(ap1->tick(delay2->tick(out1)));
+	//third tap, from second nested-allpass with leading delay
+	out3 = ap2->tick(delay4->tick(out2));
+	//out3 also goes into lpfilter... see above
+	//prescribed scale factors from Gardner
+	output = (out1 * 0.34) + (out2 * 0.14) + (out3 * 0.14);
+	// this seems to give relatively soft reverb tail:  could try this:
+	//output = (float)((out1 * 0.7) + (out2 * 0.25) + (out3 * 0.25));
+	return  output;
+}
+
+
+moorer::moorer(const MOORERDATA *pdata,double reverbtime,double damping,double srate)
+{
+	mdata		= *pdata;
+	reverb_time = reverbtime;
+	dampfac		= damping;
+	sr			= srate;
+	comb1 = comb2 = comb3 = comb4 = comb5 = comb6 = 0;
+	ap1 = ap2 = ap3 =  ap4  = 0;
+	out = 0.0;
+}
+
+moorer::~moorer()
+{
+	delete comb1;
+	delete comb2;
+	delete comb3;
+	delete comb4;
+	delete comb5;
+	delete comb6;
+	delete ap1;
+	delete ap2;		
+	delete ap3;
+	delete ap4;
+}
+
+bool moorer::create(void)
+{
+	fgain1 = 0.42 * dampfac;
+	fgain2 = 0.46 * dampfac;
+	fgain3 = 0.48 * dampfac;
+	fgain4 = 0.49 * dampfac;
+	fgain5 = 0.50 * dampfac;
+	fgain6 = 0.52 * dampfac;
+	scalefac = (1.0 / 6.0);	
+	cgain1 = exp(log(0.001)*(mdata.ctime1/reverb_time)) * (1. - fgain1);
+	cgain2 = exp(log(0.001)*(mdata.ctime2/reverb_time)) * (1. - fgain2);
+	cgain3 = exp(log(0.001)*(mdata.ctime3/reverb_time)) * (1. - fgain3);
+	cgain4 = exp(log(0.001)*(mdata.ctime4/reverb_time)) * (1. - fgain4);
+	cgain5 = exp(log(0.001)*(mdata.ctime5/reverb_time)) * (1. - fgain5);
+	cgain6 = exp(log(0.001)*(mdata.ctime6/reverb_time)) * (1. - fgain6);
+
+	comb1 = new lpcomb(to_prime(mdata.ctime1,sr),cgain1,fgain1,1.0);
+	comb2 = new lpcomb(to_prime(mdata.ctime2,sr),cgain2,fgain2,1.0);
+	comb3 = new lpcomb(to_prime(mdata.ctime3,sr),cgain3,fgain3,1.0);
+	comb4 = new lpcomb(to_prime(mdata.ctime4,sr),cgain4,fgain4,1.0);
+	comb5 = new lpcomb(to_prime(mdata.ctime5,sr),cgain5,fgain5,1.0);
+	comb6 = new lpcomb(to_prime(mdata.ctime6,sr),cgain6,fgain6,1.0);
+	ap1		= new allpassfilter((long)sr,to_prime(mdata.atime1,sr),-0.6,0);
+	ap2		= new allpassfilter((long)sr,to_prime(mdata.atime2,sr),0.4,0);
+	ap3		= new allpassfilter((long)sr,to_prime(mdata.atime3,sr),-0.61,0);
+	ap4		= new allpassfilter((long)sr,to_prime(mdata.atime4,sr),0.39,0);
+
+	if(!
+		(comb1->create()
+		&& comb2->create()
+		&& comb3->create()
+		&& comb4->create()
+		&& comb5->create()
+		&& comb6->create()
+		&& ap1->create()
+		&& ap2->create()
+		&& ap3->create()
+		&& ap4->create()
+		)) {
+		cleanup();
+		return false;
+	}
+	return true;
+}
+
+double moorer::tick(double insam)
+{
+	out =  comb1->tick(insam) +
+		  comb2->tick(insam) +
+		  comb3->tick(insam) +
+		  comb4->tick(insam) +
+		  comb5->tick(insam) +
+		  comb6->tick(insam)
+		  ;
+	out *= scalefac;
+	out = ap4->tick(ap3->tick(ap2->tick(ap1->tick(out))));
+	return out;
+}
+
+void moorer::cleanup()
+{
+	delete comb6; comb6 = 0;
+	delete comb5; comb5 = 0;
+	delete comb4; comb4 = 0;
+	delete comb3; comb3 = 0;
+	delete comb2; comb2 = 0;
+	delete comb1; comb1 = 0;
+	delete ap4; ap4 = 0;
+	delete ap3; ap3 = 0;
+	delete ap2; ap2 = 0;
+	delete ap1; ap1 = 0;
+}
+
+
+/* vmtdelay *******************************/
+vcomb4::vcomb4()
+{
+	dl_buf = 0;
+	dl_length = 0;
+	dl_input_index = 0;
+	dl_srate = 0.0;
+	
+}
+
+vcomb4::~vcomb4()
+{
+	if(dl_buf && dl_length > 0)
+		delete [] dl_buf;
+		
+}
+
+bool vcomb4::init(long srate,double length_secs)
+{
+	unsigned long len_frames = (unsigned long)(srate * length_secs );	
+	if(len_frames == 0)
+		return false;
+	/*  round upwards to allow for interpolation */
+	len_frames++;
+	/* reject init  if already created*/
+	/* TODO: more sexy error checking/reporting... */
+	if(dl_buf) 
+		return false;
+	try {
+		dl_buf = new double[len_frames];
+	}
+	catch(...){
+		return false;
+	}
+	/* for VC 6 */
+	if(dl_buf == 0)
+		return false;
+	for(unsigned long i = 0; i < len_frames;i++)
+		dl_buf[i] = 0.0;
+
+	dl_length = len_frames;
+	dl_input_index = 0;
+	dl_srate = srate;
+	gain = 0.5;
+	return true;
+}
+
+bool vcomb4::init(const MOORERDATA *p_mdata,double reverbtime,double damping,double sr)
+{
+	if(!init((long) sr,reverbtime))
+		return false;
+    return true;
+}
+
+double vcomb4::tick(double vdelay_1,double vdelay_2,double vdelay_3,double vdelay_4,double feedback,double input)
+{
+	unsigned long base_readpos1,base_readpos2,base_readpos3,base_readpos4; 
+	unsigned long next_index1,next_index2,next_index3,next_index4; 
+	double vlength,dlength,frac;
+	double readpos1,readpos2,readpos3,readpos4;
+	double* buf = dl_buf;
+	
+	dlength = (double) dl_length;					/* get maxlen, save a cast later on */
+	/* read pointer is vlength ~behind~ write pointer */
+	/* tap1*/
+	vlength = dlength -  (vdelay_1  * dl_srate);
+	vlength = vlength < dlength ? vlength : dlength;	  /* clip vdelay to max del length */
+	readpos1 = dl_input_index + vlength;
+	base_readpos1 = (unsigned long) ( readpos1);				  
+	if(base_readpos1 >= dl_length)					  /* wrap dl indices */
+		base_readpos1 -= dl_length;
+	next_index1 = base_readpos1 + 1;
+	if(next_index1 >= dl_length)
+		next_index1 -= dl_length;
+	/*tap2*/
+	vlength = dlength -  (vdelay_2  * dl_srate);
+	vlength = vlength < dlength ? vlength : dlength;
+	readpos2 = dl_input_index + vlength;
+	base_readpos2 = (unsigned long) ( readpos2);				  
+	if(base_readpos2 >= dl_length)					  /* wrap dl indices */
+		base_readpos2 -= dl_length;
+	next_index2 = base_readpos2 + 1;
+	if(next_index2 >= dl_length)
+		next_index2 -= dl_length;
+	/*tap3*/
+	vlength = dlength -  (vdelay_3  * dl_srate);
+	vlength = vlength < dlength ? vlength : dlength;
+	readpos3 = dl_input_index + vlength;
+	base_readpos3 = (unsigned long) ( readpos3);				  
+	if(base_readpos3 >= dl_length)					  /* wrap dl indices */
+		base_readpos3 -= dl_length;
+	next_index3 = base_readpos3 + 1;
+	if(next_index3 >= dl_length)
+		next_index3 -= dl_length;
+	/* tap4 */
+	vlength = dlength -  (vdelay_4  * dl_srate);
+	vlength = vlength < dlength ? vlength : dlength;
+	readpos4 = dl_input_index + vlength;
+	base_readpos4 = (unsigned long) ( readpos4);				  
+	if(base_readpos4 >= dl_length)					  /* wrap dl indices */
+		base_readpos4 -= dl_length;
+	next_index4 = base_readpos4 + 1;
+	if(next_index4 >= dl_length)
+		next_index4 -= dl_length;
+
+
+	double outsum = 0.0;
+	fb1 *= feedback;
+	fb2 *= feedback;
+	fb3 *= feedback;
+	fb4 *= feedback;
+	/* basic interp of variable delay pos */
+	frac = readpos1 - (int) readpos1;	
+	outsum += fb1 * (buf[base_readpos1]+((buf[next_index1] - buf[base_readpos1]) * frac));
+	frac = readpos2 - (int) readpos2;	
+	outsum += fb2 * (buf[base_readpos2]+((buf[next_index2] - buf[base_readpos2]) * frac));
+	frac = readpos3 - (int) readpos3;	
+	outsum += fb3 * (buf[base_readpos3]+((buf[next_index3] - buf[base_readpos3]) * frac));
+	frac = readpos4 - (int) readpos4;	
+	outsum += fb4 * (buf[base_readpos4]+((buf[next_index4] - buf[base_readpos4]) * frac));
+	/* how do we scale all this?*/
+	input += gain * (outsum);	
+	/* add in new sample + fraction of ouput, unscaled, for minimum decay at max feedback */
+	buf[dl_input_index++] = input   + outsum * 0.25;
+	if(dl_input_index == dl_length)
+		dl_input_index = 0;
+	return  input + outsum * 0.25;
+
+}

dev/externals/reverb/reverberator.h

@@ -0,0 +1,360 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+
+// reverberator.h: interface for the reverberator class.
+// RWD.9.98 now contains all my delay/allpass etc  ugs
+//RWD 03.20  corrected def of tick function to double everywhere
+//////////////////////////////////////////////////////////////////////
+
+#if !defined(AFX_REVERBERATOR_H__9DE45E43_83BA_11D1_96D4_444553540000__INCLUDED_)
+#define AFX_REVERBERATOR_H__9DE45E43_83BA_11D1_96D4_444553540000__INCLUDED_
+
+#if _MSC_VER >= 1000
+#pragma once
+#endif // _MSC_VER >= 1000
+
+#include "wavetable.h"
+
+enum {LOW_PASS,HIGH_PASS};
+enum {LOW_SHELVE,HIGH_SHELVE};
+
+unsigned int to_prime(double dur,double sr);
+
+
+
+typedef struct __tap
+{
+	double pos;
+	double val;
+}
+deltap;
+
+typedef struct nested_allpass_data{
+	unsigned int time1;
+	unsigned int time2;
+	unsigned int time3;
+	unsigned int delay1;
+	unsigned int delay2;
+	double gain1;
+	double gain2;
+	double gain3;
+}
+NESTDATA;
+
+typedef struct moorer_data{
+	double ctime1;
+	double ctime2;
+	double ctime3;
+	double ctime4;
+	double ctime5;
+	double ctime6;
+	double atime1;
+	double atime2;
+	double atime3;
+	double atime4;
+}
+MOORERDATA;
+//abstract base class for all processors
+
+class cdp_processor 
+{
+public:
+	cdp_processor() {}
+	virtual ~cdp_processor(){}
+	//not all processes will need to use this
+	virtual bool create() {return true;}
+//	virtual float tick(float insam)	{return insam;}
+    virtual double tick(double insam) = 0;
+};
+
+
+//could have a method create_taps(uint ntaps, float maxtime)  ... ?
+class tapdelay 
+{
+public:
+	tapdelay(const deltap *taps,unsigned int ctaps,double sr);
+	bool	create(void);
+	double  get_maxgain(void);
+	void set_gain(double gain);
+	double	tick(double insam);
+	virtual ~tapdelay();
+private:
+	double *tapbuf;
+	double srate;
+	double one_over_ntaps;
+	double max_gain;		/* get rms value of whole array */
+	double output_gain;
+	unsigned int ntaps;
+	unsigned int buflen;
+	unsigned int *tapptr;
+	double *tapcoeff;
+	unsigned int reader;
+	const deltap *dtaps;
+};
+
+//basic feedthrough delay (independent of srate) - could be base class?
+class delay
+{
+public:
+	delay(unsigned int len, double gain);
+	virtual ~delay();
+	bool create(void);
+	double tick(double insam);
+private:
+	double *delbuf;
+	unsigned int buflen;
+	double dgain;
+	unsigned int ptr;
+};
+
+
+
+class allpassfilter  
+{
+public:
+	allpassfilter(long sr,unsigned int time, double ingain,unsigned int predelay);
+	virtual ~allpassfilter();
+	bool create(void);    
+    void set_vfreq(double freq);
+    void set_vmod(double amount);
+	double tick(double insam);
+    double vtick(double insam);
+
+private:
+	double gain;
+	double *rvbuf1;			 //CDP allpass uses 2 buffers!
+    double vfreq,vmod;
+    fastlfo *sinelfo;
+    fastlfo *noiselfo;
+	unsigned int rvblen, writer1, reader1;
+	delay *pre_dline;
+	unsigned int pre_delay;
+    long srate;
+};
+
+///////////////// DOUBLE-NESTED ALLPASS
+///  set either or both gain2, gain3 to zero to disable the respective allpass
+////	
+class onepole;
+
+class  nested_allpass
+{
+public:
+	nested_allpass(double srate,double lpfreq,unsigned int outertime, 
+					unsigned int time1,
+					unsigned int time2,
+					double gain, 
+					double gain1,
+					double gain2, 
+					unsigned int delay1, 
+					unsigned int delay2);
+	virtual ~nested_allpass();
+	bool create(void);
+	double tick(double insam);
+    void  set_vfreq(double freq);
+    void set_vmod(double amount);
+    double vtick(double  insam);
+	void set_damping(bool damping) { damping_ = damping;}
+	bool get_damping() const { return damping_;}
+private:
+	
+	double srate_, outer_gain, ap1_gain,ap2_gain;
+	unsigned int ap1_length,ap2_length;
+	unsigned int outer_time,ap1_delay,ap2_delay;
+	allpassfilter *ap1, *ap2;
+	double *buf;
+	unsigned int writer,reader;
+	// RWD experimental - internal lf damping
+	onepole* lp; 
+	bool damping_;
+	double sr_,lpfreq_;
+
+};
+
+
+//TODO: comb filter (with lopass), and tapped delay
+
+//basic lowpass for combs
+
+class lowpass1
+{
+public:
+	lowpass1(double filtgain);
+	virtual ~lowpass1();
+	double tick(double insam);
+private:
+	double temp,gain,output;	
+};
+
+
+class onepole
+{
+public:
+	onepole(double freq,double srate,int low_high);
+	virtual ~onepole();
+	double tick(double input);
+private:
+	double a1,a2,b,lastout;
+};
+
+class tonecontrol
+{
+public:
+	tonecontrol(double frq,double dbfac,int type,double srate);
+	virtual ~tonecontrol();
+	bool create(void);
+	double tick(double input);
+private:
+	double freq, boost, sr,a0, a1, a2, b1, b2;
+	double x1,x2,y1,y2;
+	int tc_type;
+	void lshelve(void);
+	void hshelve(void);
+};
+
+class lpcomb
+{
+public:
+	lpcomb(unsigned int time, double ingain, double lpcoeff, double prescale);
+	virtual ~lpcomb();
+	bool create(void);
+	double tick(double insam);
+private:
+	double prescale, gain, lpgain;
+	double *combbuf;
+	unsigned int rvblen, reader1,writer1;
+	lowpass1 *lp;
+
+};
+
+
+class small_diffuser : public cdp_processor
+{
+public:
+	small_diffuser(unsigned int pre_delay,
+					const NESTDATA *apdata1,
+					const NESTDATA *apdata2,double fb_gain,double lp_freq,double srate);
+	virtual ~small_diffuser();
+	bool create(void);
+	double tick(double insam);
+	void set_damping(bool damping) { damping_ = damping; }
+private:
+	NESTDATA ap1_data,ap2_data;
+	nested_allpass *ap1,*ap2;
+	delay *predelay;
+	//lowpass1 *lp1;
+	onepole *lp1;
+	unsigned int predelay_time;
+	double lpgain,lpfreq,sr;
+	double out1,out2;
+	bool damping_;
+};
+
+//no public output from postdelay, so internal only
+class medium_diffuser : public cdp_processor
+{
+public:
+	medium_diffuser(double post_delay,const NESTDATA *apdata1,
+					const NESTDATA *apdata2,
+					double gain,
+					double lp_freq,
+					double srate);
+	virtual ~medium_diffuser();
+	bool create(void);
+	double tick(double insam);
+	void set_damping(bool damping) { damping_ = damping;	}
+private:
+	NESTDATA ap1_data,ap2_data;
+	nested_allpass *ap1,*ap2;
+	allpassfilter *ap3;				//no public access to this, for now...
+	onepole *lp1;
+	delay *delay1,*delay2,*postdelay;			
+	double md_gain,lpfreq,sr;
+	double out1,out2,out3;
+	unsigned int postdelay_time;
+	bool damping_;
+};
+
+
+class large_diffuser : public cdp_processor
+{
+public:
+	large_diffuser(const NESTDATA *apdata1,
+					const NESTDATA *apdata2,
+					double gain,
+					double lp_freq,
+					double srate);
+	virtual ~large_diffuser();
+	bool create(void);
+	double tick(double insam);
+	void set_damping(bool damping) { damping_ = damping; }
+private:
+	void cleanup();
+	NESTDATA ap1_data,ap2_data;
+	nested_allpass *ap1,*ap2;
+	allpassfilter *ap3,*ap4;				//no public access to this, for now...
+	onepole *lp1;
+	delay *delay1,*delay2,*delay3,*delay4;			
+	double ld_gain,lpfreq,sr;
+	double out1,out2,out3;
+	bool damping_;
+};
+
+class moorer : public cdp_processor
+{
+public:
+	moorer(const MOORERDATA *p_mdata,double reverbtime,double damping,double sr);
+	virtual ~moorer();
+	bool create(void);
+	double tick(double insam);
+private:
+	void cleanup();
+	lpcomb *comb1,*comb2,*comb3,*comb4,*comb5,*comb6;
+	allpassfilter *ap1,*ap2,*ap3,*ap4;
+	MOORERDATA mdata;
+	double cgain1,cgain2,cgain3,cgain4,cgain5,cgain6;
+	double fgain1,fgain2,fgain3,fgain4,fgain5,fgain6;	
+	double reverb_time,dampfac,sr;
+	double out,scalefac;
+};
+
+/* vmtdelay ****************/
+class vcomb4 
+{
+public:
+	vcomb4();
+	virtual ~vcomb4();
+	bool init(long srate,double length_secs);
+	bool init(const MOORERDATA *p_mdata,double reverbtime,double damping,double sr);
+	double tick(double vdelay_1,double vdelay_2,double vdelay_3,double vdelay_4,double feedback,double input);
+	void setfgains(double fg1,double fg2,double fg3,double fg4) {fb1 = fg1;fb2 = fg2;fb3 = fg3;fb4 = fg4;}
+private:
+	double*			dl_buf;
+	unsigned long	dl_length;
+	unsigned  long	dl_input_index;
+	double			dl_srate;
+	double			fb1,fb2,fb3,fb4;
+//	lowpass1		lp1,lp2,lp3,lp4;
+	double		gain;
+};
+
+
+#endif // !defined(AFX_REVERBERATOR_H__9DE45E43_83BA_11D1_96D4_444553540000__INCLUDED_)

dev/externals/reverb/rmverb.cpp

@@ -0,0 +1,949 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+#include <algorithm>
+#include <vector>
+
+using namespace std;
+// TODO: clean up "new" handling, move to VC7 and get to use exceptions!
+enum cmdargs {
+		PROGNAME,
+		INFILE,
+		OUTFILE,
+		REVERBTYPE,
+		RGAIN,
+		MIX,
+		FEEDBACK,		
+		LPABSORB,
+		LPFREQ,
+		TRTIME,
+		CONFIG
+};
+
+enum reverb_type {
+		DIFF_SMALL,
+		DIFF_MEDIUM,
+		DIFF_LARGE
+};
+
+extern "C" {
+#include "portsf.h"
+	   //NB requires stdio.h etc - time to change this?
+}
+
+#include "reverberator.h"
+//RWD.10.98 TODO: for final versions, MUST limit input delay times! can't have 108sec delays.....
+//get all reflection data
+#include "reflect.cpp"
+
+
+long readtaps(FILE *fp, deltap **taps,double sr);		
+
+const char* cdp_version = "6.0.0";
+
+void usage(void);
+
+int
+main(int argc,char *argv[])
+
+{
+	int i,ifd, ofd,reverbtype = DIFF_SMALL;
+	FILE *fp = 0;
+	long	chans,outchans = 2;
+	double	sr,trailertime;
+	cdp_processor *cdpdiff = 0;
+	NESTDATA data1,data2;
+	PSF_CHPEAK *peaks;
+	PSF_PROPS props;
+	unsigned int time1a,time1b,time1c,time2a,time2b,time2c;
+	double feedback,fb_lpfreq, lp_freq,predelay = -1.0;
+	double nyquist,pseudo_nyquist;
+	double dry_gain,wet_gain;
+	//double early_gain = 1.0;
+	double diff_gain = 1.0;
+	// nested-allpass diffusers
+	small_diffuser *s_diffuser = 0;
+	medium_diffuser *m_diffuser = 0;
+	large_diffuser *l_diffuser = 0;
+	// traioling plain allpasses for front pair de-correlation
+	allpassfilter *ap1_l = 0,*ap2_l = 0;		
+	allpassfilter *ap1_r = 0,*ap2_r = 0;
+	// decorrelation for other channels (surround)
+	allpassfilter **chan_ap = 0;
+	onepole *lp = 0;
+	// tonecontrols max -12dB slope, variable; 
+	tonecontrol *tc_l_lowcut = 0, *tc_r_lowcut = 0;	
+	tonecontrol *l_highcut = 0, *r_highcut = 0;
+	double lowcut_freq = 0.0,highcut_freq = 0.0;
+	tapdelay *tdelay = 0;	
+	long ntaps = 0;
+	deltap *ptaps = 0;
+	bool want_mono = false;
+	bool want_floats = false;
+	bool usertaps = false;
+	bool double_damping = false;
+
+	if((argc==2) && strcmp(argv[1],"--version")==0) {
+		fprintf(stdout,"%s\n",cdp_version);
+		fflush(stdout);
+		return 0;
+	}
+	if(psf_init()){
+		fprintf(stderr,"rmverb: initialisation failed\n");
+		return 1;
+	}
+    	
+	if(argc < CONFIG){
+		if(argc > 1)
+			fprintf(stderr,"rmverb: insufficient arguments\n");
+		usage();
+	}
+	
+	while((argc > 1) && argv[1][0]=='-'){
+		char *arg = argv[1];
+		switch(*(++arg)){
+		case('p'):
+			if(*(++arg) =='\0'){
+				fprintf(stderr,"rmverb: predelay requires a parameter\n");
+				usage();
+			}
+			
+			predelay = atof(arg);
+			if(predelay <0.0){
+				fprintf(stderr,"rmverb: predelay must be >= 0.0\n");
+				usage();
+			}
+			predelay *= 0.001;			//get msecs from user
+			break;
+		case('c'):
+			if(*(++arg) == '\0') {
+				fprintf(stderr,"rmverb: -c flag requires a value\n");
+				usage();
+			}
+			outchans = atoi(arg);
+			if(outchans < 1 || outchans > 16){
+				fprintf(stderr,"rmverb: impossible channel value requested\n");
+				usage();
+			}
+			if(outchans==1)
+				want_mono = true;			
+			break;
+		case('d'):
+			double_damping = true;
+			break;
+		case('e'):
+			if(*(++arg)=='\0'){
+				fprintf(stderr,"rmverb: -e flag needs a filename\n");
+				usage();
+			}
+			if((fp = fopen(arg,"r")) ==0){
+				fprintf(stderr,"rmverb: unable to open breakpoint file %s\n",arg);
+				exit(1);
+			}
+			usertaps = true;
+			break;
+			// -f is CDP-specific: read old 32bit int format as floats
+		case('f'):
+			want_floats = true;
+			if(*(++arg) != '\0')
+				fprintf(stderr,"rmverb WARNING: -f flag does not take a parameter\n");
+			break;
+		case('L'):
+			if(*(++arg) == '\0'){
+				fprintf(stderr,"rmverb: -L flag needs frequency argument\n");
+				usage();
+			}			
+			lowcut_freq = atof(arg);
+			if(lowcut_freq <= 0.0){
+				fprintf(stderr,"rmverb: Lowcut freq must be greater than zero\n");
+				usage();
+			}
+			break;
+		case('H'):
+			if(*(++arg) == '\0'){
+				fprintf(stderr,"rmverb: -H flag needs frequency argument\n");
+				usage();
+			}			
+			highcut_freq = atof(arg);
+			if(highcut_freq <= 0.0){
+				fprintf(stderr,"rmverb: Highcut freq must be greater than zero\n");
+				usage();
+			}
+			break;
+		default:
+			fprintf(stderr,"rmverb: illegal flag option %s\n",argv[1]);
+			usage();
+			break;
+		}
+		argc--;
+		argv++;
+	}
+
+	if(argc < CONFIG){
+		printf("rmverb: insufficient arguments\n");
+		usage();
+	}
+#ifdef _DEBUG
+	if(predelay > 0.0)
+		printf("got predelay = %.4lf\n",predelay);
+#endif
+
+	chan_ap = new allpassfilter*[outchans];
+	if(chan_ap==0){
+		fprintf(stderr,"rmverb: no memory for multi-channel diffusion\n");
+		exit(1);
+	}
+
+	if((ifd = psf_sndOpen(argv[INFILE],&props,0)) < 0) {
+				fprintf(stderr,"\nrmverb: cannot open input file %s", argv[INFILE]);
+				exit(1);
+	}
+	
+	if(props.format  <= PSF_FMT_UNKNOWN || props.format > PSF_AIFC){
+		fprintf(stderr,"infile is not a recognised soundfile\n");
+		psf_sndClose(ifd);
+		return 1;
+	}
+	
+	sr = (double) props.srate;
+	nyquist = sr / 2.0;
+	pseudo_nyquist = nyquist * 0.7; // filters not reliable close to nyquist
+
+	chans = props.chans;
+	if(chans > 2){
+		fprintf(stderr,"rmverb can only accept mono or stereo files\n");
+		exit(1);
+	}
+		
+	
+	reverbtype =  atoi(argv[REVERBTYPE]) - 1;
+	if((reverbtype < DIFF_SMALL) || (reverbtype > DIFF_LARGE)){
+		fprintf(stderr,"rmverb: rmsize must be 1,2 or 3\n");
+		exit(1);
+	}
+	
+	diff_gain = atof(argv[RGAIN]);
+	if(diff_gain <= 0.0 || diff_gain > 1.0){
+		fprintf(stderr,"rgain must be > 0 and <= 1.0\n");
+		exit(1);
+	}
+
+	dry_gain =  atof(argv[MIX]);			  //global output gain from diffuser
+	if(dry_gain < 0.0 || dry_gain > 1.0 ){
+		printf("reverb: mix must be  between 0.0 and 1.0\n");
+		usage();
+	}
+	wet_gain = 1.0f - dry_gain;
+	// some odd people like to have 100% wet signal...
+	wet_gain *= 2.0;					//not v scientific, but works...
+	
+	feedback = atof(argv[FEEDBACK]);			  //feedback in diffuser
+	if(feedback < 0.0 || feedback >1.0){
+		printf("rmverb: feedback must be within 0.0 - 1.0\n");
+		usage();
+	}
+	
+	fb_lpfreq = atof(argv[LPABSORB]);			 //feedback lp-filter in diffuser
+	if(fb_lpfreq < 0.0 || fb_lpfreq > pseudo_nyquist){
+		printf("rmverb: absorb must be within 0 to %f\n",pseudo_nyquist);
+		usage();
+	}
+
+	lp_freq = atof(argv[LPFREQ]);
+	if(lp_freq < 0.0 || lp_freq > pseudo_nyquist){
+		printf("rmverb: lpfreq must be within 0.0 to %.4lfHz\n",pseudo_nyquist);
+		usage();
+	}
+	trailertime = atof(argv[TRTIME]);
+	if(trailertime < 0.0){
+		fprintf(stderr,"rmverb: trtime must be >= 0.0\n");
+		usage();
+	}
+	//from -m flag; will override a stereo input too
+	if(want_mono)
+		outchans =  1;
+	// NB: now specifying two trailing allpasses for diffuser
+	// handles front pair; currently using chan_ap[] only for remaining channels
+	// values are experimental!
+	double apdelay = 0.02;
+	double apfac = 0.02 / outchans;
+	for(i = 0; i < outchans; i++){						// or use a random modifier?
+		chan_ap[i] =  new allpassfilter(sr,to_prime(apdelay - apfac * (double)i,sr),0.4,0);
+		if(chan_ap[i] ==0){
+			fprintf(stderr,"rmverb: no memory for output diffusers\n");
+			exit(1);
+		}
+		if(!chan_ap[i]->create()){
+			fprintf(stderr,"rmverb: no memory for output diffusers\n");
+			exit(1);
+		}
+	}
+    // custom early reflections from time/value text-file
+	if(usertaps){
+#ifdef _DEBUG
+		assert(fp);
+#endif
+		ntaps = readtaps(fp,&ptaps,sr);
+		if(ntaps==0){
+			fprintf(stderr,"rmverb: error reading tapfile\n");
+			exit(1);
+		}
+		printf("rmverb: loaded %ld early reflections from tapfile\n",ntaps);
+		fclose(fp);
+		fp = 0;
+	}
+
+
+	//create input low/high filters, if wanted
+	if(lowcut_freq > 0.0){
+		if(highcut_freq > 0.0 && highcut_freq <= lowcut_freq){
+			fprintf(stderr,"rmverb: Highcut frequency	 must be higher than Lowcut frequency\n");
+			usage();
+		}
+		//lowcut is based on low shelving eq filter; here fixed at 12dB
+		// but can easily add this as user param
+		tc_l_lowcut = new tonecontrol(lowcut_freq,-12.0,LOW_SHELVE,sr);
+		if(!tc_l_lowcut->create()){
+			fprintf(stderr,"rmverb: unable to create Lowcut filter\n");
+			exit(1);
+		}
+		if(chans==2){
+			tc_r_lowcut = new tonecontrol(lowcut_freq,-12.0,LOW_SHELVE,sr);
+			if(!tc_r_lowcut->create()){
+				fprintf(stderr,"rmverb: unable to create Lowcut filter\n");
+				exit(1);
+			}
+		}
+	}
+
+	if(highcut_freq > 0.0){		
+		l_highcut = new tonecontrol(highcut_freq,-12.0,HIGH_SHELVE,sr);
+		if(!l_highcut->create()){
+			fprintf(stderr,"rmverb: unable to create Highcut filter\n");
+				exit(1);
+		}
+		if(chans==2){
+			r_highcut = new tonecontrol(highcut_freq,-12.0,HIGH_SHELVE,sr);
+			if(!r_highcut->create()){
+				fprintf(stderr,"\nrmverb: unable to create Highcut filter");
+				exit(1);
+			}
+		}
+	}
+
+	//now create diffuser
+
+	switch(reverbtype){
+		case(DIFF_SMALL):
+			time1a = to_prime(0.0047, sr);
+			time1b = to_prime(0.022, sr);
+			time1c = to_prime(0.0083,sr);
+
+			time2a = to_prime(0.036,sr);
+			time2b = to_prime(0.03,sr);
+			time2c = 0;//dummy
+			data1.time1 = time1a;
+			data1.time2 = time1b;
+			data1.time3 = time1c;
+			data1.gain1 = 0.3;
+			data1.gain2 = 0.4;
+			data1.gain3 = 0.6;
+			data1.delay1 = 0;
+			data1.delay2 = 0;
+			data2.time1 = time2a;
+			data2.time2 = time2b;
+			data2.time3 = time2c;
+			data2.gain1 = 0.1;
+			data2.gain2 = 0.4;
+			data2.gain3 = 0.0;			//not used by gardner
+			data2.delay1 = 0;
+			data2.delay2 = 0;
+
+
+			break;
+		case(DIFF_MEDIUM):
+			time1a = to_prime(/*0.0047*/0.035, sr);
+			time1b = to_prime(0.0083, sr);
+			time1c = to_prime(0.022,sr);
+
+			time2a = to_prime(/*0.0292*/0.039,sr);
+			time2b = to_prime(0.0098,sr);
+			time2c = 0;//dummy
+			data1.time1 = time1a;
+			data1.time2 = time1b;
+			data1.time3 = time1c;
+			data1.gain1 = 0.3;
+			data1.gain2 = 0.7;
+			data1.gain3 = 0.5;
+			data1.delay1 = 0;
+			data1.delay2 = 0;
+			data2.time1 = time2a;
+			data2.time2 = time2b;
+			data2.time3 = time2c;
+			data2.gain1 = 0.3;
+			data2.gain2 = 0.6;
+			data2.gain3 = 0.0;			//not used by gardner
+			data2.delay1 = 0;
+			data2.delay2 = 0;
+
+
+
+			break;
+		case(DIFF_LARGE):
+			// VERY LARGE!
+			// NB: orig outer delays from Gardner generate very strong discrete echoes,
+			// and hence also some pulsing in the dense reverb
+			// for really large spaces (caverns) this is OK, otherwise we ant to control this
+			// via predelay, and earlies
+			// NOTE: this will NOT WORK when sharing a single delay line!
+			// in that case, outer delay MUST be > internal ones
+			// TODO: devise user param to control a range of delay values
+			// (maybe derived from feedback level?), to further scale room size
+			time1a = to_prime(0.03 /*0.087*/, sr);
+			time1b = to_prime(0.062, sr);
+			time1c = to_prime(0.022,sr);	//dummy: not used
+
+			time2a = to_prime(0.018/*0.12*/,sr);
+			time2b = to_prime(0.076,sr);
+			time2c = to_prime(0.03,sr);		
+			data1.time1 = time1a;
+			data1.time2 = time1b;
+			data1.time3 = time1c;
+			data1.gain1 = 0.5;
+			data1.gain2 = 0.25;
+			data1.gain3 = 0.0;			//not used
+			data1.delay1 = 0;
+			data1.delay2 = 0;
+			data2.time1 = time2a;
+			data2.time2 = time2b;
+			data2.time3 = time2c;
+			data2.gain1 = 0.5;
+			data2.gain2 = 0.25;
+			data2.gain3 = 0.25;		
+			data2.delay1 = 0;
+			data2.delay2 = 0;
+			break;
+		default:
+#ifdef _DEBUG
+			assert(false);
+#endif
+			break;
+	}
+#ifdef _DEBUG
+	//RWD: development only!
+	if(argc==CONFIG+1){
+		int got = 0;
+		double time1a,time1b,time1c,time2a,time2b,time2c,gain1a,gain1b,gain1c,gain2a,gain2b,gain2c;
+		FILE *fp = fopen(argv[CONFIG],"r");
+		if(!fp)
+			printf("rmverb: can't open datafile %s, using presets\n",argv[CONFIG]);
+		else {
+		   // time1a,gain1a,time2a,gain2a,time3a,gain3a...
+			printf("loading diffusion data...\n");
+			got = fscanf(fp,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
+			  &time1a,&gain1a,
+			  &time1b,&gain1b,
+			  &time1c,&gain1c,
+			  &time2a,&gain2a,
+			  &time2b,&gain2b,
+			  &time2c,&gain2c);
+			fclose(fp);
+			if(got != 12)
+				printf("rmverb: error reading data values\n");
+			else {
+				data1.time1 = to_prime(time1a,sr);
+				data1.time2 = to_prime(time1b,sr);
+				data1.time3 = to_prime(time1c,sr);
+				data1.gain1 = gain1a;
+				data1.gain2 = gain1b;
+				data1.gain3 = gain1c;
+
+				data2.time1 = to_prime(time2a,sr);
+				data2.time2 = to_prime(time2b,sr);
+				data2.time3 = to_prime(time2c,sr);
+				data2.gain1 = gain2a;
+				data2.gain2 = gain2b;
+				data2.gain3 = gain2c;
+			}
+		}
+	}
+	
+	else{	
+		printf("using preset parameters:\n");
+		printf("DATA 1: time1 = %u\n\ttime2 = %u\n\ttime3 = %u\n",data1.time1,data1.time2,data1.time3);
+		printf("DATA 2: time1 = %u\n\ttime2 = %u\n\ttime3 = %u\n",data2.time1,data2.time2,data2.time3);
+		printf("feedback = %.4lf,fb_lpfreq = %.4lf\n",feedback,fb_lpfreq);
+		printf("dry gain = %.4lf, wet gain = %.4lf\n",dry_gain,wet_gain);
+	}
+#endif	
+
+	
+	switch(reverbtype){
+	case(DIFF_SMALL):	
+		printf("rmverb: using small room model\n");
+		s_diffuser = new small_diffuser((unsigned int)(0.024 * sr),&data1,&data2,feedback,fb_lpfreq,sr);
+		s_diffuser->set_damping(double_damping);
+		cdpdiff =  dynamic_cast<cdp_processor *>(s_diffuser);
+		if(!usertaps){
+			ptaps = smalltaps;
+			ntaps = sizeof(smalltaps) / sizeof(deltap);
+		}
+		break;
+	case(DIFF_MEDIUM):
+		printf("rmverb: using medium room model\n");
+		m_diffuser = new medium_diffuser(0.108,&data1,&data2,feedback,fb_lpfreq,sr);
+		m_diffuser->set_damping(double_damping);
+		cdpdiff =  dynamic_cast<cdp_processor *>(m_diffuser);
+		if(!usertaps){
+			ptaps = mediumtaps;
+			ntaps = sizeof(mediumtaps) / sizeof(deltap);
+		}
+		//ptaps = longtaps;
+		//ntaps = sizeof(longtaps) / sizeof(deltap);
+
+		break;
+	case(DIFF_LARGE):
+		printf("using large room model\n");
+		//fb_lpfreq = 2600.0; // Gardner fixed val
+		l_diffuser = new large_diffuser(&data1,&data2,feedback,fb_lpfreq,sr);
+		l_diffuser->set_damping(double_damping);
+		cdpdiff =  dynamic_cast<cdp_processor *>(l_diffuser);
+		if(!usertaps){
+			ptaps = largetaps;
+			ntaps = sizeof(largetaps) / sizeof(deltap);
+		}
+		break;
+	default:
+#ifdef _DEBUG
+		assert(false);
+#endif
+		break;
+	}
+
+	if(!cdpdiff){
+		puts("rmverb: no memory for reverberator\n");
+		exit(1);
+	}
+	if(!cdpdiff->create()){
+		puts("rmverb: no memory for reverberator buffers\n");
+		exit(1);
+	}
+
+
+	/*********** GLOBAL COMPONENTS *************/
+
+	if(lp_freq > 0.0)
+		lp		= new onepole(lp_freq,sr,LOW_PASS);
+	
+	//if user prescribes a predelay, adjust all taptimes
+	if(predelay >= 0.0){
+		double current_delay = ptaps[0].pos;
+		double adjust = predelay - current_delay;
+
+		for(int i = 0; i < ntaps; i++)
+			ptaps[i].pos += adjust;
+
+	}
+	else
+		predelay = ptaps[0].pos;	//to report to user
+	tdelay = new tapdelay(ptaps,ntaps,sr);
+
+	if(!tdelay->create()) {
+		fputs("rmverb: no memory for early reflections\n",stderr);
+		return 1;
+	}
+
+	//max_gain = tdelay->get_maxgain();
+	//if(max_gain > 1.0)
+	//	tdelay->set_gain(max_gain * early_gain);
+
+
+	// final allpass chain: 0.02,0.014,0.009,0.006
+	ap1_l		= new allpassfilter(sr,to_prime(0.02,sr),-0.6,0);
+	ap2_l		= new allpassfilter(sr,to_prime(0.014,sr),0.4,0);
+//	ap3_l		= new allpassfilter(to_prime(0.009,sr),-0.6,0);
+//	ap4_l		= new allpassfilter(to_prime(0.006,sr),0.4,0);
+	if(!(ap1_l->create()
+		&& ap2_l->create()
+//		&& ap3_l->create()
+//		&& ap4_l->create()
+		)){
+		fprintf(stderr,"rmverb: no memory for allpass chain\n");
+		exit(1);
+	}
+	if(outchans >=2){
+		// slightly different params to make stereo!
+		// tune front pair precisely; further chans use formula
+		ap1_r		= new allpassfilter(sr,to_prime(0.025,sr),-0.6,0);
+		ap2_r		= new allpassfilter(sr,to_prime(0.0145,sr),0.4,0);
+//		ap3_r		= new allpassfilter(to_prime(0.0095,sr),-0.6,0);
+//		ap4_r		= new allpassfilter(to_prime(0.0065,sr),0.4,0);
+		if(!(ap1_r->create()
+			&& ap2_r->create()
+//			&& ap3_r->create()
+//			&& ap4_r->create()
+			)){
+			fprintf(stderr,"rmverb: no memory for allpass chain\n");
+			exit(1);
+		}
+	}
+
+	
+
+	printf("using %ld early reflections: predelay = %.2lf msecs\n",ntaps,predelay * 1000.0);
+    int clipfloats = 1;
+    int minheader = 0;
+	if(want_floats) {
+		props.samptype =PSF_SAMP_IEEE_FLOAT;
+        clipfloats = 0;        
+    }
+	props.chans = outchans;
+    
+	if((ofd = psf_sndCreate(argv[OUTFILE],&props,clipfloats,minheader,PSF_CREATE_RDWR)) < 0) {
+        fprintf(stderr,"rmverb: cannot open output file %s\n",argv[OUTFILE]);
+        return 1;
+	}
+
+/****
+ 																		|wet
+ *****   in-> ---o-->pre_lp -->earlies-o--->diffuser-->*diffgain-->--+--*---+----out
+				 |				       |							 |		|
+				 |					   |>--------*-(earlies)->-------|		*--dry
+				 |							     ^earlygain				    |
+				 |__________________________________________________________|
+
+ ****/
+	printf("\ninfile chans = %ld, outfile chans = %ld",chans,outchans);
+	printf("\n");
+	long outframes = 0;
+	long step = (long)(0.25 * sr);
+    float insamps[2], outsamps[16];
+	double l_out,r_out;
+//    float  finsamp, foutsamp;
+	for(;;){
+		int rv;
+		double l_ip,r_ip, out,l_direct,r_direct,mono_direct,earlies = 0.0f;
+        
+		//read mono/left
+		if((rv = psf_sndReadFloatFrames(ifd,insamps,1)) < 0){
+			fprintf(stderr,"rmverb: error reading file\n"); 
+			exit(1);			
+		}		
+		if(rv==0) break;		// end of inputfile	- handle any trailertime
+		l_ip = (double) insamps[0];
+		//apply any conditioning to direct signal
+		if(tc_l_lowcut)
+			l_ip =  tc_l_lowcut->tick(l_ip);
+		if(l_highcut)
+			l_ip =  l_highcut->tick(l_ip);
+		l_direct = l_ip;
+		mono_direct = l_direct;
+		r_direct = l_direct;
+		//handle R channnel if active
+		if(chans==2){
+            r_ip = (double) insamps[1];
+			//apply any conditioning to direct signal			
+			if(tc_r_lowcut)
+				r_ip =  tc_r_lowcut->tick(r_ip);
+			if(r_highcut)
+				r_ip =  r_highcut->tick(r_ip);
+			r_direct = r_ip;
+			if(want_mono)	//input merged sig to reverbreator
+				mono_direct = (l_direct + r_direct) * 0.5;			 
+		}
+		// mono_direct = (mixed) mono input to  reverb		
+		//possibly also filter it...
+		if(lp)
+			mono_direct =  lp->tick(mono_direct);					//input lowpass
+		earlies = tdelay->tick(mono_direct);				//early reflections	
+		//send (filtered) mono output from reverberator
+
+		// TODO: find formula for diffgain...
+		out = earlies + (cdpdiff->tick(earlies + mono_direct) * diff_gain * 0.707);		 //the dense reverberation		
+		// output:: use 2 allpasses per chan to generate stereo reverb, wider diffusion
+		l_out = wet_gain *  ap2_l->tick(ap1_l->tick(out));
+
+		// old method - 1 allpass per channel
+		//l_out = chan_ap[0]->tick(out) * wet_gain;
+
+		if(outchans == 1)			
+			l_out += (mono_direct * dry_gain);
+		else 
+			l_out += l_direct * dry_gain;
+        outsamps[0] = (float) l_out;
+		
+		if(outchans>=2){
+			r_out = wet_gain *  ap2_r->tick(ap1_r->tick(out));			
+			// old version
+			//r_out = chan_ap[1]->tick(out) * wet_gain;
+			r_out += r_direct * dry_gain;
+			outsamps[1] = (float) r_out;
+		}
+
+		//now any further channels; reduced level of direct sig
+		for(i=2;i < outchans; i++){
+            out = earlies + (cdpdiff->tick(earlies + (mono_direct * 0.3)) * diff_gain * 0.707);
+			l_out = wet_gain * chan_ap[i]->tick(out);
+            outsamps[i] = (float) l_out;
+			
+		}
+        if(psf_sndWriteFloatFrames(ofd,outsamps,1) != 1){
+            fprintf(stderr,"rmverb: error writing output file\n"); 
+			return 1;
+        }
+		outframes++;
+		if((outframes % step) == 0)
+			fprintf(stdout,"%.2f\r",(double)outframes/(double)sr);
+	}
+	
+	int trtime  = (int)( trailertime * sr);
+	for(i = 0; i < 	trtime; i++){		
+		double tr_l_ip,tr_r_ip = 0.0f,tr_out,tr_l_direct,tr_r_direct,
+			tr_mono_direct,tr_earlies = 0.0f,tr_l_out,tr_r_out;
+		//need last active samps from input filter(s)
+		tr_l_ip = 0.0;
+		if(tc_l_lowcut)
+			tr_l_ip =  tc_l_lowcut->tick(tr_l_ip);
+		if(l_highcut)
+			tr_l_ip =  l_highcut->tick(tr_l_ip);
+		tr_l_direct = tr_l_ip;
+		tr_mono_direct = tr_l_direct;
+		tr_r_direct = tr_l_direct;
+		//handle R channnel if active
+		if(chans==2){
+			tr_r_ip = 0.0;	   // inout is zero, except if using input filters
+			// get any samples from input conditioning filters
+			if(tc_r_lowcut)
+				tr_r_ip =  tc_r_lowcut->tick(0.0);						   
+			if(r_highcut)
+				tr_r_ip =  r_highcut->tick(tr_r_ip);
+			tr_r_direct = tr_r_ip;
+			if(want_mono)
+				tr_mono_direct = (tr_l_direct + tr_r_direct) *0.5f;			
+		}
+
+		// mono_direct = (mixed) mono input to  reverb
+		//possibly also filter it...
+		if(lp)
+			tr_mono_direct =  lp->tick(tr_mono_direct);		
+		tr_earlies = tdelay->tick(tr_mono_direct);
+		//send (filtered) mono output from reverberator
+		tr_out = tr_earlies + (cdpdiff->tick(tr_earlies + tr_mono_direct) * diff_gain * 0.707);		
+		// output:: use 2 allpasses to generate stereo reverb, further diffusion
+
+		// new version
+		tr_l_out = wet_gain *  ap2_l->tick(ap1_l->tick(tr_out));
+		// old version		
+		//tr_l_out = chan_ap[0]->tick(tr_out) * wet_gain;
+
+		if(want_mono)
+			tr_l_out += (tr_mono_direct * dry_gain);
+		else			
+			tr_l_out += tr_l_direct * dry_gain;
+
+        outsamps[0] = (float) tr_l_out;
+		
+		if(outchans>=2){						
+			tr_r_out = wet_gain * ap2_r->tick(ap1_r->tick(tr_out));
+			// old version:
+			//tr_r_out = chan_ap[1]->tick(tr_out) * wet_gain;
+
+			tr_r_out += tr_r_direct * dry_gain;	   // still need this cos of input filters
+            outsamps[1] = (float) tr_r_out;
+		}
+		for(int j = 2; j < outchans; j++){			
+			tr_r_out = chan_ap[j]->tick(tr_out) * wet_gain;			
+            outsamps[j] = (float) tr_r_out;
+			
+		}
+        if(psf_sndWriteFloatFrames(ofd,outsamps,1) != 1){
+            fprintf(stderr,"rmverb: error writing output file\n"); 
+			return 1;
+        }
+		outframes++;
+		if((outframes % step)==0)
+			//inform(step,sr);
+			fprintf(stdout,"%.2f\r",(double)outframes/(double)sr);
+	}
+	fprintf(stdout,"%.2f\n",(double)outframes/(double)sr);
+	//stopwatch(0);
+    
+	peaks = (PSF_CHPEAK *) malloc(sizeof(PSF_CHPEAK) * outchans);
+    if(peaks && psf_sndReadPeaks(ofd,peaks,NULL)) {
+        printf("\nPEAK data:\n");
+        for(i=0;i < outchans;i++)
+            printf("Channel %d: %.4f at frame %u: %.4lf secs\n",i+1,
+                   peaks[i].val,peaks[i].pos, (double) peaks[i].pos / (double)props.srate);
+	}
+	else {
+        fputs("rmverb: warning: unable to read PEAK data\n",stderr);
+    }
+	
+	psf_sndClose(ifd);
+	
+	if(psf_sndClose(ofd) < 0)
+		fprintf(stderr,"rmverb: error closing outfile\n");
+    printf("\n");
+    
+    psf_finish();
+	free(peaks);
+	
+	delete tdelay;
+	
+	if(usertaps)
+		delete [] ptaps;
+	//can I call delete on base-class ptr?
+	if(s_diffuser)
+		delete s_diffuser;
+	if(m_diffuser)
+		delete m_diffuser;
+	if(l_diffuser)
+		delete l_diffuser;
+	if(lp)
+		delete lp;
+	if(outchans > 1){
+		for(i=0;i < outchans; i++)
+			delete chan_ap[i];
+		delete [] chan_ap;
+	}	
+	delete tc_l_lowcut;	
+	delete tc_r_lowcut;	
+	delete l_highcut;
+	delete r_highcut;
+	return 0;
+}
+
+
+void
+usage(void)
+{
+	fprintf(stderr,"\nrmverb: Multi-channel reverb with room simulation\n\tVersion 1.3 Release 5 CDP 1998,1999,2011");
+	fprintf(stderr,
+	"\nusage:\nrmverb [flags] infile outfile rmsize rgain mix fback absorb lpfreq trtime"
+	"\nflags    : any or all of the following"	
+	"\n   -cN   : create outfile with N channels (1 <= N <= 16 : default =2)"
+	"\n   -d    : apply double lowpass damping (internal to nested allpasses)"
+	"\n            (see <absorb>. NB: reduces reverb time: "
+	"\n                  increase feedback to compensate)"
+	"\n   -eFILE: read early reflections data from breakpoint textfile FILE"
+	"\n   -f    : force floating-point output (default: infile format)"	
+	"\n   -HN   : apply 12dB Highcut filter with cutoff freq N Hz) to infile"
+	"\n   -LN   : apply 12dB Lowcut filter with cutoff freq N Hz to infile"
+	"\n           (NB: currently fixed at 12dB/oct - could be made variable)\n"
+	"\n   -pN   : force reverb predelay to N msecs (shifts early reflections)"
+	"\nrmsize   : 1,2 or 3: small, medium or large"
+	"\nrgain    : set level of dense reverb (0.0 < egain <= 1.0)"
+	"\nmix      : dry/wet balance (source and reverb). 1.0<-- mix -->=0.0"
+	"\nfback    : reverb feedback level: controls decay time. 0.0 <= fback <= 1.0"
+	"\nabsorb   : cutoff freq (Hz) of internal lowpass filters (models air absorption)\n"
+	"                (typ: 2500 for large room --- 4200 for small room "
+	"\nlpfreq   : lowpass filter cutoff freq(Hz) applied at input to reverb"	
+	"\n           NB: to disable either filter (absorb,lpfreq), use the value 0"
+	"\ntrtime   : trailer time added to outfile for reverb tail (secs)\n");	
+	exit(0);
+}
+
+long readtaps(FILE *fp, deltap **taps,double sr)
+{
+	vector<deltap> vtaps;
+	deltap thistap;
+	char line[256];
+	int ntaps = 0;
+	int linecount = 1;
+	bool error = false;
+	double time= 0.0, current_time = 0.0, val = 0.0;
+
+	double sample_duration = 1.0 / sr;
+	if(fp==0 || sr <= 0.0){
+		*taps = 0;
+		return 0;
+	}
+	while(fgets(line,256,fp))	{
+		int got;
+		if(line[0] == '\n' || line[0]== ';'){	//allow comment lines; this does not check for leading white-space...
+			linecount++;
+			continue;
+		}
+
+		if((got = sscanf(line,"%lf%lf",&time,&val))<2){			
+			fprintf(stderr,"\nerror in reflections file: line %d",linecount);
+			error =true;
+			break;
+		}
+		if(time < 0.0){
+			fprintf(stderr,"\nerror in tapfile: line %d: bad time value",linecount);
+			error = true;
+			break;
+		}
+		//if (first) val = 0.0, or VERY close
+		if(time < sample_duration)	{ //non-zero time must be at least one sample on!			
+			fprintf(stderr,"\nWARNING: very small taptime %.4lf treated as zero",time);
+			time = 0.0;			
+		}
+
+		if(current_time != 0.0 && time==current_time){
+			fprintf(stderr,"\nWARNING: duplicate time in line %d: ignored",linecount);
+			linecount++;
+			continue;
+		}
+		if(time < current_time){
+			fprintf(stderr,"\nerror in tapfile: time out of sequence: line %d",linecount);
+			error = true;
+			break;
+		}
+		current_time = time;
+		thistap.pos = time;
+		if(val <=0.0 || val > 1.0){
+			fprintf(stderr,"\nerror: bad amplitude in tapfile: line %d",linecount);
+			error = true;
+			break;
+		}
+		thistap.val = val;
+		vtaps.push_back(thistap);
+		linecount++;
+	}
+	
+	ntaps = vtaps.size();
+	if(ntaps==0){
+		fprintf(stderr,"\ntapfile contains no data!");
+		error = true;
+	}
+	if(error){
+		*taps = 0;
+		return 0;
+	}
+	*taps = new deltap[ntaps];
+	if(taps==0)
+		return 0;
+	vector<deltap>::iterator I =  vtaps.begin(); 
+	for(int i = 0; I != vtaps.end();i++, I++){
+		(*taps)[i].pos = I->pos;
+		(*taps)[i].val = I->val;
+	}
+
+	return ntaps;
+}

dev/externals/reverb/roomresp.cpp

@@ -0,0 +1,713 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+/***********************************************************************
+Motorola
+
+Motorola does not assume any liability ...
+Author:   Tom Zudock ([email protected])
+
+see:
+ http://collaboration.cmc.ec.gc.ca/science/rpn/biblio/ddj/Website/articles/DDJ/1996/9612/9612d/9612d.htm
+ 
+Presented in Dr Dobbs Journal Vol21:12, December 1996.
+ 
+Filename: roomresp.cpp
+
+This program calculates the room response for a simple room.  The output is not
+particularly delightful to listen to (it could sound phasey and
+have ringing due to the simple model used), but this program is the
+foundation for more complex virtual audio models which will eliminate
+these anomalies.
+
+Example form of the InputParametersTextFile:
+
+0.5               # Input gain for data file (float)
+1                 # Open Path = 0 Closed Path = 1 (int)
+6    5    4       # Room Dimensions (X,Y,Z) (float)
+0.85              # Reflectivity of walls (float)
+5                 # Number of rooms per dimension, 1 or greater
+2    3    2       # Listener coord (X,Y,Z) (float)
+2    2    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+3    2    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+4    3    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+3    4    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+2    4    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+
+***********************************************************************/
+
+// included header files
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+//#include <process.h>
+#include <string.h> 
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+#include <vector>
+#include <algorithm>
+using namespace std;
+
+#ifndef _MAX
+#define _MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+
+// macro defined constants
+//#define Vs 334.0
+static const double Vs = 334.0;
+static const double srate = 20000;
+#define PI 3.14159265359
+#define RefDist 1.0
+enum args {
+		REFLECTIVITY,
+		REFLECTIONS,
+		R_X,
+		R_Y,
+		R_Z,
+		S_X,
+		S_Y,
+		S_Z,
+		L_X,
+		L_Y,
+		L_Z
+};
+
+// structure for a point in 3 dimensional space
+struct point {
+	double X;
+	double Y;
+	double Z;
+	point(){
+		X = Y = Z = 0.0;
+	}
+	virtual ~point() {}
+};
+
+// structure for a point paired with a time duration
+struct position {
+	point Coord;
+	double Time;
+	position *NextSource;
+};
+
+const char* cdp_version = "5.0.0";
+
+bool input_cmds(char *args[],int argc,point &room,point &source,point &listener,double &reflcoef,int &NR);
+
+void
+CalcImpResp(
+    long *,
+    double *,
+    point &,
+    point &,
+    point &,
+    double,
+    double,
+    int);
+
+position *
+InputParameters(
+    FILE *,
+    point &,
+ //   point &,
+    point &,
+    double &,
+    double &,
+   int &);
+
+void
+InputParamsError(
+    int);
+
+//RWD stuff
+typedef struct deltap {
+		long time;
+		double val;
+	} DELTAP;
+
+bool operator<(const DELTAP &a,const DELTAP &b)
+{
+	return (a.time < b.time);
+}
+
+void usage()					//RWD.1..12.98 left out flag -nNORM for now...
+{
+	printf("\n*********  ROOMRESP.EXE: by Tom Zudock  CDP build 1998,1999 ******"
+		   "\n  GENERATE ROOM-RESPONSE EARLY REFLECTIONS DATAFILE\n"
+		   "\nusage: roomresp [-aMAXAMP][-rRES] txtout.dat liveness nrefs "
+		   "\n                          roomL roomW roomH"
+		   "\n                             sourceL sourceW sourceH"
+		   "\n                                listenerL listenerW listenerH");
+	printf("\ntxtout.dat: text outfile containing early reflections in brkpoint format"
+			"\n             suitable for input to rmverb, reverb or tdelay"
+			"\nMAXAMP   : peak amp for data (0.0 < MAXAMP <= 1.0; default 1.0)"
+			"\nRES      : time resolution for reflections in msecs"
+			"\n             (0.1 <= RES <= 2; default 0.1)"
+			"\nliveness : degree of reflection from each surface (0 to 1; typ: 0.95)"
+			"\nnrefs    : number of reflections from each surface(nrefs > 0; typ: 2 to 5)"
+			"\n           (WARNING: high values will create extremely long data files!)"
+			"\nROOM PARAMETERS:"
+			"\nroomL roomW roomH            : room size: (Length, Width, Height)"
+			"\nsourceL sourceW ourceH       : position of sound source (as above)"
+			"\nlistenerL listenerW listenerH: position of listener (as above)"
+			"\n                               all dimensions are in metres"
+			"\nNB: first output time is non-zero."
+			);
+}
+
+
+/***********************************************************************
+main
+
+The main program is responsible for the primary program flow.
+It contains the primary processing loop which calls subroutines
+to handle the specialized functionality.
+***********************************************************************/
+int main (int argc, char **argv) {
+	
+//	FILE        *PositionFile = 0;  // input text parameters file
+	//RWD to print earlies as text
+	FILE		*earliesFile = 0;
+#ifdef TESTING
+	FILE		*rawfile = 0;				///just for testing...raw output from impresp funtionm
+#endif
+	int			i,NR = 0;							// num of mirrored rooms per dimension
+	long		*TapDelay = 0;			// pointer to array of tap delays	
+	double		Fs = 0.0;							// the sampling frequency
+	double		*TapGain = 0;				// pointer to array of tap gains
+	double		ReflCoef = 0.0;				// the reflectivity of the walls
+//	position	*CurrentSource = 0;	// current source in positions list 
+	point		Room;						// coords for the room size
+	point		Source;					// coords of current sound source
+	point		Listener;				// coords of the listener
+#ifdef _DEBUG
+	bool		write_c_struct = false;
+#endif
+	bool		do_normalize = true;
+	double		normfac = 1.0;
+	double      min_separation = 0.0001;
+
+	if((argc==2) && strcmp(argv[1],"--version")==0) {
+		fprintf(stdout,"%s\n",cdp_version);
+		fflush(stdout);
+		return 0;
+	}
+	if(argc==1){
+		usage();
+		exit(0);
+	}
+	if (argc < 13){
+		fprintf(stderr,"\ninsufficient arguments\n");
+  		usage();
+		exit(-1);
+	}
+  
+	while(argc > 0 && argv[1][0]=='-'){
+		double val;
+		switch(argv[1][1]){
+#ifdef _DEBUG
+//my private facility...
+		case('c'):
+			write_c_struct = true;
+			printf("\nwriting C-format file");
+			break;
+#endif
+		case('a'):
+			if(argv[1][2]=='\0'){
+				fprintf(stderr,"\n-a flag requires a value\n");
+				usage();
+				exit(1);
+			}
+			val = atof(&(argv[1][2]));
+#ifdef _DEBUG
+			if(val <0.0 || val > 1.0){
+#else
+			if(val <=0.0 || val > 1.0){
+#endif
+				fprintf(stderr,"\n normalization value %lf out of range\n",val);
+				usage();
+				exit(1);
+			}
+#ifdef _DEBUG
+			if(val==0.0) {
+				do_normalize = false;
+				normfac = 1.0;
+				printf("\nskipping normalization");
+			}
+			else
+#endif
+				normfac = val;
+			break;
+		case('r'):
+			if(argv[1][2]=='\0'){
+				fprintf(stderr,"\n-r flag requires a value\n");
+				usage();
+				exit(1);
+			}
+			val = atof(&(argv[1][2]));			
+			if(val < 0.1 || val > 2.0) {
+				fprintf(stderr,"\nvalue %.4lf for RES out of range\n",val);
+				usage();
+				exit(1);				
+			}
+			else
+				min_separation = val * 0.001;
+			break;
+
+
+
+		default:		
+			fprintf(stderr,"\nbad flag option\n");
+			exit(1);
+			break;
+		}
+		argc--;
+		argv++;
+	}
+
+	//step along arglist:
+	argc--; argv++;
+
+	
+	if(argc != 12){
+		fprintf(stderr,"\nincorrect number of numeric arguments\n");
+		usage();
+		exit(1);
+	}	
+	if((earliesFile = fopen(argv[0],"w"))==NULL){
+		printf ("\nUnable to open %s for output\n",argv[2]);
+		exit(-1);
+	}
+#ifdef TESTING
+	if((rawfile = fopen("rawdata.txt","w"))==NULL){
+		printf ("Unable to open rawdata.txt for output\n");
+		exit(-1);
+	}
+#endif
+	//step along again...
+	argc--; argv++;
+
+	//NB this sets the efective time-resolution for each tap, ie 1.0 / fs
+	//so 1000 is a minimum, for 1msec resolution
+
+	Fs = srate;
+
+	if(!input_cmds(argv,argc,Room,Source,Listener,ReflCoef,NR))
+		exit(1);	//already got the errmsg
+  
+	if((TapDelay = new long[NR*NR*NR])==NULL) {
+    	printf("Error: Unable to allocate memory\n");
+    	exit(-1);
+    }
+    if((TapGain = new double[NR*NR*NR])==NULL) {
+    	printf("Error: Unable to allocate memory\n");
+    	exit(-1);
+    }
+	
+  
+	CalcImpResp(
+	    TapDelay,
+		TapGain,
+		Room,
+		/*CurrentSource->Coord,*/Source,
+		Listener,
+		ReflCoef,
+		Fs,
+		NR);
+
+	//sort deltaps, eliminate duplicates,  write to file
+	vector<DELTAP> vtaps;
+	DELTAP thistap,prevtap;
+	prevtap.time = 0;
+	prevtap.val = 0.0;
+	for(i = 0; i < NR*NR*NR;i++) {
+		thistap.time = TapDelay[i];
+		thistap.val  =TapGain[i];
+#ifdef TESTING
+		fprintf(rawfile,"%ld\t%.6lf\n",thistap.time,thistap.val);
+#endif
+		//exclude tapval of 1.0: = direct sig at listening point
+		if(thistap.val > 0.0 && thistap.val < 1.0)			// get a problem in release build otherwise....
+			vtaps.push_back(thistap);	
+	}
+
+	//now we can sort them
+
+	sort(vtaps.begin(),vtaps.end());
+	
+	vector<DELTAP>::iterator it;
+	it = vtaps.begin();
+	prevtap.time = it->time;
+	prevtap.val =  it->val;
+
+	//print all taps as we got them, first!
+
+	double maxval = 0.0;
+	for(;it != vtaps.end();it++){
+#ifdef NOTDEF		
+		fprintf(earliesFile,"\n\t%.4lf\t%.6lf",(double)it->time / Fs,it->val);
+#endif
+		maxval = _MAX(it->val,maxval);
+		
+	}
+#ifdef _DEBUG
+	printf("\noverall maxamp = %.4lf",maxval);
+#endif
+
+	//print first line, formatted for code
+	//normalize all vals to max amplitude!
+	
+	if(do_normalize){
+		normfac  /= prevtap.val;
+#ifdef _DEBUG
+		printf("\nnormalizing by %.4lf",normfac);
+#endif
+	}
+	it = vtaps.begin();
+	
+	
+	//normfac= floor(normfac);
+#ifdef _DEBUG
+	printf("\ninitial time=  %ld, initial amp = %.4lf",it->time,it->val);
+
+	//RWD private facility!
+	if(write_c_struct){
+		fprintf(earliesFile,"\n\n{\n\t{%.4lf,%.6lf}",(double)prevtap.time / Fs,prevtap.val * normfac);
+		int count = 1;
+	
+		for(;it != vtaps.end();it++){
+			//get maxamp of any taps sharing a time (to 5 sig figs) - or should I add them somehow?
+			if(fabs(it->time - prevtap.time) < min_separation){
+				it->val = _MAX(it->val,prevtap.val);
+				continue;
+			}
+			fprintf(earliesFile,",\n\t{%.5lf,%.6lf}",(double)it->time / Fs,it->val * normfac);
+			prevtap = *it;
+			count++;
+		}
+		fprintf(earliesFile,"\n\t};\n\nNUMTAPS = %d\n",count);
+	}
+	else {
+#endif
+		int count = 1;
+		fprintf(earliesFile,"%.5lf\t%.6lf\n",(double)prevtap.time / Fs, prevtap.val * normfac);
+		for(;it != vtaps.end();it++){
+			if(fabs(it->time - prevtap.time) < min_separation){
+				it->val = _MAX(it->val,prevtap.val);
+				continue;
+			}
+			fprintf(earliesFile,"%.5lf\t%.6lf\n",(double)it->time / Fs, it->val * normfac);
+			prevtap = *it;
+			count++;
+		}
+		printf("\n\nNUMTAPS = %d\n",count);
+#ifdef _DEBUG
+	}
+#endif
+	fclose(earliesFile);
+#ifdef TESTING
+	fclose(rawfile);
+#endif
+	delete TapDelay,
+	delete TapGain,
+	printf ("done\n");
+	return(0);
+}               
+
+/***********************************************************************
+CalcImpResp
+
+This subroutine calculates the time delays (in samples) and
+attenuations for the early reflections in the room.
+***********************************************************************/
+void
+CalcImpResp(
+  long * TapDelay,
+  double * TapGain,
+  point & Room,
+  point & Source,
+  point & Listener,
+  double ReflCoef,
+  double Fs,
+  int NR)						
+{
+	double Dist = 0.0;                       // distance travelled by sound ray
+	double ReflOrd = 0.0;                    // reflection order of sound ray
+	//RWD.11.98 need volatile, or VC++ optimizer does more bad things!
+	volatile point MirrSource;                  // mirrored source x,y,z coords
+	//RWD
+	int NRovr2 = NR / 2;
+	int index = 0;
+  for (int i=-NRovr2;i<=NRovr2;i++) {    // loop through all 3 dimensions
+    for (int j=-NRovr2;j<=NRovr2;j++) {
+      for (int k=-NRovr2;k<=NRovr2;k++) {
+        // calc x,y,z sound source coords in mirrored room		  
+        MirrSource.X = (i)*Room.X
+					   +fabs(((i)%2)*Room.X)
+                       +pow(-1,(i))*Source.X;                       
+        MirrSource.Y = (j)*Room.Y
+                       +fabs(((j)%2)*Room.Y)
+                       +pow(-1,(j))*Source.Y;
+        MirrSource.Z = (k)*Room.Z
+                       +fabs(((k)%2)*Room.Z)
+                       +pow(-1,(k))*Source.Z;
+        // calculate distance to listener
+        Dist = sqrt(pow(MirrSource.X-Listener.X,2)
+                    +pow(MirrSource.Y-Listener.Y,2)
+                    +pow(MirrSource.Z-Listener.Z,2));
+        // calculate delayed arrival time of reflection
+		//RWD problem with release build...
+		
+		index = (i+NRovr2)*NR*NR+(j+NRovr2)*NR+(k+NRovr2);
+		//this also stopped the optimizer doing bad things....
+		//if(index%NR==0)
+		//	printf("\nindex = %d\tX= %.4lf\tY=%.4lf\tZ=%.4lf\tDist = %.4lf",
+		//		index,MirrSource.X,MirrSource.Y,MirrSource.Z,Dist);
+																					
+        TapDelay[index]=Dist/Vs*Fs;		
+        ReflOrd = abs(i)+abs(j)+abs(k);
+        // calculate attenuation for the reflection
+		//RWD div/zero avoidance:
+		if(Dist==0.0)
+			TapGain[index] = 1.0;
+		else
+			TapGain[index]= pow(RefDist/Dist,2.0)*pow(ReflCoef,ReflOrd);		
+      }
+    }
+  }
+}
+
+/***********************************************************************
+InputParameters
+
+This subroutine inputs parameters from a text file.  An example of
+the format of the parameters file is below.  It also checks the
+range of some of the parameters and creates a linked list representing
+how the sound source changes position in the modeled room.
+
+0.5               # Input gain for data file (float)
+1                 # Open Path = 0 Closed Path = 1 (int)
+6    5    4       # Room Dimensions (X,Y,Z) (float)
+0.85              # Reflectivity of walls (float)
+5                 # Number of rooms per dimension, 1 or greater
+2    3    2       # Listener coord (X,Y,Z) (float)
+2    2    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+3    2    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+4    3    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+3    4    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+2    4    2    1  # Source coord and time at coord (X,Y,Z,t) (float)
+
+***********************************************************************/
+position*
+InputParameters(
+  FILE *PositionFile,
+  point &Room,
+//  point &Source,
+  point &Listener,
+  double &ReflCoef,
+  double &InGain,
+  int &NR)
+{
+	position *First = new position,*Previous=First,*Next=Previous;
+	int ClosedPath,LineCount=1;
+	char StrLine[256];
+
+	#define ReadParms(arg1)                           \
+	if ((fgets(StrLine,256,PositionFile)) != NULL) {  \
+		if (arg1)                                       \
+			LineCount++;                                  \
+		else                                            \
+			InputParamsError(LineCount);                  \
+	}                                                 \
+	else                                              \
+			InputParamsError(0);
+		
+	
+#ifdef _DEBUG
+	assert(PositionFile);
+#endif
+	ReadParms(sscanf(StrLine,"%le",&InGain)==1)
+	ReadParms(sscanf(StrLine,"%d",&ClosedPath)==1)
+	ReadParms(sscanf(StrLine,"%le %le %le",&Room.X,&Room.Y,&Room.Z)==3)
+	ReadParms(sscanf(StrLine,"%le",&ReflCoef)==1)
+	ReadParms(sscanf(StrLine,"%d",&NR)==1)
+	ReadParms(sscanf(StrLine,"%le %le %le",&Listener.X,&Listener.Y,&Listener.Z)==3)
+
+	if (NR%2 != 1) {
+		printf ("Error: The number of reflected rooms per dimension must be odd (for symmetry)\n");
+		exit(-1);
+	}
+	//RWD!!! was bitwise OR:  '|'
+	if ((Listener.Y > Room.Y)||(Listener.X > Room.X)) { 					// Check listener position
+		printf ("Error: Listener position puts listener outside of room \n");
+		exit(-1);
+	}		
+  // Choose a default position for listener
+	//RWD check this! am I always getting these setings????
+	First->Coord.X = Room.X/2;
+	First->Coord.Y = Room.Y/2;
+	First->Coord.Z = Room.Z/2;
+	First->Time = 0;
+	First->NextSource=NULL;
+	while (feof(PositionFile)==0) {                                              
+		if ((fgets(StrLine,256,PositionFile)) != NULL) {
+			if(sscanf(StrLine,"%le %le %le %le",&(Next->Coord.X),&(Next->Coord.Y),&(Next->Coord.Z),&(Next->Time))==4) {
+				Next->NextSource = new position;
+				Previous = Next;
+				Next = Next->NextSource;
+				LineCount++;
+			}
+			else {
+				InputParamsError(LineCount);
+			}
+		}
+	}
+	delete (Previous->NextSource);
+	if (ClosedPath)
+		Previous->NextSource=First;
+	else
+		Previous->NextSource=NULL;
+	return(First);
+}
+
+//RWD get params from cmdline: refl,nrefs,rX,rY,rZ,sX,sY,sZ,lX,lY,lZ
+//fix ingain at 1.0
+bool input_cmds(char *args[],int argc,point &room,point &source,point &listener,double &reflcoef,int &NR)
+{
+	double dval;
+	int ival;
+#ifdef _DEBUG
+	assert(argc==11);
+#endif
+	if(argc != 11)
+		return false;
+
+	dval = atof(args[REFLECTIVITY]);
+	if(dval <= 0.0 || dval > 1.0){
+		fprintf(stderr,"\nbad value for liveliness\n");
+		return false;
+	}
+	reflcoef = dval;
+
+	ival = atoi(args[REFLECTIONS]);
+	if(ival < 1){
+		fprintf(stderr,"\nbad value for nrefs\n");
+		return false;
+	}
+
+	NR = 1 + (2* ival);
+
+	room.X = atof(args[R_X]);
+	if(room.X <= 0.0){
+		fprintf(stderr,"\nbad size for roomX\n");
+		return false;
+	}
+	room.Y = atof(args[R_Y]);
+	if(room.Y <= 0.0){
+		fprintf(stderr,"\nbad size for roomY\n");
+		return false;
+	}
+	room.Z = atof(args[R_Z]);
+	if(room.Z <= 0.0){
+		fprintf(stderr,"\nbad size for roomZ\n");
+		return false;
+	}
+	listener.X = atof(args[L_X]);
+	if(listener.X < 0.0){
+		fprintf(stderr,"\nbad size for listenerX\n");
+		return false;
+	}
+	listener.Y = atof(args[L_Y]);
+	if(listener.Y < 0.0){
+		fprintf(stderr,"\nbad size for listenerY\n");
+		return false;
+	}
+	listener.Z = atof(args[L_Z]);
+	if(listener.Z < 0.0){
+		fprintf(stderr,"\nbad size for listenerZ\n");
+		return false;
+	}
+
+	if ((listener.Y > room.Y)||(listener.X > room.X)|| (listener.Z > room.Z)) { // Check listener position
+		printf ("\nError: Listener position puts listener outside room\n");
+		return false;
+	}
+
+	source.X = atof(args[S_X]);
+	if(source.X <= 0.0){
+		fprintf(stderr,"\nbad size for sourceX\n");
+		return false;
+	}
+
+	source.Y = atof(args[S_Y]);
+	if(source.X <= 0.0){
+		fprintf(stderr,"\nbad size for sourceY\n");
+		return false;
+	}
+
+	source.Z = atof(args[S_Z]);
+	if(source.X <= 0.0){
+		fprintf(stderr,"\nbad size for sourceX\n");
+		return false;
+	}
+
+	if ((source.Y > room.Y)||(source.X > room.X)|| (source.Z > room.Z)) { // Check source position
+		printf ("\nError: Source position puts source outside room\n");
+		return false;
+	}
+
+
+	return true;
+}
+
+
+
+/***********************************************************************
+InputParamsError
+
+This subroutine outputs an text error message to standard output
+when an error is encountered in the InputParams subroutine.  The 
+error text that is output is selected by the line number of the
+input parameters file where the error occured.
+***********************************************************************/
+void
+InputParamsError(
+  int ErrNum)
+{
+	char ErrorStr[][256] = {
+	"Error: EOF reached before all parameters input\n",
+	"Error: Line %d requires 1 parameter the gain for the input file\n",
+	"Error: Line %d requires 1 parameter (0 or 1): 0=closed path, 1=open path\n",
+	"Error: Line %d requires 3 parameters (X,Y,Z): the room dimensions\n",
+	"Error: Line %d requires 1 parameter the wall refelctivity\n",
+	"Error: Line %d requires 1 parameter the number of mirrored rooms per dimension\n",
+	"Error: Line %d requires 3 parameters (X,Y,Z): the listeners position\n",
+	"Error: Line %d requires 4 parameters (X,Y,Z,t): sound position position and time\n"};
+	if (ErrNum < 7)
+		printf(ErrorStr[ErrNum],ErrNum);
+	else
+		printf(ErrorStr[7],ErrNum);
+	exit(-1);
+
+}
+
+ 

dev/externals/reverb/smalldiff.cpp

@@ -0,0 +1,95 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+//smalldiff.cpp
+//implements diffuse reverberator for small rooms, based on gardner nested allpasses
+#include "reverberator.h"
+
+small_diffuser::small_diffuser(unsigned int pre_delay, 
+							   const NESTDATA *apdata1,
+							   const NESTDATA *apdata2,double fb_gain,double lp_coeff){
+
+	  ap1_data = *apdata1;
+	  ap2_data = *apdata2;
+	  predelay_time = pre_delay;
+	  lpgain = fb_gain;
+	  lpcoeff = lp_coeff;
+	  ap1 = ap2 = 0;
+	  lp1 = 0;
+	  predelay = 0;
+	  out1 = out2 = 0.0f;
+}
+
+small_diffuser::~small_diffuser()
+{
+	delete ap1;
+	delete ap2;
+	if(predelay)
+		delete predelay;
+	delete lp1;
+}
+
+bool small_diffuser::create(void)
+{
+	if(ap1_data.time1 == 0 || ap1_data.gain1 <= 0.0)
+		return false;
+	if(ap1_data.time2 ==0 || ap1_data.gain2 <= 0.0)
+		return false;
+	if(ap2_data.time1 == 0 || ap2_data.gain1 < 0.0)
+		return false;
+	if(ap2_data.time2 ==0 || ap2_data.gain2 < 0.0)
+		return false;
+	ap1 = new nested_allpass(ap1_data.time1,ap1_data.time2,ap1_data.time3,
+							ap1_data.gain1,ap1_data.gain2,ap1_data.gain3,ap1_data.delay1,ap1_data.delay2);
+	if(!ap1->create())
+		return false;
+	if(ap2_data.gain1 != 0.0){	 //allow ap to eliminate second block
+		ap2 = new nested_allpass(ap2_data.time1,ap2_data.time2,ap2_data.time3,
+							ap2_data.gain1,ap2_data.gain2,ap2_data.gain3,ap2_data.delay1,ap2_data.delay2);
+		if(!ap2->create())
+			return false;
+	}
+	if(predelay_time > 0){
+		predelay = new delay(predelay_time,1.0);
+		if(!predelay->create())
+			return false;
+	}
+	if(lpcoeff > 0.0)
+		lp1 = new lowpass1(lpcoeff);
+	return true;
+}
+
+
+float small_diffuser::tick(float insam) 
+{
+	float ip = insam;
+	double temp;
+	if(lp1)
+		temp = lp1->tick((out2 + out1)*0.5f);
+	else
+		temp = (out2 + out1)*0.5f;
+	ip += (float)(lpgain * temp);	
+	if(predelay)
+		ip = predelay->tick(ip);
+	out1 = ap1->tick(ip);
+	if(ap2)
+		out2 = ap2->tick(out1);
+	return (out1 + out2)  /** 0.5f*/;
+}

dev/externals/reverb/tdelaymain.cpp

@@ -0,0 +1,688 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+//tdelay.cpp
+//mono/stereo tapped delay line
+//RWD,CDP 1998
+//VERSION: initial, 0.01
+//usage: tdelay [-f] infile outfile feedback tapfile.txt [trailertime]
+//		-f : floatsam output
+//		feedback (double) - obvious...
+// tapfile.txt	lines of at least two values: time  amplitude  [pan]
+// times must be increasing, non-identical
+// amplitude +- 1.0
+// pan -1.0--1.0: -1 = hard Left, 1= hard Right, 0.0 = Centre
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+extern "C" {
+#include <sfsys.h>
+}
+#include <props.h>
+#include "reverberator.h"
+#include <vector>
+#include <algorithm>
+using namespace std;
+
+typedef struct stereo_tap {
+	double time;
+	double amp;
+	double pan;
+} STAP ;
+
+
+void usage(void);
+#define ROOT2	(1.4142136)
+enum pan_direction {SIGNAL_TO_LEFT,SIGNAL_TO_RIGHT};
+
+enum args {INFILE=1,OUTFILE,TAPGAIN,FEEDBACK,MIX,TAPFILE,TRTIME};
+
+void pancalc(double position,double *leftgain,double *rightgain);
+
+const char* cdp_version = "5.0.0";
+
+
+int main(int argc, char **argv)
+{
+	int ifd,ofd,inchans;
+	float feedback = 0.0f;
+	float dryfac = 0.0f,wetfac = 1.0f;
+	double trailertime = 0.0;
+	double current_time = 0.0;
+	deltap *l_taps = 0,*r_taps = 0;
+	tapdelay *l_dline = 0, *r_dline = 0;
+	unsigned int i,nchans = 1,ntaps = 0,trailersamps = 0;
+	bool is_stereo = false;
+	bool floatsams_wanted = false;
+	bool mix_input = true;
+	SFPROPS props;
+	FILE *tapfile = 0;
+	double sample_duration,tapgain;
+	double max_lgain = 1.0,max_rgain= 1.0;
+	float max_gain = 1.0f;
+	double sr; 
+	CHPEAK *peaks;
+	
+	if((argc==2) && strcmp(argv[1],"--version")==0) {
+		fprintf(stdout,"%s\n",cdp_version);
+		fflush(stdout);
+		return 0;
+	}
+
+	if(argc < TRTIME){
+		usage();
+		exit(1);
+	}
+	if(sflinit("tdelay")){
+		sfperror("tdelay: initialisation");
+		exit(1);
+	}
+	
+	while(argv[1][0] == '-'){
+		switch(argv[1][1]){
+		//case('s'):
+		//	is_stereo = true;
+		//	break;
+		case('f'):
+			floatsams_wanted = true;
+			break;
+		default:
+			fprintf(stderr,"\nillegal flag option %s",argv[1]);
+			usage();
+			exit(1);
+			break;
+		}
+		argv++;
+		argc--;
+	}
+	if(argc < TRTIME){
+		usage();
+		exit(1);
+	}
+
+	if((ifd = sndopenEx(argv[INFILE],0,CDP_OPEN_RDONLY)) < 0) {
+		fprintf(stderr,"reverb: cannot open input file %s: %s\n", argv[1],rsferrstr);
+		exit(1);
+	}
+	if(!snd_headread(ifd,&props)){
+		fprintf(stderr,"\nerror reading infile header: %s\n",rsferrstr);
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	if(props.type != wt_wave){
+		fprintf(stderr,"\ninfile is not a waveform file");
+		sndcloseEx(ifd);
+		exit(1);
+	}
+
+	inchans = props.chans;
+	if(inchans > 2){
+		fprintf(stderr,"\nSorry - cannot process files with more than two channels!");
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	sample_duration = 1.0 / props.srate;
+
+	if(floatsams_wanted)
+		props.samptype =  FLOAT32;
+	
+
+	tapgain = atof(argv[TAPGAIN]);
+	if(tapgain <= 0.0){
+		fprintf(stderr,"\nError: tapgain must be > 0.0\n");
+		exit(1);
+	}
+
+
+	feedback = (float) atof(argv[FEEDBACK]);
+	if(feedback < -1.0f || feedback > 1.0f){
+		fprintf(stderr,"\nfeedback out of range	-1.0 to +1.0");
+		usage();
+		sndcloseEx(ifd);
+		exit(1);
+	}
+
+	dryfac = (float) atof(argv[MIX]);
+	if(dryfac < 0.0f || dryfac >= 1.0f){
+		fprintf(stderr,"\nmix value %s out of range",argv[4]);
+		usage();
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	wetfac = 1.0f - dryfac;
+
+	if(argc==TRTIME+1){
+		trailertime = atof(argv[TRTIME]);
+		if(trailertime < 0.0){
+			fprintf(stderr,"\ntrailertime must be >= 0.0");
+			exit(1);
+		}
+		trailersamps = (unsigned int) (trailertime * props.srate);
+	}
+	tapfile = fopen(argv[TAPFILE],"r");
+	if(!tapfile){
+		fprintf(stderr,"\nunable to open tapfile %s",argv[TAPFILE]);
+		usage();
+		sndcloseEx(ifd);
+		exit(1);
+	}
+
+	//read all lines, into taps array
+	vector<STAP> ltaps;
+	STAP thistap;
+	double this_time,this_amp,this_pan;
+	char line[255];
+
+	int linecount = 1;
+	this_time = 0.0; this_amp = 0.0;
+	bool read_error = false;
+	vector<STAP>::iterator I;
+
+	//RWD TODO: take a tap at time = 0.0 as spec of direct sound : amp and pan
+	
+	while(fgets(line,255,tapfile) != NULL){
+		this_pan = 0.0;
+		int got;
+		if(line[0] == '\n' || line[0]== ';'){	//allow comment lines; this does not check for leading white-space...
+			linecount++;
+			continue;
+		}
+		if((got = sscanf(line,"%lf%lf%lf",&this_time,&this_amp,&this_pan)) < 2){
+			fprintf(stderr,"\nerror in tapfile: line %d: insufficient parameters",linecount);
+			read_error = true;
+			break;
+		}
+		//time= 0 no good for a taptime!
+		if(this_time < 0.0){
+			fprintf(stderr,"\nerror in tapfile: line %d: bad time value",linecount);
+			read_error = true;
+			break;
+		}
+		if(this_time < sample_duration)	{ //non-zero time must be at least one sample on!			
+			fprintf(stderr,"\nWARNING: very small taptime %.4lf treated as zero - ignoring mix value",this_time);
+			this_time = 0.0;
+			mix_input = false;
+		}
+		if(current_time != 0.0 && this_time==current_time){
+			fprintf(stderr,"\nWARNING: duplicate time in line %d: ignored",linecount);
+			linecount++;
+			continue;
+		}		
+		if(this_time < current_time){
+			fprintf(stderr,"\nerror in tapfile: time out of sequence: line %d",linecount);
+			read_error = true;
+			break;
+		}
+		current_time = this_time;
+		thistap.time = this_time;
+		if(this_amp < 0.0 || this_amp > 1.0){
+			fprintf(stderr,"\nerror in tapfile: line %d: bad amplitude value",linecount);
+			read_error = true;
+			break;
+		}
+		
+		thistap.amp = this_amp;
+		if(got==3){
+			thistap.pan = this_pan;
+			is_stereo = true;
+			nchans = 2;
+		}
+		else
+			thistap.pan = 0.0;
+		ltaps.push_back(thistap);
+		linecount++;
+	}
+
+	if(read_error){
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	ntaps = ltaps.size();
+	if(ntaps==0){
+		printf("\nfile contains no data!");
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	
+	if(!mix_input)
+		printf("\nzero taptime used: taking input mix control from tapfile");
+#ifdef _DEBUG
+	printf("\nread %d taps",ntaps);
+	
+	for(I = ltaps.begin(); I != ltaps.end();I++)
+		printf("\n%.4lf\t%.4lf\t%.4lf",I->time,I->amp,I->pan);	//or whatever.....
+
+	printf("\nfeedback =  %.4f",feedback);
+	printf("\ndryfac = %.4f, wetfac = %.4f",dryfac,wetfac);
+
+#endif
+	if(is_stereo)
+		printf("\ncreating stereo outfile");
+	else
+		printf("\ncreating mono outfile");
+	
+	
+	l_taps = new deltap[ntaps];
+	
+	if(is_stereo){
+		r_taps = new deltap[ntaps];
+		//and create stereo delay-lines
+		int j;
+		for(j = 0,I = ltaps.begin(); I != ltaps.end();j++,I++){
+		//printf("\n%.4lf\t%.4lf\t%.4lf",I->time,I->amp,I->pan);	//or whatever.....
+			double l_gain,r_gain;
+			pancalc(I->pan,&l_gain,&r_gain);
+			l_taps[j].pos = r_taps[j].pos = I->time;
+			l_taps[j].val = I->amp * l_gain;
+			r_taps[j].val = I->amp * r_gain;
+		}
+	}
+	//else create a mono delayline
+	else {
+		int j;
+		for(j = 0,I = ltaps.begin(); I != ltaps.end();j++,I++){				
+			l_taps[j].pos = I->time;
+			l_taps[j].val = I->amp;
+		}
+	}
+
+
+	l_dline  = new tapdelay(l_taps,ntaps,props.srate);
+	if(!l_dline->create()){
+		fputs("\nunable to  create delay line - probably no memory",stderr);
+		delete [] l_taps;
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	max_lgain = l_dline->get_maxgain();
+	if(is_stereo){
+		r_dline = new tapdelay(r_taps,ntaps,props.srate);
+		if(!r_dline->create()){
+			fputs("\nunable to  create stereo delay line - probably no memory",stderr);
+			delete [] l_taps;
+			delete [] r_taps;
+			delete l_dline;
+			sndcloseEx(ifd);
+			exit(1);
+		}
+		max_rgain = r_dline->get_maxgain();
+	}
+	if(is_stereo)
+		max_gain = (float) min(max_lgain,max_rgain);
+	else
+		max_gain = (float) max_lgain;
+
+	/* recale max_gain by user param, eg:*/
+	max_gain *= (float) tapgain;
+
+
+	l_dline->set_gain((double)max_gain);
+	if(is_stereo)
+		r_dline->set_gain((double)max_gain);
+
+
+	props.chans = nchans;
+	peaks = (CHPEAK *) malloc(sizeof(CHPEAK) * nchans);
+	if(peaks==NULL){
+		fputs("No memory for PEAK data",stderr);
+		exit(1);
+	}
+	for(i=0; i < nchans; i++){
+		peaks[i].value = 0.0f;
+		peaks[i].position = 0;
+	}
+	//OK, now open outfile!
+
+	ofd = sndcreat_ex(argv[OUTFILE], -1,&props,SFILE_CDP,CDP_CREATE_NORMAL); 
+	//ofd = sndcreat_formatted(argv[OUTFILE],-1,outsamp_type,nchans,props.srate,CDP_CREATE_NORMAL);
+	if(ofd < 0){
+		printf("\nunable to open outfile %s",argv[OUTFILE]);
+		delete [] l_taps;
+		delete l_dline;
+		if(r_taps)
+			delete [] r_taps;
+		sndcloseEx(ifd);
+		exit(1);
+	}
+	printf("\n");
+
+	float l_op = 0.0f;		  
+	float r_op = 0.0f;
+	float l_out = 0.0f;
+	float r_out = 0.0f;
+	long outsamps = 0;
+	long step = (long) (0.25 * props.srate);
+	//RWD.10.98 BIG TODO: optimize all this!
+	sr = (double) props.srate;
+	if(is_stereo){
+		for(;;){
+			int rv;
+			float ip,l_ip,r_ip;										
+			if((rv = fgetfloatEx(&l_ip,ifd,0)) < 0){ 
+					fprintf(stderr,"\nerror reading infile data");
+					delete [] l_taps;
+					delete l_dline;
+					if(r_taps)
+						delete [] r_taps;
+					sndcloseEx(ifd);
+					exit(1);				
+			}
+			if(!rv) 
+				break;
+
+			ip = l_ip;
+			if(inchans==2){
+				//mix stereo ip to mono for delay purposers
+				if((rv = fgetfloatEx(&r_ip,ifd,0)) < 0){ 
+					fprintf(stderr,"\nerror reading infile data");
+					delete [] l_taps;
+					delete l_dline;
+					if(r_taps)
+						delete [] r_taps;
+					sndcloseEx(ifd);
+					exit(1);				
+				}
+				if(!rv) 
+				break;
+				ip = (l_ip + r_ip) * 0.5f;
+			}
+			
+		
+			
+			if(feedback != 0.0f) {
+				l_op = l_dline->tick(ip + (l_op * feedback));
+				r_op = r_dline->tick(ip + (r_op * feedback));
+			}
+			else{
+				l_op = l_dline->tick(ip);
+				r_op = r_dline->tick(ip);
+			}
+
+			//mix mono or stereo input with panned output:
+			if(mix_input){
+				l_out = l_op * wetfac + l_ip * dryfac;
+				r_out = r_op * wetfac;
+				if(inchans==2)
+					r_out  += r_ip * dryfac;
+				else
+					r_out += l_ip * dryfac;
+			}
+			else{
+				l_out = l_op;
+				r_out = r_op;
+			}
+			/* rescale */
+			l_out *=  max_gain;
+			r_out *=  max_gain;
+			if((float)fabs((double)l_out) > peaks[0].value) {
+				peaks[0].value = (float)fabs((double)l_out);
+				peaks[0].position = outsamps;
+			}
+			if((float)fabs((double)r_out) > peaks[1].value) {
+				peaks[1].value = (float)fabs((double)r_out);
+				peaks[1].position = outsamps;
+			}
+
+
+			if(fputfloatEx(&l_out,ofd) < 1){
+				fprintf(stderr,"\nerror writing to outfile: %s",sferrstr());
+				delete [] l_taps;
+				delete l_dline;
+				if(r_taps)
+					delete [] r_taps;
+				sndcloseEx(ifd);
+				exit(1);
+			}
+			if(fputfloatEx(&r_out,ofd) < 1){
+				fprintf(stderr,"\nerror writing to outfile: %s",sferrstr());
+				delete [] l_taps;
+				delete l_dline;
+				if(r_taps)
+					delete [] r_taps;
+				sndcloseEx(ifd);
+				exit(1);
+			}
+			outsamps++;
+			if((outsamps % step) == 0)
+				//inform(step,props.srate);
+				fprintf(stdout,"%.2f\r",(double)outsamps/sr);
+		}
+		for(i = 0; i < trailersamps;i++){
+			if(feedback != 0.0f){
+				l_op = l_dline->tick(feedback * l_op);
+				r_op = r_dline->tick(feedback * r_op);
+			}
+			else {
+				l_op = l_dline->tick(0.0f);
+				r_op = r_dline->tick(0.0f);
+			}
+			if(mix_input){			
+				l_out = l_op * wetfac;
+				r_out = r_op * wetfac;
+			}
+			else {
+				l_out = l_op;
+				r_out = r_op;
+			}
+			/* rescale */
+			l_out *=  max_gain;
+			r_out *=  max_gain;
+			if((float)fabs((double)l_out) > peaks[0].value) {
+				peaks[0].value = (float)fabs((double)l_out);
+				peaks[0].position = outsamps;
+			}
+			if((float)fabs((double)r_out) > peaks[1].value) {
+				peaks[1].value = (float)fabs((double)r_out);
+				peaks[1].position = outsamps;
+			}
+			if(fputfloatEx(&l_out,ofd) < 1){
+				fprintf(stderr,"\nerror writing to outfile: %s",sferrstr());
+				delete [] l_taps;
+				delete l_dline;
+				if(r_taps)
+					delete [] r_taps;
+				sndcloseEx(ifd);
+				exit(1);
+			}
+			if(fputfloatEx(&r_out,ofd) < 1){
+				fprintf(stderr,"\nerror writing to outfile: %s",sferrstr());
+				delete [] l_taps;
+				delete l_dline;
+				if(r_taps)
+					delete [] r_taps;
+				sndcloseEx(ifd);
+				exit(1);
+			}
+			outsamps++;
+			if((outsamps % step) == 0)
+				//inform(step,props.srate);
+				fprintf(stdout,"%.2f\r",(double)outsamps/sr);
+		}
+	}
+
+	else{
+		for(;;){
+			int rv;
+			float ip,l_ip,r_ip;			
+			if((rv = fgetfloatEx(&l_ip,ifd,0)) < 0){ 
+				fprintf(stderr,"\nerror reading infile data");
+				delete [] l_taps;
+				delete l_dline;
+				sndcloseEx(ifd);
+				exit(1);
+
+			}
+			if(!rv) 
+				break;
+			if(inchans==2){
+				//mix stereo ip to mono for delay purposers
+				if((rv = fgetfloatEx(&r_ip,ifd,0)) < 0){ 
+					fprintf(stderr,"\nerror reading infile data");
+					delete [] l_taps;
+					delete l_dline;
+					sndcloseEx(ifd);
+					exit(1);				
+				}
+				if(!rv) 
+				break;
+				ip = (l_ip + r_ip) * 0.5f;
+			}
+			else
+				ip = l_ip;
+			if(feedback != 0.0f)
+				l_op = l_dline->tick(ip + l_op * feedback);
+			else
+				l_op = l_dline->tick(ip);
+
+			if(mix_input)
+				l_out = l_op * wetfac + l_ip * dryfac;
+			else
+				l_out = l_op;
+			/* rescale */
+			l_out *=  max_gain;
+			if((float)fabs((double)l_out) > peaks[0].value) {
+				peaks[0].value = (float)fabs((double)l_out);
+				peaks[0].position = outsamps;
+			}
+			
+			if(fputfloatEx(&l_out,ofd) < 1){
+				fprintf(stderr,"\nerror writing to outfile");
+				delete [] l_taps;
+				delete l_dline;
+				sndcloseEx(ifd);
+				exit(1);
+			}
+			outsamps++;
+			if((outsamps % step) == 0)
+				//inform(step,props.srate);
+				fprintf(stdout,"%.2f\r",(double)outsamps/sr);
+		}
+		for(i = 0; i < trailersamps;i++){						
+			l_op = l_dline->tick(feedback * l_op);
+			if(mix_input)
+				l_out = l_op * wetfac;
+			else
+				l_out = l_op;
+			if((float)fabs((double)l_out) > peaks[0].value) {
+				peaks[0].value = (float)fabs((double)l_out);
+				peaks[0].position = outsamps;
+			}
+			if(fputfloatEx(&l_out,ofd) < 1){
+				fprintf(stderr,"\nerror writing to outfile");
+				delete [] l_taps;
+				delete l_dline;
+				sndcloseEx(ifd);
+				exit(1);
+			}
+			outsamps++;
+			if((outsamps % step) == 0)
+				//inform(step,props.srate);
+				fprintf(stdout,"%.2f\r",(double)outsamps/sr);
+		}
+	}
+	fprintf(stdout,"%.2f\n",(double)outsamps/sr);
+	printf("\nPEAK data:\n");
+	for(i=0;i < nchans;i++)
+		printf("Channel %u: %.4f at frame %d: %.4lf secs\n",i+1,
+			peaks[i].value,peaks[i].position, (double) peaks[i].position / (double)props.srate);
+	sndputpeaks(ofd,nchans,peaks);
+	if(sndcloseEx(ofd) < 0){
+		fprintf(stderr,"\nwarning: error closing outfile");
+
+	}
+	delete [] l_taps;
+	delete l_dline;
+	if(is_stereo) {
+		delete [] r_taps;
+		delete r_dline;
+	}
+	
+	sndcloseEx(ifd);
+	return 0;
+}
+
+	
+void usage(void)
+{
+	printf("\n*******  STEREO MULTI-TAPPED DELAY WITH PANNING v1.0 1998 : CDP Release 4 ******\n");
+	printf("\nusage:\n\ttapdelay [-f] infile outfile tapgain feedback mix taps.txt [trailtime]\n");
+	printf("\n\t-f\t\twrite floating-point outfile");
+	printf("\n\t\t\t(default: outfile has same format as infile)");
+	printf("\n\ttapgain\t\tgain factor applied to output from delayline");
+	printf("\n\t\t\t(tapgain > 0.0, typ 0.25)");
+	printf("\n\tfeedback\trange: -1 <= feeedback <= 1.0");
+	printf("\n\tmix\t\tproportion of source mixed with delay output.\n\t\t\trange: 0.0 <= mix < 1.0");
+	printf("\n\ttrailtime\textra time in seconds (beyond infile dur) ");
+	printf("\n\t\t\t  for delays to play out.");
+	printf("\n\tStereo inputs are mixed to mono at input to the delay line.\n");
+	printf("\n\ttaps.txt consists of list of taps in the format:");
+	printf("\n\ttime amp [pan]");
+	printf("\n\n\tAll values are floating-point, one tap per line.");
+	printf("\n\tTimes (seconds) must be increasing. Duplicate times are ignored.");
+	printf("\n\tA zero time (no delay) overrides the mix parameter,");
+	printf("\n\t  and determines the level and pan of the (mono) input.");	
+	printf("\n\tAmp values must be in the range 0.0 to 1.0");
+	printf("\n\tEmpty lines are permitted, as are lines starting with  a semni-colon,");	
+	printf("\n\t  which may be used for comments.");
+	printf("\n\tIf a pan value is used in any line, outfile will be stereo.");
+	printf("\n\tPan values are nominally in the range -1 to +1: 0 = centre (mono)");
+	printf("\n\t within which constant-power panning is used.");
+	printf("\n\tValues beyond these limits result in attenuation according to the");
+	printf("\n\t   inverse-square law, to suggest distance beyond the speakers.\n");
+	
+}
+
+//TW's constant-power pan routine, even does inv-square reductions beyond the speakers
+void pancalc(double position,double *leftgain,double *rightgain)
+{
+	int dirflag;
+	double temp;
+	double relpos;
+	double reldist, invsquare;
+
+	if(position < 0.0)
+		dirflag = SIGNAL_TO_LEFT;		/* signal on left */
+	else
+		dirflag = SIGNAL_TO_RIGHT;
+
+	if(position < 0) 
+		relpos = -position;
+	else 
+		relpos = position;
+	if(relpos <= 1.0){		/* between the speakers */
+		temp = 1.0 + (relpos * relpos);
+		reldist = ROOT2 / sqrt(temp);
+		temp = (position + 1.0) / 2.0;
+		*rightgain = temp * reldist;
+		*leftgain = (1.0 - temp ) * reldist;
+	} else {				/* outside the speakers */
+		temp = (relpos * relpos) + 1.0;			 /*RWD: NB same in both cases! */
+		reldist  = sqrt(temp) / ROOT2;   /* relative distance to source */
+		invsquare = 1.0 / (reldist * reldist);
+		if(dirflag == SIGNAL_TO_LEFT){
+			*leftgain = invsquare;
+			*rightgain = 0.0;
+		} else {   /* SIGNAL_TO_RIGHT */
+			*rightgain = invsquare;
+			*leftgain = 0;
+		}
+	}
+}

dev/externals/reverb/wavetable.cpp

@@ -0,0 +1,809 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+// Jan 2009  extended random lfo to support gaussian distribution
+
+
+//////////////////////////////////////////////////////////////////////
+
+#include <math.h>
+#include <memory.h>
+#include <time.h>
+#include "wavetable.h"
+#include <assert.h>
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+#ifndef PI
+#define PI (3.141592634)
+#endif
+
+#ifndef TWOPI
+#define TWOPI (2.0 * PI)
+#endif
+
+#define LAMBDA (10.0)
+
+
+/* rand stuff from Csound*/
+static const long RANDINT_MAX = 2147483647L;        // 2^^31 - 1 */
+static const long BIPOLAR = 0x4fffffff;
+static const long RIA = 16807;
+static const long RIM = 0x7fffffffL;
+static const double dv2_31 = 1.0 / 2147483648.0;
+static const long MAXLEN = 0x1000000L;
+static const long PHMASK = 0x0ffffffL;
+static const double FMAXLEN = (double)(0x1000000L);
+// for tpdf - bit of a hack, but seems close enough!
+static const long TPDFSHIFT = 0x18000000;
+// for gauss, trial and error!
+static const long BISHIFT = 0x16000000;
+
+
+
+fastlfo::fastlfo()
+{
+	param.freq = 0.0;
+	param.modrange = 0.0;
+	m_cur_WaveType = LFO_SINE;
+	curfreq = param.freq = 0.0;
+	curphase = 0.0;
+	incr = 0.0;
+	tf = 0;
+	kicvt = 0;
+	offset= 0.0;
+	b_sync = false;
+	tpdf = false;
+    gauss = false;
+    /* for csrand */
+    csseed = 12345;
+}
+
+double      fastlfo::csrand()  
+{ 
+    long rval;
+    rval =  randint31(csseed);
+    csseed = rval;
+    return (double) rval / (double) RANDINT_MAX;
+}
+
+
+long fastlfo::init(double srate, double normphase, long seedval, unsigned long ksmps)
+{
+	long seed;
+
+	m_srate = srate / ksmps;
+	twopiovrsr = TWOPI / m_srate;	
+	m_inv_srate = 1.0 / m_srate;
+	set_WaveType(LFO_SINE);
+	phs = (long)(MAXLEN * normphase);
+	curphase = TWOPI * normphase;
+	if(seedval == 0)
+		seed = (long) time(NULL);
+	else
+		seed = seedval;
+	rand = randint31(seed);
+	rand = randint31(rand);	
+    num1 = 0;
+    // get 2nd randon val for diff calc
+    // TODO: check this! why << 1 ?
+	rand = randint31(rand);
+	//num2 = (double)(rand<<1) * dv2_31;
+    num2 = (double)((long)((unsigned)rand<<1)-BIPOLAR) * dv2_31;
+
+	dfdmax = (num2 - num1) / FMAXLEN;
+	kicvt = (long)((double)MAXLEN / m_srate);
+	lastval = 0.0;
+	ksamps = ksmps;
+	if(ksamps==0)
+		ksamps++;
+	kcount = 1;
+	b_sync = false;
+	curphs = (long)(normphase * MAXLEN);
+	return seed;
+}
+
+
+// we keep current wavetype; no reset of random values ;
+/* phase ranged 0 -1 */
+void fastlfo::reset(double phase)
+{ 
+	curphase  = TWOPI * phase;
+	phs =  (long)(MAXLEN * phase);
+	curphs = phs; // assume phase offset = 0 here? 
+	lastval = 0.0;
+	b_sync = false;
+	kcount = 1;
+}
+
+
+
+#define OSC_WRAPPHASE  if(curphase >= TWOPI) curphase -= TWOPI;	\
+						if(curphase < 0.0) curphase += TWOPI
+// input range 0-1
+void fastlfo::sync_phase(double phase,double phaseoffset,double phaseincr)
+{
+	if(b_sync) {
+		curphase = phase * TWOPI;
+		offset = phaseoffset * TWOPI;
+		incr = phaseincr * TWOPI;
+		phs_incr = (long)(phaseincr*ksamps * MAXLEN);
+		phs = (long)(phase * MAXLEN);
+		phs &= PHMASK;
+		phs_offset = (long)(phaseoffset * MAXLEN);
+		curphs = phs + phs_offset;
+		curphs &= PHMASK;
+	}
+}
+
+inline double fastlfo::sinetick(void)
+{
+	double val;
+	double thisphase = curphase + offset;	
+	while(thisphase >= TWOPI)
+		thisphase -= TWOPI;
+	if(kcount==ksamps){
+		if(!b_sync){
+			val = sin(thisphase);
+			if(curfreq != param.freq){
+				curfreq = param.freq;
+				incr = twopiovrsr * param.freq;		
+			}
+		}
+		else
+			val = sin(thisphase);
+		curphase += incr;
+		//OSC_WRAPPHASE;
+		if(curphase >= TWOPI) 
+			curphase -= TWOPI;	
+		if(curphase < 0.0) 
+			curphase += TWOPI;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return lastval;
+}
+
+inline double fastlfo::tritick(void)
+{
+	double val;
+	double thisphase = curphase + offset;	
+	while(thisphase >= TWOPI)
+		thisphase -= TWOPI;
+	if(kcount==ksamps){
+		if(!b_sync){
+			if(curfreq != param.freq){
+				curfreq = param.freq;
+				incr = twopiovrsr * param.freq;		
+			}
+		}
+		if(thisphase <= PI)	{
+			val =   (4.0 * (thisphase * (1.0 / TWOPI) )) - 1.0;
+		}
+		else {
+			val =  3.0 - 4.0 * (thisphase * (1.0 / TWOPI) );
+		}
+		curphase += incr;
+		OSC_WRAPPHASE;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+	
+}
+
+inline double fastlfo::squaretick(void)
+{
+	double val;
+	double thisphase = curphase + offset;	
+	while(thisphase >= TWOPI)
+		thisphase -= TWOPI;
+	if(kcount==ksamps){
+		if(!b_sync){
+			if(curfreq != param.freq){
+				curfreq = param.freq;
+				incr = twopiovrsr * param.freq;		
+			}
+		}
+		if(thisphase <= PI)
+			val = 1.0;
+		else
+			val = -1;
+		curphase += incr;
+		OSC_WRAPPHASE;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+	
+}
+
+double fastlfo::sawuptick(void)
+{
+	double val;
+	double thisphase = curphase + offset;	
+	while(thisphase >= TWOPI)
+		thisphase -= TWOPI;
+
+	if(kcount==ksamps){
+		if(!b_sync){
+			if(curfreq != param.freq){
+				curfreq = param.freq;
+				incr = twopiovrsr * param.freq;		
+			}
+		}
+		val =  (2.0 * (thisphase * (1.0 / TWOPI) )) - 1.0;
+		curphase += incr;
+		OSC_WRAPPHASE;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+	
+}
+
+inline double fastlfo::sawdowntick(void)
+{
+	double val;
+	double thisphase = curphase + offset;	
+	while(thisphase >= TWOPI)
+		thisphase -= TWOPI;
+
+	if(kcount==ksamps){
+		if(!b_sync){
+			if(curfreq != param.freq){
+				curfreq = param.freq;
+				incr = twopiovrsr * param.freq;		
+			}
+		}
+		val =  1.0 - 2.0 * (thisphase * (1.0 / TWOPI) );
+		curphase += incr;
+		OSC_WRAPPHASE;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+
+// some redundant code with phs and curphs here, but keep both for now!
+
+
+inline double fastlfo::randomtick(void)
+{
+	double val;
+	double dval = 0.0;
+    int i;
+
+	if(kcount==ksamps){
+		if(b_sync){
+			val = num1 +(double)(curphs * dfdmax);
+			phs += phs_incr;
+			curphs += phs_incr;
+			if(curphs >= MAXLEN){
+				long r =  randint31(rand);
+				if(tpdf){
+					long rr = randint31(r);
+                    long ri;
+                    unsigned long sum = r + rr;
+                    ri = (long)(sum/2);                       
+                    r =  ri - TPDFSHIFT;
+				}
+                else if(gauss){
+                    int nrands = 16;
+                    long ri = 0;
+                    long rr[16];
+                    rr[0] = rand;
+                    for(i = 1; i < nrands; i++)
+                        rr[i] = randint31(rr[i-1]);
+                    for(i=0;i < nrands;i++)
+                        ri += rr[i]>>4;                        
+                    r = ri - BISHIFT;    
+                }
+                else if(biexp){
+                    long rr = r;
+                    double dbiexp = 0.0;
+                                    
+                    dbiexp = 2.0 * (double)rr * dv2_31;
+                    while(dbiexp==0.0 || dbiexp == 2.0){                    
+                        rr = randint31(rr);                                            
+                    }                
+                    if(dbiexp > 1.0)
+                        dval = -log(2.0-dbiexp) / LAMBDA;
+                    if(dbiexp < 1.0)
+                        dval = log(dbiexp) / LAMBDA;
+                    if( dval > 0.9)
+                        dval = 0.9;
+                    if(dval < -0.9)
+                        dval = -0.9;                                        
+                    r =   ((unsigned long)(dval/dv2_31) +BIPOLAR)>>1;                                       
+                }
+				curphs &= PHMASK;
+				rand = r;
+				num1 = num2;
+				num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;
+				dfdmax = (num2-num1) / FMAXLEN;
+			}
+			if(phs >= MAXLEN)
+				phs &= PHMASK;
+		}
+		else{           
+            val = num1 +(double)(phs * dfdmax);
+            phs +=(long) (param.freq * kicvt);
+            if(phs >= MAXLEN){
+                long r =  randint31(rand);
+                assert(r>0);
+                if(tpdf){
+                    long rr = randint31(r);
+                    long ri;
+                    unsigned long sum = r + rr;
+                    ri = (long)(sum/2);                   
+                    r =  ri - TPDFSHIFT;                     
+                }
+                else if(gauss){
+                    int nrands = 16;
+                    long ri = 0;
+                    long rr[16];
+                    rr[0] = rand;
+                    for(i = 1;i < nrands;i++)
+                        rr[i]= randint31(rr[i-1]);
+                    for(i=0;i < nrands;i++)
+                        ri += rr[i]>>4;                        
+                    r = ri - BISHIFT;                        
+                }
+                else if(biexp){
+                    long rr = r;
+                    double dbiexp = 0.0;
+                                    
+                    dbiexp = 2.0 * (double)rr * dv2_31;
+                    while(dbiexp==0.0 || dbiexp == 2.0){                    
+                        rr = randint31(rr);                                            
+                    }                
+                    if(dbiexp > 1.0)
+                        dval = -log(2.0-dbiexp) / LAMBDA;
+                    if(dbiexp < 1.0)
+                        dval = log(dbiexp) / LAMBDA;
+                    if( dval > 0.9)
+                        dval = 0.9;
+                    if(dval < -0.9)
+                        dval = -0.9;                                        
+                    r =   ((unsigned long)(dval/dv2_31) +BIPOLAR)>>1;
+                    /* assert(r>0);
+                      assert(r < RANDINT_MAX);
+                    */
+                }            
+                phs &= PHMASK;
+                rand = r;
+                num1 = num2;
+                num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;
+                /* NB makes uniform rand not really uniform, more a sort of rounded distr
+                 * so maybe not ideal for strict plain uniform distr - still while though!
+                */
+                dfdmax = (num2-num1) / FMAXLEN;                            
+            }
+		}
+        if(gauss)
+            val *= 1.85;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+inline double fastlfo::rand_tpdf_tick(void)
+{
+	double val;
+
+	if(kcount==ksamps){
+		if(b_sync){
+			val = num1 +(double)(curphs * dfdmax);
+			phs += phs_incr;
+			curphs += phs_incr;
+			if(curphs >= MAXLEN){
+				long r =  randint31(rand);
+				
+				long rr = randint31(r);
+                long ri;
+                long sum = r + rr;
+                ri = (long)(sum/2);                       
+                r =  ri - TPDFSHIFT;
+				                              
+				curphs &= PHMASK;
+				rand = r;
+				num1 = num2;
+				num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;
+				dfdmax = (num2-num1) / FMAXLEN;
+			}
+			if(phs >= MAXLEN)
+				phs &= PHMASK;
+		}
+		else{           
+            val = num1 +(double)(phs * dfdmax);
+            /* need this test for audiorate noise; but is it the right test for everywhere? */
+            /* time-varying noise rate? */
+            if(ksamps > 1)
+                phs +=(long) (param.freq * kicvt + 0.5);
+            else
+                phs = MAXLEN;
+            if(phs >= MAXLEN){
+                long r =  randint31(rand);                                
+                long rr = randint31(r);
+                long rsum;
+#ifdef _DEBUG
+                double d1,d2,dtpdf;
+                d1 = ((2.0 *(double)r)  / (double) RANDINT_MAX) - 1.0;
+                d2 = ((2.0 *(double)rr) / (double) RANDINT_MAX) - 1.0;
+                dtpdf = d1 + d2;
+                dtpdf *= 0.5;
+#endif
+                rand = rr;
+                /* bipolar value  =  2*r - 1  or 2* (r - 0.5)*/
+                long bipolard2 = 0x40000000;
+                r  = (r - bipolard2);
+                rr = (rr- bipolard2);
+                rsum = r + rr;
+                double dsum = (double)rsum  / (double) RANDINT_MAX ;                                                                                          
+                phs &= PHMASK; 
+                num1 = num2;
+                num2 = dsum;
+                dfdmax = (num2-num1) / FMAXLEN;                            
+            }
+		}        
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+inline double fastlfo::rand_gauss_tick(void)
+{
+	double val;
+    int i;
+
+	if(kcount==ksamps){
+		if(b_sync){
+			val = num1 +(double)(curphs * dfdmax);
+			phs += phs_incr;
+			curphs += phs_incr;
+			if(curphs >= MAXLEN){
+				long r =  randint31(rand);				                
+                int nrands = 16;
+                long ri = 0;
+                long rr[16];
+
+                rr[0] = rand;
+                for(i = 1; i < nrands; i++)
+                    rr[i] = randint31(rr[i-1]);
+                for(i=0;i < nrands;i++)
+                    ri += rr[i]>>4;                        
+                r = ri - BISHIFT;    
+                                
+				curphs &= PHMASK;
+				rand = r;
+				num1 = num2;
+				num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;
+				dfdmax = (num2-num1) / FMAXLEN;
+			}
+			if(phs >= MAXLEN)
+				phs &= PHMASK;
+		}
+		else{           
+            val = num1 +(double)(phs * dfdmax);
+            if(ksamps > 1)
+                phs +=(long) (param.freq * kicvt + 0.5);
+            else
+                phs = MAXLEN;
+            
+            if(phs >= MAXLEN){
+                long r =  randint31(rand);                                                
+                int nrands = 16;
+                long ri = 0;
+                long rr[16];
+
+                rr[0] = rand;
+                for(i = 1;i < nrands;i++)
+                    rr[i]= randint31(rr[i-1]);
+                for(i=0;i < nrands;i++)
+                    ri += rr[i]>>4;                        
+                r = ri - BISHIFT;                        
+                                         
+                phs &= PHMASK;
+                rand = r;
+                num1 = num2;
+                num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;                
+                dfdmax = (num2-num1) / FMAXLEN;                            
+            }
+		}
+        
+        val *= 1.85;
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+inline double fastlfo::rand_exp_tick(void)
+{
+	double val;
+	double dval = 0.0;
+
+	if(kcount==ksamps){
+		if(b_sync){
+			val = num1 +(double)(curphs * dfdmax);
+			phs += phs_incr;
+			curphs += phs_incr;
+			if(curphs >= MAXLEN){
+				long r =  randint31(rand);				                
+                long rr = r;
+                double dbiexp = 0.0;
+                                    
+                dbiexp = 2.0 * (double)rr * dv2_31;
+                while(dbiexp==0.0 || dbiexp == 2.0){                    
+                    rr = randint31(rr);                                            
+                }                
+                if(dbiexp > 1.0)
+                    dval = -log(2.0-dbiexp) / LAMBDA;
+                if(dbiexp < 1.0)
+                    dval = log(dbiexp) / LAMBDA;
+                if( dval > 0.9)
+                    dval = 0.9;
+                if(dval < -0.9)
+                    dval = -0.9;                                        
+                r =   ((unsigned long)(dval/dv2_31) +BIPOLAR)>>1;                                       
+                
+				curphs &= PHMASK;
+				rand = r;
+				num1 = num2;
+				num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;
+				dfdmax = (num2-num1) / FMAXLEN;
+			}
+			if(phs >= MAXLEN)
+				phs &= PHMASK;
+		}
+		else{           
+            val = num1 +(double)(phs * dfdmax);
+            if(ksamps > 1)
+                phs +=(long) (param.freq * kicvt + 0.5);
+            else
+                phs = MAXLEN;            
+            if(phs >= MAXLEN){
+                long r =  randint31(rand);                                
+                long rr = r;
+                double dbiexp = 0.0;
+                                    
+                dbiexp = 2.0 * (double)rr * dv2_31;
+                while(dbiexp==0.0 || dbiexp == 2.0){                    
+                    rr = randint31(rr);                                            
+                }                
+                if(dbiexp > 1.0)
+                    dval = -log(2.0-dbiexp) / LAMBDA;
+                if(dbiexp < 1.0)
+                    dval = log(dbiexp) / LAMBDA;
+                if( dval > 0.9)
+                    dval = 0.9;
+                if(dval < -0.9)
+                    dval = -0.9;                                        
+                r =   ((unsigned long)(dval/dv2_31) +BIPOLAR)>>1;
+                                                                  
+                phs &= PHMASK;
+                rand = r;
+                num1 = num2;
+                num2 = (double)((long)((unsigned)r<<1)-BIPOLAR) * dv2_31;
+                
+                dfdmax = (num2-num1) / FMAXLEN;                            
+            }
+		}        
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+
+inline double fastlfo::randhtick(void)
+{
+	double val;
+	
+	if(kcount==ksamps){
+		if(b_sync){
+			val = num1;
+			phs += phs_incr;
+			curphs += phs_incr;
+			if(curphs >= MAXLEN){
+				long r = randint31(rand);
+				if(tpdf){
+					long rr = randint31(r);
+                    long ri;
+                    unsigned long sum = r + rr;
+                    ri = (long)(sum/2);                       
+                    r =  ri - TPDFSHIFT;
+				}
+				rand = r;
+				num1 = (double)((long)((unsigned)r<<1) - BIPOLAR) * dv2_31;
+				curphs &= PHMASK;
+			}
+			if(phs >= MAXLEN)
+				phs &= PHMASK;
+		}
+		else {
+			val = num1;
+            if(ksamps > 1)
+                phs +=(long) (param.freq * kicvt + 0.5);
+            else
+                phs = MAXLEN;
+			
+			if(phs >= MAXLEN){
+				long r = randint31(rand);
+				if(tpdf){
+					long rr = randint31(r);
+                    long ri;
+                    unsigned long sum = r + rr;
+                    ri = (long)(sum/2);                       
+                    r =  ri - TPDFSHIFT;
+				}
+				phs &= PHMASK;
+				rand = r;
+				num1 = (double)((long)((unsigned)r<<1) - BIPOLAR) * dv2_31;
+			}
+		}
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+
+inline double fastlfo::randh_tpdf_tick(void)
+{
+	double val;
+	
+	if(kcount==ksamps){
+		if(b_sync){
+			val = num1;
+			phs += phs_incr;
+			curphs += phs_incr;
+			if(curphs >= MAXLEN){
+				long r = randint31(rand);
+				
+				long rr = randint31(r);
+                long ri;
+                unsigned long sum = r + rr;
+                ri = (long)(sum/2);                       
+                r =  ri - TPDFSHIFT;
+				
+				rand = r;
+				num1 = (double)((long)((unsigned)r<<1) - BIPOLAR) * dv2_31;
+				curphs &= PHMASK;
+			}
+			if(phs >= MAXLEN)
+				phs &= PHMASK;
+		}
+		else {
+			val = num1;			
+			if(ksamps > 1)
+                phs +=(long) (param.freq * kicvt + 0.5);
+            else
+                phs = MAXLEN;
+			if(phs >= MAXLEN){
+				long r = randint31(rand);
+				
+				long rr = randint31(r);
+                long ri;
+                unsigned long sum = r + rr;
+                ri = (long)(sum/2);                       
+                r =  ri - TPDFSHIFT;
+				
+				phs &= PHMASK;
+				rand = r;
+				num1 = (double)((long)((unsigned)r<<1) - BIPOLAR) * dv2_31;
+			}
+		}
+		lastval = val * param.modrange;
+		kcount = 1;
+	}
+	else
+		kcount++;
+	return 	lastval;
+}
+
+/* return positive value between 0 and RANDINT_MAX */
+inline long fastlfo::randint31(long seed31)
+{
+	unsigned long rilo, rihi;
+
+    rilo = RIA * (long)(seed31 & 0xFFFF);
+    rihi = RIA * (long)((unsigned long)seed31 >> 16);
+    rilo += (rihi & 0x7FFF) << 16;
+    if (rilo > (unsigned long) RIM) {
+      rilo &= RIM;
+      ++rilo;
+    }
+    rilo += rihi >> 15;
+    if (rilo > (unsigned long) RIM) {
+      rilo &= RIM;
+      ++rilo;
+    }
+    return (long)rilo;
+}
+
+bool fastlfo::set_WaveType(LfoWaveType type)		
+{	
+	switch(type){
+	case(LFO_SINE):
+            tf = &fastlfo::sinetick;
+		break;
+	case(LFO_TRIANGLE):
+		tf = &fastlfo::tritick;
+		break;
+	case(LFO_SQUARE):
+		tf = &fastlfo::squaretick;
+		break;
+	case(LFO_SAWUP):
+		tf = &fastlfo::sawuptick;
+		break;
+	case(LFO_SAWDOWN):
+		tf = &fastlfo::sawdowntick;
+		break;
+	case(LFO_RANDOM):
+		tf = &fastlfo::randomtick;
+		break;
+       case(LFO_RAND_TPDF):
+		tf = &fastlfo::rand_tpdf_tick;
+		break;
+       case(LFO_RAND_GAUSS):
+		tf = &fastlfo::rand_gauss_tick;
+		break;
+       case(LFO_RAND_EXP):
+		tf = &fastlfo::rand_exp_tick;
+		break;
+	case(LFO_RND_SH):
+		tf = &fastlfo::randhtick;
+		break;
+       case(LFO_RANDH_TPDF):
+		tf = &fastlfo::randh_tpdf_tick;
+		break;
+	default:
+		return false;
+	}
+	m_cur_WaveType = type;	
+	return true;
+}
+

dev/externals/reverb/wavetable.h

@@ -0,0 +1,126 @@
+/*
+ * Copyright (c) 1983-2013 Richard Dobson and Composers Desktop Project Ltd
+ * http://people.bath.ac.uk/masrwd
+ * http://www.composersdesktop.com
+ * This file is part of the CDP System.
+ * The CDP System is free software; you can redistribute it
+ * and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation; either 
+ * version 2.1 of the License, or (at your option) any later version. 
+ *
+ * The CDP System is distributed in the hope that it will be useful, 
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of 
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
+ * See the GNU Lesser General Public License for more details.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with the CDP System; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+ 
+
+// wavetable.h: interface for the fast_lfo class.
+// Feb 2009 added gaussian option to random
+//FWD 3.20: eliminate warning of unused private var (phsmask)
+
+//////////////////////////////////////////////////////////////////////
+
+
+#if !defined(AFX_WAVETABLE_H__2F961225_3FF4_11D2_96D4_444553540000__INCLUDED_)
+#define AFX_WAVETABLE_H__2F961225_3FF4_11D2_96D4_444553540000__INCLUDED_
+
+#if _MSC_VER >= 1000
+#pragma once
+#endif // _MSC_VER >= 1000
+
+typedef enum lfowavetype {LFO_SINE,LFO_TRIANGLE,LFO_SAWUP,LFO_SAWDOWN,LFO_SQUARE,LFO_RANDOM,LFO_RND_SH,
+LFO_RAND_TPDF,LFO_RAND_GAUSS,LFO_RAND_EXP, LFO_RANDH_TPDF,
+/*LFO_RANDH_GAUSS,LFO_RANDH_EXP, */
+LFO_NUM_WAVES } LfoWaveType;
+
+enum {PVS_LFO_FREQ_LORANGE, PVS_LFO_FREQ_HIRANGE};
+
+
+typedef struct lfoparam
+{
+	double freq;
+	double modrange;
+} LfoParam;
+
+
+
+class fastlfo {
+	typedef double (fastlfo:: *tickfunc)(void);
+public:
+
+	fastlfo();
+	virtual ~fastlfo() {}
+	long init(double srate, double normphase = 0.0,long seedval = 0,unsigned long ksmps = 1);
+	double		tick(void)							{ return (this->*tf)(); }
+	tickfunc	tf;
+	void		set_freq(double freq)				{ param.freq = freq;	}
+	double		get_freq(void)				const	{ return param.freq;	}
+	void		set_mod(double mod)					{ param.modrange = mod;	}
+	double		get_mod(void)				const	{ return param.modrange;}
+	void		reset(double phase);
+	bool		set_WaveType(LfoWaveType type);
+	void		set_tpdf()              		    { tpdf = true; gauss = biexp = false;}
+    void        set_white()                         {  tpdf = gauss = biexp = false;}
+    void        set_biexp()                         {  tpdf = gauss = false; biexp = true;}
+    void        set_flat()                          { tpdf = gauss = biexp = false;}
+    void        set_gaussian()                      { tpdf = biexp = false; gauss = true;}
+	void		sync_phase(double phase, double offset, double phaseincr);
+	void		set_sync(bool onoff)				{ b_sync = onoff;}
+    double      csrand();                   
+private:
+	double	sinetick(void);
+	double	tritick(void);
+	double	squaretick(void);
+	double	sawuptick(void);
+	double	sawdowntick(void);
+	double	randomtick(void);
+	double	randhtick(void);
+    double  rand_tpdf_tick(void);
+    double  rand_gauss_tick(void);
+    double  rand_exp_tick(void);
+    double  randh_tpdf_tick(void);
+ //   double  randh_gauss_tick(void);
+ //   double  randh_exp_tick(void);
+
+	long	randint31(long seed31);
+
+	double	m_srate, m_inv_srate;
+	double twopiovrsr;
+	double curfreq;
+	double curphase;
+	double incr;
+	LfoParam	param;
+	LfoWaveType  m_cur_WaveType;  
+	/* vars for random oscs, with a little help from Csound!*/
+	long		phs;
+	long		kicvt;
+//	long		phsmask;   // currently not used
+	long		rand;
+	double		num1,num2;
+	double		dfdmax;
+	/* for krate output */
+	unsigned long		ksamps;
+	unsigned long		kcount;
+	double		lastval;
+	bool	tpdf;
+    bool gauss;
+    bool biexp;
+	// for lfo sync
+	bool b_sync;
+	double		offset;
+	// for random oscs
+	long curphs;
+	long phs_incr;
+	long phs_offset;
+    /* for csrand */
+    long csseed;
+};
+
+#endif // !defined(AFX_WAVETABLE_H__2F961225_3FF4_11D2_96D4_444553540000__INCLUDED_)
+
+