initial commit

dev/externals/CMakeLists.txt

@@ -0,0 +1,8 @@
+add_subdirectory(portsf)
+add_subdirectory(fastconv)
+add_subdirectory(reverb)
+add_subdirectory(mctools)
+
+##if(USE_LOCAL_PORTAUDIO)
+  add_subdirectory(paprogs)
+##endif()

dev/externals/fastconv/CMakeLists.txt

@@ -0,0 +1,21 @@
+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")
+  else()
+    set(CMAKE_C_FLAGS "-O3 -Wall -Dlinux -Dunix -fomit-frame-pointer -funroll-loops")
+  endif()
+endif()
+
+link_directories(../lib)
+
+include_directories(../include)
+
+add_executable(fconv fconv.cpp genrespframe2.cpp mxfftd.c)
+
+target_link_libraries(fconv portsf sfsys ${EXTRA_LIBRARIES})
+
+my_install(fconv)
+

dev/externals/fastconv/fconv-fftw.cpp

@@ -0,0 +1,768 @@
+/*
+ * 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
+ *
+ */
+
+/* fastconv.cpp: */
+/*  self-contained prog for fast convolution (reverb, FIR filtering, bformat etc) 
+ *  
+*/
+
+/* um, there is nothing C++ish about this code apart from use of new, delete! */
+
+/*TODO: control wet/dry mix with something... */
+/* auto-rescale arbitrary impulse responses? */
+/* Feb 2013: rebuilt with updated portsf */
+/* August 2013 epanded usage message re text file. */
+extern "C"
+{
+
+#include <portsf.h>
+
+}
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
+#include <sys/types.h>
+#include <sys/timeb.h>
+
+
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+
+#ifdef unix
+#include <ctype.h>
+int stricmp(const char *a, const char *b);
+int strnicmp(const char *a, const char *b, const int length);
+#endif
+
+
+#include <rfftw.h>
+
+
+#ifndef WAVE_FORMAT_IEEE_FLOAT
+#define WAVE_FORMAT_IEEE_FLOAT (0x0003)
+#endif
+
+//#define FFTTEST
+
+/* convert from mono text impulse file */
+long genimpframe1(double *insbuf, double*** outbuf, long npoints, double scalefac);
+/* convert from multichan soundfile */
+long genimpframe2(int ifd,double*** outframe, double* rms,long imchans,double scalefac);
+void mc_split(double* inbuf,double** out,long insize,long chans);
+void mc_interl(double** in,double* out,long insize,long chans);
+void complexmult(double *frame,const double *impulse,long length);
+
+#ifdef FFTTEST
+extern "C"{
+void fft_(double *, double *,int,int,int,int);
+void fftmx(double *,double *,int,int,int,int,int,int *,double *,double *,double *,double *,int *,int[]);
+void reals_(double *,double *,int,int);
+}
+#endif
+
+#define DEFAULT_BUFLEN (32768)
+
+enum {ARG_NAME,ARG_INFILE,ARG_IMPFILE,ARG_OUTFILE,ARG_NARGS};
+
+void usage(const char *progname)
+{
+
+    printf("\n\n%s v1.2 RWD,CDP July 2010,2013",progname);
+	printf("\nMulti-channel FFT-based convolver");
+
+	printf("\nUsage: \n    %s [-aX][-f] infile impulsefile outfile [dry]\n"		               
+           "   -aX        : scale output amplitude by X\n"
+           "   -f         : write output as floats (no clipping)\n"
+           "  infile      : input soundfile to be processed.\n"
+           "  impulsefile : m/c soundfile or mono text file containing impulse response,\n"
+           "                  e.g. reverb or FIR filter.\n"
+           "                Text file name must have extension .txt.\n"
+		   "                  File must contain 1 column of floating point values,\n"
+		   "                  typically in the range -1.0 to 1.0.\n"
+           "              Supported channel combinations:\n"
+           "               (a) mono infile, N-channel impulsefile;\n"
+           "               (b) channels are the same;\n"
+           "               (c) mono impulsefile, N-channel infile.\n"
+           "  [dry]       :  set dry/wet mix (e.g. for reverb)\n"
+           "                 Range: 0.0 - 1.0,  default = 0.0\n"
+           "                 (uses sin/cos law for constant power mix)\n"
+           "Note: some recorded reverb impulses effectively include the direct signal.\n"
+           "In such cases  <dry>  need not be used\n"
+           "Where impulsefile is filter response (FIR), optimum length is power-of-two - 1.\n",progname
+		   );		   
+}
+
+
+double
+timer()
+{
+	struct timeb now;
+	double secs, ticks;	
+	ftime(&now);
+	ticks = (double)now.millitm/(1000.0);
+	secs = (double) now.time;
+
+	return secs + ticks;
+}
+
+
+void			
+stopwatch(int flag) 
+{
+	static double start;		    
+	long mins=0,hours=0;
+	double end, secs=0;
+
+	if(flag)
+		start = timer();	
+	else 
+	{
+		end    = timer(); 
+		secs   = end - start;
+		mins   = (long)(secs/60.0);
+		secs  -= mins*60.0; 
+		hours  = mins/60;
+		mins  -= hours*60;
+
+		fprintf(stderr,"\nElapsed time: ");
+		if(hours > 0)
+			fprintf(stderr,"%ld h ", hours);
+		if(mins > 0)
+			fprintf(stderr,"%ld m ", mins);
+		fprintf(stderr,"%2.3lf s\n\n",secs);
+	}
+}
+
+/* how do we want this to work?
+
+  (1)  multi-chan impulse:  EITHER: mono infile  OR infile with same chan count
+  (2)  mono impulse:  expanded to infile chancount
+     Therefore: need to open both infiles before deciding action
+*/
+#define PI2 (1.570796327)
+
+int main(int argc,char **argv)
+{
+
+	long fftlen = 1,imlen = 0,chans = 0,srate;
+	long buflen = 0;
+	double scalefac = 1.0f;
+	double ampfac = 1.0;
+    double Ninv = 1.0;
+	int insndfile = -1,inimpfile = -1,outsndfile = -1; 
+    int dorms = 0;
+    int nameoffset = 0;
+	PSF_PROPS inprops,outprops, impulseprops;
+    psf_format format;
+    PSF_CHPEAK  *fpeaks = NULL;
+	MYLONG peaktime;
+	int i,jj,minheader = 0,do_timer = 1;
+    int writefloats = 0;
+    long framesneeded;
+	double oneovrsr;
+    double *inmonobuf = 0;
+	double *insbuf=0,*outsbuf = 0;
+    double **inbuf = 0, **outbuf = 0;    
+	double **imbuf = 0;
+	double **overlapbuf = 0;
+    double rms = 0.0;
+    double dry = 0.0;
+    double wet = 1.0;
+#ifdef FFTTEST
+    double *anal = 0;
+#endif
+    rfftwnd_plan* forward_plan = 0;
+    rfftwnd_plan* inverse_plan = 0;
+	
+    if(argv[0][0]=='.' && argv[0][1]=='/'){
+        nameoffset  = 2;
+    }
+	
+    if(argc==2){
+        if(strcmp(argv[1],"--version")==0){
+            printf("1.2.0.\n");
+            return 0;
+        }
+    }
+
+	if(psf_init()){
+		puts("unable to start portsf\n");
+		return 1;
+	}
+
+
+	if(argc < ARG_NARGS){		 
+		usage(argv[0]+nameoffset);
+		return(1);
+	}
+
+    
+
+    while(argc > 1 && argv[1][0]=='-'){				
+		switch(argv[1][1]){
+        case 'a':
+            scalefac =  atof(&(argv[1][2]));
+            if(scalefac <=0.0){
+                printf("Error: scalefac must be positive!\n");
+                return 1;
+            }
+            break;        
+        case 'f':
+            writefloats = 1;
+            break;
+        default:
+            break;
+        }
+        argv++;
+        argc--;
+    }		
+	/* 2 legal possibilities: infile and outfile, or -I used with infile only */
+	if(argc< ARG_NARGS){
+		fprintf(stderr,"\nInsufficient arguments.");
+		usage(argv[0]+nameoffset);
+		return(1);
+	}
+    if(argc==ARG_NARGS+1){
+        double dryarg;
+        dryarg = atof(argv[ARG_NARGS]);
+        if(dryarg < 0.0 || dryarg > 1.0){
+            printf("dry value out of range.\n");
+            return 0;
+        }
+        if(dryarg==1.0){
+            printf("Warning: dry=1 copies input!\n");
+            wet = 0.0;
+        }
+        else{
+            dry = sin(PI2 * dryarg);
+            wet = cos(PI2 * dryarg);
+        }
+    }
+
+		
+
+/* TODO:  where -F[file] is combined with -A, it signifies create analysis file
+           compatible with impulse file (e.g 50% overlap, etc) */
+
+	
+	/* open infile, check props */	
+		
+	if((insndfile = psf_sndOpen(argv[ARG_INFILE],&inprops, 0)) < 0){
+		fprintf(stderr,"\nUnable to open input soundfile %s",argv[1]);
+        delete []imbuf;
+		return 1;
+	}
+    srate = inprops.srate;
+	if(srate <=0){
+		fprintf(stderr,"\nBad srate found: corrupted file?\n");
+        delete []imbuf;
+		return 1;
+	}
+	chans = inprops.chans;
+    framesneeded = psf_sndSize(insndfile);
+    if(framesneeded <= 0){
+        printf("Error in input file - no data!\n");
+        psf_sndClose(insndfile);
+        return 1;
+    }
+    /* open impulse file */
+    /* check for soundfile */
+      
+    format = psf_getFormatExt(argv[ARG_IMPFILE]);
+    if(format==PSF_FMT_UNKNOWN){  /* must be a text file */
+        FILE *fp = 0;
+        char tmp[80];
+        char* dot;
+        int size = 0,got = 0;
+        int read = 0;
+
+        dot = strrchr(argv[ARG_IMPFILE],'.');
+        if(dot==NULL || stricmp(dot,".txt") != 0){
+            fprintf(stderr,"Error: impulse text file must have .txt extension.\n");
+            return 1;
+        }
+        fp = fopen(argv[ARG_IMPFILE],"r");
+        if(fp==0){
+            printf("Cannot open impulse text file %s\n.",argv[ARG_IMPFILE]);
+            return 1;
+        }
+        /* count lines! */
+        while(fgets(tmp,80,fp) != NULL)
+            size++;
+        if(size==0){
+            printf("Impulse textfile %s has no data!.\n",argv[ARG_IMPFILE]);
+            return 1;
+        }
+        rewind(fp);
+            
+        inmonobuf = new double[size+1];
+        for(i=0;i < size;i++)  {
+            read = fscanf(fp,"%lf",&inmonobuf[i]);
+            if(read==0){
+                i--;
+                continue;
+            }
+            if(read==EOF)
+                break;
+            got++;
+        }
+        imlen = got;
+        impulseprops.chans = 1;
+        fclose(fp);
+    }
+    else{
+        if((inimpfile = psf_sndOpen(argv[ARG_IMPFILE],&impulseprops, 0))< 0){
+            fprintf(stderr,"\nUnable to open impulse file %s",argv[ARG_IMPFILE]);       
+            return 0;
+        }
+        //printf("impulse file sr = %d\n",impulseprops.srate);
+        if(srate != impulseprops.srate){
+            printf("Error: files must have same sample rate");
+            delete []imbuf;
+            return 1;
+        }
+        long  len = psf_sndSize(inimpfile);
+        if(len <= 0){
+            printf("Error in impulse soundfile - no data!\n");
+            return 1;
+        }
+        if(impulseprops.chans > 1){
+            if(! (chans == 1 || chans == impulseprops.chans)){
+                fprintf(stderr,"\nChannel mismatch.\n    Infile must be mono or have same channels as irfile.\n");
+                return(1);
+            }
+            chans = impulseprops.chans;
+        }    
+    }
+
+    imbuf = new double*[chans];
+
+    /* if impulse file is mono, we will need to copy data to the other buffers */
+    // allocate  impulseprops.chans buffers ; may need more
+
+    if(inimpfile >=0){
+	    if((imlen = genimpframe2(inimpfile,&imbuf,&rms,impulseprops.chans, scalefac)) == 0){
+	    	fprintf(stderr,"Error reading impulse file\n");
+	    	return(1);
+	    }
+    }
+    
+    else if(imlen){        
+        genimpframe1(inmonobuf,&imbuf,imlen,scalefac);       
+    }
+    printf("mean rms level = %.4lf (%.2lf dB)\n",rms, 20.0 * log10(rms));
+    
+
+    framesneeded +=  imlen;        /* can  close outfile, when this length reached */
+
+    while(fftlen < imlen*2)			/* convolution style - double-length */
+		fftlen <<= 1;
+    double norm = sqrt(2.0);       /* relative: rms of 0dBFS sine is 0.707 = 1/root2 */
+    Ninv = 1.0 / sqrt(imlen*2);
+    Ninv *= norm;
+    Ninv /= imlen;
+    // take simple avg of  rms for adjustment factor.
+    // may not adequately represent some responses, e.g. with strong attack/earlies, soft reverb tail */
+    double rmsdif2 =  (-20.0 * log10(rms)) * 0.5;
+    double rmsadjust = pow(10, rmsdif2/20.0);
+#ifdef _DEBUG
+    printf("rescaling factor = %.6e\n",Ninv);   
+    printf("rmsadjust = %.4lf\n",rmsadjust);
+#endif
+   Ninv /= rmsadjust; 
+
+
+    /* copy buffers if required */
+    for(i = impulseprops.chans;i < chans;i++){
+        imbuf[i] = new double[fftlen+2];
+        memcpy(imbuf[i],imbuf[0],(fftlen+2)* sizeof(double));
+    }         
+    oneovrsr = 1.0 / (double) srate;
+    
+    /*make sure buflen is at least fftlen */
+//	if(fftlen > buflen)
+//		buflen = fftlen;
+
+    
+	/* main i/o buffers */
+	
+	insbuf     = new double [fftlen  * chans];
+    outsbuf    = new double [fftlen  * chans];
+	inbuf      = new double*[chans];
+    outbuf     = new double*[chans];  /* overlap-add buffers */
+    overlapbuf = new double*[chans];
+    
+    for(i=0;i < chans;i++){
+        inbuf[i]  = new double[fftlen+2];   // main in-place buf 
+        outbuf[i] = new double[fftlen+2];        
+	    /* channel buffers */	
+        overlapbuf[i] = new double[fftlen/2];
+        memset(overlapbuf[i],0,sizeof(double)*(fftlen/2));        
+    }
+#ifdef FFTTEST
+    anal = new double[fftlen+2];
+#endif
+    
+    
+    forward_plan = new rfftwnd_plan[chans];
+    inverse_plan = new rfftwnd_plan[chans];
+    int insize = (int) fftlen;
+    for(i=0;i < chans;i++){		    
+	    memset(inbuf[i], 0, sizeof(double) * (fftlen+2));
+        memset(outbuf[i],0,sizeof(double) * (fftlen+2));
+        
+        forward_plan[i] = rfftwnd_create_plan_specific(1,&insize, 
+		FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE,
+		inbuf[i],1,NULL,1);
+        inverse_plan[i] = rfftwnd_create_plan_specific(1,&insize, 
+		FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE,
+		inbuf[i],1,NULL,1);
+    }
+
+    
+	/* use generic init function */
+/*bool phasevocoder::init(long outsrate,long fftlen,long winlen,long decfac,float scalefac,
+						pvoc_scaletype stype,pvoc_windowtype wtype,pvocmode mode)
+*/
+	
+	/*create output sfile */
+	psf_stype smptype;
+	psf_format outformat;
+	/* will it be aiff, wav, etc? */
+	outformat = psf_getFormatExt(argv[ARG_OUTFILE]);
+	if(outformat==PSF_FMT_UNKNOWN){
+		fprintf(stderr,"\nOutfile name has unrecognized extension.");
+        delete []imbuf;
+		return(1);
+	}
+		
+		
+	smptype = inprops.samptype;
+    if(writefloats)
+        smptype= PSF_SAMP_IEEE_FLOAT;
+
+	/*the one gotcha: if output is floats and format is aiff, change to aifc */
+	if(smptype==PSF_SAMP_IEEE_FLOAT && outformat==PSF_AIFF){
+		fprintf(stderr,"Warning: AIFF output written as AIFC for float samples\n");
+		outformat = PSF_AIFC;
+	}
+
+	outprops          = inprops;
+	outprops.chans    = chans;
+	outprops.srate    = srate;
+	outprops.format   = outformat;
+	outprops.samptype = smptype;
+	outprops.chformat = STDWAVE;
+    /* if irfile is BFormat, need to set outfile fmt likewise */
+    if(impulseprops.chformat==MC_BFMT)
+        outprops.chformat = MC_BFMT;
+    /* I suppose people will want automatic decoding too...grrr */
+
+    fpeaks = new PSF_CHPEAK[chans];
+	if(fpeaks==NULL){
+		puts("no memory for fpeak data buffer\n");
+		return 1;
+	}
+
+	/*last arg: no clipping of f/p samples: we use PEAK chunk */
+	if((outsndfile = psf_sndCreate(argv[ARG_OUTFILE],&outprops,0,minheader,PSF_CREATE_RDWR)) <0 ){
+		fprintf(stderr,"\nUnable to open outfile %s",argv[ARG_OUTFILE]);
+        delete []imbuf;
+		return(1);
+	}
+				
+
+	
+	printf("\n");
+	
+	
+	stopwatch(1);
+
+	long written,thisblock,framesread;
+    long frameswritten = 0;
+	double intime= 0.0,outtime = 0.0;
+    int last = 0;
+				
+	while((framesread = psf_sndReadDoubleFrames(insndfile,insbuf,imlen)) > 0){        
+		written = thisblock =  0;
+
+		for(i = framesread * inprops.chans; i< fftlen * inprops.chans; i++)
+			insbuf[i] = 0.0f;
+        framesread = imlen;
+		memset(inbuf[0],0,(fftlen+2) * sizeof(double));
+        if(chans == inprops.chans)  {
+            /* must clean buffers! */
+		    for(i=0;i < chans;i++)
+                memset(inbuf[i],0,(fftlen+2) * sizeof(double));
+            mc_split(insbuf,inbuf,imlen * chans, chans);
+
+        }
+        else{
+            for(i=0;i < chans;i++) {
+                memset(inbuf[i],0,(fftlen+2) * sizeof(double));
+                memcpy(inbuf[i],insbuf,imlen * sizeof(double));
+                memset(outbuf[i],0,sizeof(double) * (fftlen+2));
+
+            }
+        }
+        if(impulseprops.chans==1){
+            for(jj = 0; jj < chans;jj++){
+#ifdef FFTTEST
+                int zz;
+                memcpy(anal,inbuf[jj],(fftlen+2) * sizeof(double));
+                fft_(anal,anal+1,1,fftlen/2,1,-2);
+	            reals_(anal,anal+1,fftlen/2,-2);
+                for(zz=0;zz < fftlen+2;zz++)
+                    anal[zz] *= 0.001;
+#endif
+                rfftwnd_one_real_to_complex(forward_plan[jj],inbuf[jj],NULL);			    				
+                complexmult(inbuf[jj],imbuf[0],fftlen+2);                
+			    rfftwnd_one_complex_to_real(inverse_plan[jj],(fftw_complex * )inbuf[jj],NULL);				
+            }           
+        }
+        else{
+            for(jj = 0; jj < chans;jj++){
+				rfftwnd_one_real_to_complex(forward_plan[jj],inbuf[jj],NULL);                   				
+                complexmult(inbuf[jj],imbuf[jj],fftlen+2);                
+				rfftwnd_one_complex_to_real(inverse_plan[jj],(fftw_complex * )inbuf[jj],NULL);    				    			    
+            }            
+        }
+        
+        /* overlap-add  for each channel */
+        /* TODO: verify  use of imlen here - should it be fftlen/2 -1 ? */
+        for(jj=0;jj < chans;jj++){
+            for(i=0;i < imlen;i++) {
+                outbuf[jj][i] = inbuf[jj][i] + overlapbuf[jj][i];
+                overlapbuf[jj][i] = inbuf[jj][i+imlen];
+            }
+        }
+		mc_interl(outbuf,outsbuf,imlen,chans);
+        if(inprops.chans == chans){
+            for(i=0;i < framesread; i++) {
+                for(jj=0;jj < chans;jj++){
+                    long  outindex = i*chans + jj;            
+                    outsbuf[outindex] *= Ninv;
+                    outsbuf[outindex] *= wet;
+                    outsbuf[outindex] += dry * insbuf[outindex];
+                }
+            }
+        }
+        /* elso mono input */
+        else {
+            for(i=0;i < framesread; i++) {
+                for(jj=0;jj < chans;jj++){
+                    long  outindex = i*chans + jj;
+                    double inval = dry *  insbuf[i];
+                    outsbuf[outindex] *= Ninv;
+                    outsbuf[outindex] *= wet;
+                    outsbuf[outindex] += inval;
+                }
+            }
+        }
+
+        if((written = psf_sndWriteDoubleFrames(outsndfile,outsbuf,framesread)) != framesread){
+		    fprintf(stderr,"\nerror writing to outfile");
+		    return(1);		               
+        }
+        frameswritten += written;
+
+		if(do_timer){
+			intime += (double)framesread * oneovrsr;
+			printf("Input time: %.2lf\r",intime);
+		}    
+    }
+    /* complete tail */
+
+    if(frameswritten < framesneeded){
+        // TODO: imlen? see above
+        mc_interl(overlapbuf,outsbuf,imlen,chans);
+        long towrite = framesneeded - frameswritten; 
+        for(i=0;i < towrite * chans; i++) {
+            outsbuf[i] *= Ninv;
+            outsbuf[i] *= wet;
+        }
+        if((written = psf_sndWriteDoubleFrames(outsndfile,outsbuf,towrite)) != towrite){
+	        fprintf(stderr,"\nerror writing to outfile");
+	        return(1);
+	    }
+    }
+
+	stopwatch(0);
+    if(psf_sndReadPeaks( outsndfile,fpeaks,&peaktime)){
+        double peakmax = 0.0;
+		printf("PEAK values:\n");
+		for(i=0; i < chans; i++)	{
+            peakmax = fpeaks[i].val > peakmax ? fpeaks[i].val : peakmax;
+			if(fpeaks[i].val < 0.0001f)
+				printf("CH %d:\t%e (%.2lf dB)\tat frame %10lu :\t%.4lf secs\n",i+1,
+				fpeaks[i].val,20.0*log10(fpeaks[i].val),fpeaks[i].pos,(double) (fpeaks[i].pos) / (double)srate);
+			else
+				printf("CH %d:\t%.4lf (%.2lf dB)\tat frame %10lu :\t%.4lf secs\n",i+1,
+				fpeaks[i].val,20.0 * log10(fpeaks[i].val),fpeaks[i].pos,(double) (fpeaks[i].pos) / (double)srate);		
+		}
+        if(peakmax > 1.0)
+            printf("Overflows!  Rescale by %.10lf\n",1.0 / peakmax);
+	}
+
+	if(insndfile >=0)
+		psf_sndClose(insndfile);
+	if(outsndfile >=0)
+		psf_sndClose(outsndfile);
+    if(inimpfile >=0)
+	    psf_sndClose(inimpfile);
+	psf_finish();
+    delete [] fpeaks;
+	if(insbuf)
+		delete [] insbuf;
+	if(outsbuf)
+		delete [] outsbuf;
+    
+
+    for(i=0;i < chans;i++){
+        if(inbuf){       
+		    delete [] inbuf[i];
+        }           
+        if(outbuf){
+		    delete [] outbuf[i];
+        }       	
+        if(imbuf){        
+	        delete [] imbuf[i];
+        }
+        if(overlapbuf)
+            delete [] overlapbuf[i];
+        if(forward_plan)
+            rfftwnd_destroy_plan(forward_plan[i]);
+        if(inverse_plan)
+            rfftwnd_destroy_plan(inverse_plan[i]);
+    }
+    delete [] outbuf;
+    delete [] inbuf;
+    delete [] imbuf;
+    delete [] overlapbuf;
+	return 0;
+}
+
+
+// insize is raw samplecount,buflen is insize/chans
+void mc_split(double* inbuf,double** out,long insize,long chans)
+{
+    long i,j,buflen = insize/chans;
+    double* pinbuf;
+
+    
+    for(j=0;j < chans;j++){
+        pinbuf = inbuf+j;
+        for(i=0;i < buflen;i++){
+            out[j][i] = *pinbuf;
+            pinbuf += chans;
+        }
+    }
+}
+
+
+/* insize is m/c frame count */
+void mc_interl(double** in,double* out,long insize,long chans)
+{
+    long i,j;
+    double* poutbuf;
+
+    for(j = 0;j < chans;j++){
+        poutbuf = out+j;
+        for(i=0;i < insize;i++){
+            *poutbuf = in[j][i];
+            poutbuf += chans;
+        }
+    }
+}
+
+/* OR:  apply scalefac to impulse responses */
+void complexmult(double *frame,const double *impulse,long length)
+{
+	double re,im;
+	
+	int i,j;
+	
+
+	for(i=0,j = 1;i < length;i+=2,j+=2){
+		re = frame[i] * impulse[i] - frame[j] * impulse[j];
+		im = frame[i] * impulse[j] + frame[j]* impulse[i];
+		frame[i] = re;
+		frame[j] = im;
+	}
+}
+
+#ifdef unix
+int stricmp(const char *a, const char *b)
+{
+	while(*a != '\0' && *b != '\0') {
+		int ca = islower(*a) ? toupper(*a) : *a;
+		int cb = islower(*b) ? toupper(*b) : *b;
+        
+		if(ca < cb)
+			return -1;
+		if(ca > cb)
+			return 1;
+        
+		a++;
+		b++;
+	}
+	if(*a == '\0' && *b == '\0')
+		return 0;
+	if(*a != '\0')
+		return 1;
+	return -1;
+}
+
+int
+strnicmp(const char *a, const char *b, const int length)
+{
+	int len = length;
+    
+	while(*a != '\0' && *b != '\0') {
+		int ca = islower(*a) ? toupper(*a) : *a;
+		int cb = islower(*b) ? toupper(*b) : *b;
+        
+		if(len-- < 1)
+			return 0;
+        
+		if(ca < cb)
+			return -1;
+		if(ca > cb)
+			return 1;
+        
+		a++;
+		b++;
+	}
+	if(*a == '\0' && *b == '\0')
+		return 0;
+	if(*a != '\0')
+		return 1;
+	return -1;
+}
+#endif
+
+
+
+
+
+

dev/externals/fastconv/fconv.cpp

@@ -0,0 +1,752 @@
+/*
+ * 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
+ *
+ */
+
+/* fastconv.cpp: */
+/*  self-contained prog for fast convolution (reverb, FIR filtering, bformat etc) 
+ *  
+*/
+
+/* um, there is nothing C++ish about this code apart from use of new, delete! */
+
+/*TODO: control wet/dry mix with something... */
+/* auto-rescale arbitrary impulse responses? */
+/* Feb 2013: rebuilt with updated portsf */
+/* August 2013 epanded usage message re text file. */
+extern "C"
+{
+
+#include <portsf.h>
+
+}
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
+#include <sys/types.h>
+#include <sys/timeb.h>
+
+
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+
+#ifdef unix
+#include <ctype.h>
+int stricmp(const char *a, const char *b);
+int strnicmp(const char *a, const char *b, const int length);
+#endif
+
+#ifndef WAVE_FORMAT_IEEE_FLOAT
+#define WAVE_FORMAT_IEEE_FLOAT (0x0003)
+#endif
+
+/* convert from mono text impulse file */
+int genimpframe1(double *insbuf, double*** outbuf, int npoints, double scalefac);
+/* convert from multichan soundfile */
+int genimpframe2(int ifd,double*** outframe, double* rms,int imchans,double scalefac);
+void mc_split(double* inbuf,double** out,int insize,int chans);
+void mc_interl(double** in,double* out,int insize,int chans);
+void complexmult(double *frame,const double *impulse,int length);
+
+
+extern "C"{
+void fft_(double *, double *,int,int,int,int);
+void fftmx(double *,double *,int,int,int,int,int,int *,double *,double *,double *,double *,int *,int[]);
+void reals_(double *,double *,int,int);
+}
+
+
+#define DEFAULT_BUFLEN (32768)
+
+enum {ARG_NAME,ARG_INFILE,ARG_IMPFILE,ARG_OUTFILE,ARG_NARGS};
+
+void usage(const char *progname)
+{
+
+    printf("\n\n%s v1.2 RWD,CDP July 2010,2013",progname);
+	printf("\nMulti-channel FFT-based convolver");
+
+	printf("\nUsage: \n    %s [-aX][-f] infile impulsefile outfile [dry]\n"		               
+           "   -aX        : scale output amplitude by X\n"
+           "   -f         : write output as floats (no clipping)\n"
+           "  infile      : input soundfile to be processed.\n"
+           "  impulsefile : m/c soundfile or mono text file containing impulse response,\n"
+           "                  e.g. reverb or FIR filter.\n"
+           "                Text file name must have extension .txt.\n"
+		   "                  File must contain 1 column of floating point values,\n"
+		   "                  typically in the range -1.0 to 1.0.\n"
+           "              Supported channel combinations:\n"
+           "               (a) mono infile, N-channel impulsefile;\n"
+           "               (b) channels are the same;\n"
+           "               (c) mono impulsefile, N-channel infile.\n"
+           "  [dry]       :  set dry/wet mix (e.g. for reverb)\n"
+           "                 Range: 0.0 - 1.0,  default = 0.0\n"
+           "                 (uses sin/cos law for constant power mix)\n"
+           "Note: some recorded reverb impulses effectively include the direct signal.\n"
+           "In such cases  <dry>  need not be used\n"
+           "Where impulsefile is filter response (FIR), optimum length is power-of-two - 1.\n",progname
+		   );		   
+}
+
+
+double
+timer()
+{
+	struct timeb now;
+	double secs, ticks;	
+	ftime(&now);
+	ticks = (double)now.millitm/(1000.0);
+	secs = (double) now.time;
+
+	return secs + ticks;
+}
+
+
+void			
+stopwatch(int flag) 
+{
+	static double start;		    
+	int mins=0,hours=0;
+	double end, secs=0;
+
+	if(flag)
+		start = timer();	
+	else 
+	{
+		end    = timer(); 
+		secs   = end - start;
+		mins   = (int)(secs/60.0);
+		secs  -= mins*60.0; 
+		hours  = mins/60;
+		mins  -= hours*60;
+
+		fprintf(stderr,"\nElapsed time: ");
+		if(hours > 0)
+			fprintf(stderr,"%d h ", hours);
+		if(mins > 0)
+			fprintf(stderr,"%d m ", mins);
+		fprintf(stderr,"%2.3lf s\n\n",secs);
+	}
+}
+
+/* how do we want this to work?
+
+  (1)  multi-chan impulse:  EITHER: mono infile  OR infile with same chan count
+  (2)  mono impulse:  expanded to infile chancount
+     Therefore: need to open both infiles before deciding action
+*/
+#define PI2 (1.570796327)
+
+int main(int argc,char **argv)
+{
+
+	int fftlen = 1,imlen = 0,chans = 0,srate;
+	double scalefac = 1.0f;
+    double Ninv = 1.0;
+	int insndfile = -1,inimpfile = -1,outsndfile = -1; 
+    int nameoffset = 0;
+	PSF_PROPS inprops,outprops, impulseprops;
+    psf_format format;
+    PSF_CHPEAK  *fpeaks = NULL;
+	MYLONG peaktime;
+	int i,jj,minheader = 0,do_timer = 1;
+    int writefloats = 0;
+    int framesneeded;
+	double oneovrsr;
+    double *inmonobuf = 0;
+	double *insbuf=0,*outsbuf = 0;
+    double **inbuf = 0, **outbuf = 0;    
+	double **imbuf = 0;
+	double **overlapbuf = 0;
+    double rms = 0.0;
+    double dry = 0.0;
+    double wet = 1.0;
+#ifdef FFTTEST
+    double *anal = 0;
+#endif
+    	
+    if(argv[0][0]=='.' && argv[0][1]=='/'){
+        nameoffset  = 2;
+    }
+	
+    if(argc==2){
+        if(strcmp(argv[1],"--version")==0){
+            printf("1.2.0.\n");
+            return 0;
+        }
+    }
+
+	if(psf_init()){
+		puts("unable to start portsf\n");
+		return 1;
+	}
+
+
+	if(argc < ARG_NARGS){		 
+		usage(argv[0]+nameoffset);
+		return(1);
+	}
+
+    
+
+    while(argc > 1 && argv[1][0]=='-'){				
+		switch(argv[1][1]){
+        case 'a':
+            scalefac =  atof(&(argv[1][2]));
+            if(scalefac <=0.0){
+                printf("Error: scalefac must be positive!\n");
+                return 1;
+            }
+            break;        
+        case 'f':
+            writefloats = 1;
+            break;
+        default:
+            break;
+        }
+        argv++;
+        argc--;
+    }		
+	/* 2 legal possibilities: infile and outfile, or -I used with infile only */
+	if(argc< ARG_NARGS){
+		fprintf(stderr,"\nInsufficient arguments.");
+		usage(argv[0]+nameoffset);
+		return(1);
+	}
+    if(argc==ARG_NARGS+1){
+        double dryarg;
+        dryarg = atof(argv[ARG_NARGS]);
+        if(dryarg < 0.0 || dryarg > 1.0){
+            printf("dry value out of range.\n");
+            return 0;
+        }
+        if(dryarg==1.0){
+            printf("Warning: dry=1 copies input!\n");
+            wet = 0.0;
+        }
+        else{
+            dry = sin(PI2 * dryarg);
+            wet = cos(PI2 * dryarg);
+        }
+    }
+
+		
+
+/* TODO:  where -F[file] is combined with -A, it signifies create analysis file
+           compatible with impulse file (e.g 50% overlap, etc) */
+
+	
+	/* open infile, check props */	
+		
+	if((insndfile = psf_sndOpen(argv[ARG_INFILE],&inprops, 0)) < 0){
+		fprintf(stderr,"\nUnable to open input soundfile %s",argv[1]);
+        delete []imbuf;
+		return 1;
+	}
+    srate = inprops.srate;
+	if(srate <=0){
+		fprintf(stderr,"\nBad srate found: corrupted file?\n");
+        delete []imbuf;
+		return 1;
+	}
+	chans = inprops.chans;
+    framesneeded = psf_sndSize(insndfile);
+    if(framesneeded <= 0){
+        printf("Error in input file - no data!\n");
+        psf_sndClose(insndfile);
+        return 1;
+    }
+    /* open impulse file */
+    /* check for soundfile */
+      
+    format = psf_getFormatExt(argv[ARG_IMPFILE]);
+    if(format==PSF_FMT_UNKNOWN){  /* must be a text file */
+        FILE *fp = 0;
+        char tmp[80];
+        char* dot;
+        int size = 0,got = 0;
+        int read = 0;
+
+        dot = strrchr(argv[ARG_IMPFILE],'.');
+        if(dot==NULL || stricmp(dot,".txt") != 0){
+            fprintf(stderr,"Error: impulse text file must have .txt extension.\n");
+            return 1;
+        }
+        fp = fopen(argv[ARG_IMPFILE],"r");
+        if(fp==0){
+            printf("Cannot open impulse text file %s\n.",argv[ARG_IMPFILE]);
+            return 1;
+        }
+        /* count lines! */
+        while(fgets(tmp,80,fp) != NULL)
+            size++;
+        if(size==0){
+            printf("Impulse textfile %s has no data!.\n",argv[ARG_IMPFILE]);
+            return 1;
+        }
+        rewind(fp);
+            
+        inmonobuf = new double[size+1];
+        for(i=0;i < size;i++)  {
+            read = fscanf(fp,"%lf",&inmonobuf[i]);
+            if(read==0){
+                i--;
+                continue;
+            }
+            if(read==EOF)
+                break;
+            got++;
+        }
+        imlen = got;
+        impulseprops.chans = 1;
+        fclose(fp);
+    }
+    else{
+        if((inimpfile = psf_sndOpen(argv[ARG_IMPFILE],&impulseprops, 0))< 0){
+            fprintf(stderr,"\nUnable to open impulse file %s",argv[ARG_IMPFILE]);       
+            return 0;
+        }
+        //printf("impulse file sr = %d\n",impulseprops.srate);
+        if(srate != impulseprops.srate){
+            printf("Error: files must have same sample rate");
+            delete []imbuf;
+            return 1;
+        }
+        int  len = psf_sndSize(inimpfile);
+        if(len <= 0){
+            printf("Error in impulse soundfile - no data!\n");
+            return 1;
+        }
+        if(impulseprops.chans > 1){
+            if(! (chans == 1 || chans == impulseprops.chans)){
+                fprintf(stderr,"\nChannel mismatch.\n    Infile must be mono or have same channels as irfile.\n");
+                return(1);
+            }
+            chans = impulseprops.chans;
+        }    
+    }
+
+    imbuf = new double*[chans];
+
+    /* if impulse file is mono, we will need to copy data to the other buffers */
+    // allocate  impulseprops.chans buffers ; may need more
+
+    if(inimpfile >=0){
+	    if((imlen = genimpframe2(inimpfile,&imbuf,&rms,impulseprops.chans, scalefac)) == 0){
+	    	fprintf(stderr,"Error reading impulse file\n");
+	    	return(1);
+	    }
+    }
+    
+    else if(imlen){        
+        genimpframe1(inmonobuf,&imbuf,imlen,scalefac);       
+    }
+    printf("mean rms level = %.4lf (%.2lf dB)\n",rms, 20.0 * log10(rms));
+    
+
+    framesneeded +=  imlen;        /* can  close outfile, when this length reached */
+
+    while(fftlen < imlen*2)			/* convolution style - double-length */
+		fftlen <<= 1;
+    double norm = sqrt(2.0);       /* relative: rms of 0dBFS sine is 0.707 = 1/root2 */
+    // scale factor: most of this sheer guesswork!
+    Ninv = fftlen;
+    Ninv /= sqrt(imlen*2);
+    Ninv *= norm;
+    Ninv /= imlen;
+    
+    // take simple avg of  rms for adjustment factor.
+    // may not adequately represent some responses, e.g. with strong attack/earlies, soft reverb tail */
+    double rmsdif2 =  (-20.0 * log10(rms)) * 0.5;
+    double rmsadjust = pow(10, rmsdif2/20.0);
+#ifdef _DEBUG
+    printf("rescaling factor = %.6e\n",Ninv);   
+    printf("rmsadjust = %.4lf\n",rmsadjust);
+#endif
+   Ninv /= rmsadjust; 
+
+
+    /* copy buffers if required */
+    for(i = impulseprops.chans;i < chans;i++){
+        imbuf[i] = new double[fftlen+2];
+        memcpy(imbuf[i],imbuf[0],(fftlen+2)* sizeof(double));
+    }         
+    oneovrsr = 1.0 / (double) srate;
+    
+   
+	/* main i/o buffers */
+	
+	insbuf     = new double [fftlen  * chans];
+    outsbuf    = new double [fftlen  * chans];
+	inbuf      = new double*[chans];
+    outbuf     = new double*[chans];  /* overlap-add buffers */
+    overlapbuf = new double*[chans];
+    
+    for(i=0;i < chans;i++){
+        inbuf[i]  = new double[fftlen+2];   // main in-place buf 
+        outbuf[i] = new double[fftlen+2];        
+	    /* channel buffers */	
+        overlapbuf[i] = new double[fftlen/2];
+        memset(overlapbuf[i],0,sizeof(double)*(fftlen/2));        
+    }
+#ifdef FFTTEST
+    anal = new double[fftlen+2];
+#endif
+    
+    for(i=0;i < chans;i++){		    
+	    memset(inbuf[i], 0, sizeof(double) * (fftlen+2));
+        memset(outbuf[i],0,sizeof(double) * (fftlen+2));
+    }
+
+    
+	/* use generic init function */
+/*bool phasevocoder::init(int outsrate,int fftlen,int winlen,int decfac,float scalefac,
+						pvoc_scaletype stype,pvoc_windowtype wtype,pvocmode mode)
+*/
+	
+	/*create output sfile */
+	psf_stype smptype;
+	psf_format outformat;
+	/* will it be aiff, wav, etc? */
+	outformat = psf_getFormatExt(argv[ARG_OUTFILE]);
+	if(outformat==PSF_FMT_UNKNOWN){
+		fprintf(stderr,"\nOutfile name has unrecognized extension.");
+        delete []imbuf;
+		return(1);
+	}
+		
+		
+	smptype = inprops.samptype;
+    if(writefloats)
+        smptype= PSF_SAMP_IEEE_FLOAT;
+
+	/*the one gotcha: if output is floats and format is aiff, change to aifc */
+	if(smptype==PSF_SAMP_IEEE_FLOAT && outformat==PSF_AIFF){
+		fprintf(stderr,"Warning: AIFF output written as AIFC for float samples\n");
+		outformat = PSF_AIFC;
+	}
+
+	outprops          = inprops;
+	outprops.chans    = chans;
+	outprops.srate    = srate;
+	outprops.format   = outformat;
+	outprops.samptype = smptype;
+	outprops.chformat = STDWAVE;
+    /* if irfile is BFormat, need to set outfile fmt likewise */
+    if(impulseprops.chformat==MC_BFMT)
+        outprops.chformat = MC_BFMT;
+    /* I suppose people will want automatic decoding too...grrr */
+
+    fpeaks = new PSF_CHPEAK[chans];
+	if(fpeaks==NULL){
+		puts("no memory for fpeak data buffer\n");
+		return 1;
+	}
+
+	/*last arg: no clipping of f/p samples: we use PEAK chunk */
+	if((outsndfile = psf_sndCreate(argv[ARG_OUTFILE],&outprops,0,minheader,PSF_CREATE_RDWR)) <0 ){
+		fprintf(stderr,"\nUnable to open outfile %s",argv[ARG_OUTFILE]);
+        delete []imbuf;
+		return(1);
+	}
+				
+
+	
+	printf("\n");
+	
+	
+	stopwatch(1);
+
+	int written,thisblock,framesread;
+    int frameswritten = 0;
+	double intime= 0.0;
+				
+	while((framesread = psf_sndReadDoubleFrames(insndfile,insbuf,imlen)) > 0){        
+		written = thisblock =  0;
+
+		for(i = framesread * inprops.chans; i< fftlen * inprops.chans; i++)
+			insbuf[i] = 0.0f;
+        framesread = imlen;
+		memset(inbuf[0],0,(fftlen+2) * sizeof(double));
+        if(chans == inprops.chans)  {
+            /* must clean buffers! */
+		    for(i=0;i < chans;i++)
+                memset(inbuf[i],0,(fftlen+2) * sizeof(double));
+            mc_split(insbuf,inbuf,imlen * chans, chans);
+
+        }
+        else{
+            for(i=0;i < chans;i++) {
+                memset(inbuf[i],0,(fftlen+2) * sizeof(double));
+                memcpy(inbuf[i],insbuf,imlen * sizeof(double));
+                memset(outbuf[i],0,sizeof(double) * (fftlen+2));
+
+            }
+        }
+        if(impulseprops.chans==1){
+			
+            for(jj = 0; jj < chans;jj++){
+#ifdef FFTTEST
+                int zz;
+                memcpy(anal,inbuf[jj],(fftlen+2) * sizeof(double));
+                fft_(anal,anal+1,1,fftlen/2,1,-2);
+	            reals_(anal,anal+1,fftlen/2,-2);
+                for(zz=0;zz < fftlen+2;zz++)
+                    anal[zz] *= 0.001;
+#endif
+                //rfftwnd_one_real_to_complex(forward_plan[jj],inbuf[jj],NULL);	
+				double *danal = inbuf[jj];
+				fft_(danal,danal+1,1,fftlen/2,1,-2);
+				reals_(danal,danal+1,fftlen/2,-2);
+                complexmult(inbuf[jj],imbuf[0],fftlen+2);                
+			    //rfftwnd_one_complex_to_real(inverse_plan[jj],(fftw_complex * )inbuf[jj],NULL);
+				reals_(danal,danal+1,fftlen/2,2);
+				fft_(danal,danal+1,1,fftlen/2,1,2);
+            }           
+        }
+        else{
+            for(jj = 0; jj < chans;jj++){
+				//rfftwnd_one_real_to_complex(forward_plan[jj],inbuf[jj],NULL);                   				
+                double *danal = inbuf[jj];
+				fft_(danal,danal+1,1,fftlen/2,1,-2);
+				reals_(danal,danal+1,fftlen/2,-2);
+				complexmult(inbuf[jj],imbuf[jj],fftlen+2);                
+				//rfftwnd_one_complex_to_real(inverse_plan[jj],(fftw_complex * )inbuf[jj],NULL); 
+				reals_(danal,danal+1,fftlen/2,2);
+				fft_(danal,danal+1,1,fftlen/2,1,2);
+            }            
+        }
+        
+        /* overlap-add  for each channel */
+        /* TODO: verify  use of imlen here - should it be fftlen/2 -1 ? */
+        for(jj=0;jj < chans;jj++){
+            for(i=0;i < imlen;i++) {
+                outbuf[jj][i] = inbuf[jj][i] + overlapbuf[jj][i];
+                overlapbuf[jj][i] = inbuf[jj][i+imlen];
+            }
+        }
+		mc_interl(outbuf,outsbuf,imlen,chans);
+
+        if(inprops.chans == chans){
+            for(i=0;i < framesread; i++) {
+                for(jj=0;jj < chans;jj++){
+                    int  outindex = i*chans + jj;            
+                    outsbuf[outindex] *= Ninv;
+                    outsbuf[outindex] *= wet;
+                    outsbuf[outindex] += dry * insbuf[outindex];
+                }
+            }
+        }
+        /* elso mono input */
+        else {
+            for(i=0;i < framesread; i++) {
+                for(jj=0;jj < chans;jj++){
+                    int  outindex = i*chans + jj;
+                    double inval = dry *  insbuf[i];
+                    outsbuf[outindex] *= Ninv;
+                    outsbuf[outindex] *= wet;
+                    outsbuf[outindex] += inval;
+                }
+            }
+        }
+
+        if((written = psf_sndWriteDoubleFrames(outsndfile,outsbuf,framesread)) != framesread){
+		    fprintf(stderr,"\nerror writing to outfile");
+		    return(1);		               
+        }
+        frameswritten += written;
+
+		if(do_timer){
+			intime += (double)framesread * oneovrsr;
+			printf("Input time: %.2lf\r",intime);
+		}    
+    }
+    /* complete tail */
+
+    if(frameswritten < framesneeded){
+        // TODO: imlen? see above
+        mc_interl(overlapbuf,outsbuf,imlen,chans);
+        int towrite = framesneeded - frameswritten; 
+        for(i=0;i < towrite * chans; i++) {
+            outsbuf[i] *= Ninv;
+            outsbuf[i] *= wet;
+        }
+        if((written = psf_sndWriteDoubleFrames(outsndfile,outsbuf,towrite)) != towrite){
+	        fprintf(stderr,"\nerror writing to outfile");
+	        return(1);
+	    }
+    }
+
+	stopwatch(0);
+    if(psf_sndReadPeaks( outsndfile,fpeaks,&peaktime)){
+        double peakmax = 0.0;
+		printf("PEAK values:\n");
+		for(i=0; i < chans; i++)	{
+            peakmax = fpeaks[i].val > peakmax ? fpeaks[i].val : peakmax;
+			if(fpeaks[i].val < 0.0001f)
+				printf("CH %d:\t%e (%.2lf dB)\tat frame %10u :\t%.4lf secs\n",i+1,
+				fpeaks[i].val,20.0*log10(fpeaks[i].val),fpeaks[i].pos,(double) (fpeaks[i].pos) / (double)srate);
+			else
+				printf("CH %d:\t%.4lf (%.2lf dB)\tat frame %10u :\t%.4lf secs\n",i+1,
+				fpeaks[i].val,20.0 * log10(fpeaks[i].val),fpeaks[i].pos,(double) (fpeaks[i].pos) / (double)srate);		
+		}
+        if(peakmax > 1.0)
+            printf("Overflows!  Rescale by %.10lf\n",1.0 / peakmax);
+	}
+
+	if(insndfile >=0)
+		psf_sndClose(insndfile);
+	if(outsndfile >=0)
+		psf_sndClose(outsndfile);
+    if(inimpfile >=0)
+	    psf_sndClose(inimpfile);
+	psf_finish();
+    delete [] fpeaks;
+	if(insbuf)
+		delete [] insbuf;
+	if(outsbuf)
+		delete [] outsbuf;
+    
+
+    for(i=0;i < chans;i++){
+        if(inbuf){       
+		    delete [] inbuf[i];
+        }           
+        if(outbuf){
+		    delete [] outbuf[i];
+        }       	
+        if(imbuf){        
+	        delete [] imbuf[i];
+        }
+        if(overlapbuf)
+            delete [] overlapbuf[i];
+    }
+    delete [] outbuf;
+    delete [] inbuf;
+    delete [] imbuf;
+    delete [] overlapbuf;
+	return 0;
+}
+
+
+// insize is raw samplecount,buflen is insize/chans
+void mc_split(double* inbuf,double** out,int insize,int chans)
+{
+    int i,j,buflen = insize/chans;
+    double* pinbuf;
+
+    
+    for(j=0;j < chans;j++){
+        pinbuf = inbuf+j;
+        for(i=0;i < buflen;i++){
+            out[j][i] = *pinbuf;
+            pinbuf += chans;
+        }
+    }
+}
+
+
+/* insize is m/c frame count */
+void mc_interl(double** in,double* out,int insize,int chans)
+{
+    int i,j;
+    double* poutbuf;
+
+    for(j = 0;j < chans;j++){
+        poutbuf = out+j;
+        for(i=0;i < insize;i++){
+            *poutbuf = in[j][i];
+            poutbuf += chans;
+        }
+    }
+}
+
+/* OR:  apply scalefac to impulse responses */
+void complexmult(double *frame,const double *impulse,int length)
+{
+	double re,im;
+	
+	int i,j;
+	
+
+	for(i=0,j = 1;i < length;i+=2,j+=2){
+		re = frame[i] * impulse[i] - frame[j] * impulse[j];
+		im = frame[i] * impulse[j] + frame[j]* impulse[i];
+		frame[i] = re;
+		frame[j] = im;
+	}
+}
+
+#ifdef unix
+int stricmp(const char *a, const char *b)
+{
+	while(*a != '\0' && *b != '\0') {
+		int ca = islower(*a) ? toupper(*a) : *a;
+		int cb = islower(*b) ? toupper(*b) : *b;
+        
+		if(ca < cb)
+			return -1;
+		if(ca > cb)
+			return 1;
+        
+		a++;
+		b++;
+	}
+	if(*a == '\0' && *b == '\0')
+		return 0;
+	if(*a != '\0')
+		return 1;
+	return -1;
+}
+
+int
+strnicmp(const char *a, const char *b, const int length)
+{
+	int len = length;
+    
+	while(*a != '\0' && *b != '\0') {
+		int ca = islower(*a) ? toupper(*a) : *a;
+		int cb = islower(*b) ? toupper(*b) : *b;
+        
+		if(len-- < 1)
+			return 0;
+        
+		if(ca < cb)
+			return -1;
+		if(ca > cb)
+			return 1;
+        
+		a++;
+		b++;
+	}
+	if(*a == '\0' && *b == '\0')
+		return 0;
+	if(*a != '\0')
+		return 1;
+	return -1;
+}
+#endif
+
+
+
+
+
+

dev/externals/fastconv/genrespframe2.cpp

@@ -0,0 +1,235 @@
+/*
+ * 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
+ *
+ */
+ 
+/* genrespframe2.cpp */
+/* generate m/c pvoc frames containing impulse response */
+extern "C"
+{
+#include <portsf.h>
+}
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+
+
+#ifdef _DEBUG
+#include <assert.h>
+#endif
+
+/* deinterleave input buffer of insize samps to array of chans buffers */
+void mc_split(double* inbuf,double** out,int insize,int chans);
+extern "C"{
+	void fft_(double *, double *,int,int,int,int);
+	void fftmx(double *,double *,int,int,int,int,int,int *,double *,double *,double *,double *,int *,int[]);
+	void reals_(double *,double *,int,int);
+}
+
+/*  read impulse data from soundfile */
+/* return 0 for error, or size of impulse */
+/*we assume pvsys initialized by caller */
+int genimpframe2(int ifd, double*** outbufs, double* rms,int imchans, double scalefac)
+{
+	int i,j;	
+	double *insbuf = 0;
+    double **inbuf = *outbufs;
+	int fftsize = 1;
+    int inframes;
+    double*  rmsfac = 0;
+    double* ampsum = 0;
+    double  maxrmsfac = 0.0;
+#ifdef _DEBUG
+    double *re,*im,*mag;
+#endif
+	inframes = psf_sndSize(ifd);  // m/c frames
+	if(inframes <= 0)
+        return 0;
+	while(fftsize < inframes*2)			/* convolution style - double-length */
+		fftsize <<= 1;
+	printf("impulse length = %d frames, fftsize = %d\n",inframes,fftsize);
+    insbuf = new double[inframes * imchans];
+#ifdef _DEBUG
+    re = new double[fftsize/2+1];
+    im = new double[fftsize/2+1];
+    mag = new double[fftsize/2+1];
+#endif
+    for(i=0;i < imchans;i++){	
+	    inbuf[i] = new double[fftsize+2];
+	    if(inbuf[i]==NULL){
+		    puts("no memory for file buffer!\n");                    
+		    return 0;
+	    }
+	    memset(inbuf[i],0,(fftsize+2)* sizeof(double));
+    }
+
+
+	if(inframes != psf_sndReadDoubleFrames(ifd,insbuf,inframes)){
+		fprintf(stderr,"Error reading impulse data\n");
+		
+		delete [] insbuf;
+        // and do other cleanup!
+        
+		return 0;
+	}
+    rmsfac = new double[imchans];
+    ampsum = new double[imchans];
+    for(i=0;i< imchans;i++) {
+        rmsfac[i] = 0.0;
+        ampsum[i] = 0.0;
+    }
+    /* do user scaling first */
+    for(i = 0;i < inframes * imchans;i++)
+        insbuf[i] *= scalefac;
+    /* now try to adjust rms  */
+
+    for(i = 0;i < inframes;i++) {
+        for(j=0;j < imchans;j++){
+            double val = insbuf[i*imchans + j];
+            ampsum[j] += fabs(val);
+            rmsfac[j] += val*val; 
+        }
+    }
+    for(j=0;j < imchans;j++){
+        //    rmsfac = sqrt(rmsfac);
+        rmsfac[j] /= inframes;
+        rmsfac[j] = sqrt(rmsfac[j]);
+        ampsum[j] /= inframes;
+        if(rmsfac[j] > maxrmsfac)
+            maxrmsfac = rmsfac[j];
+#ifdef _DEBUG
+        if(ampsum[j] > maxampsum)
+            maxampsum = ampsum[j];
+#endif
+    }
+    /* do the rescaling! */
+          
+#ifdef _DEBUG        
+        printf("ampsum = %.4f\n",maxampsum);
+#endif
+    
+    // now deinterleave to each inbuf
+    mc_split(insbuf,inbuf,inframes * imchans,imchans);
+    for(i=0;i < imchans;i++){
+		double *anal = inbuf[i];
+	    //rfftwnd_one_real_to_complex(forward_plan[i],inbuf[i],NULL);
+		fft_(anal,anal+1,1,fftsize/2,1,-2);
+		reals_(anal,anal+1,fftsize/2,-2);
+    }
+#ifdef _DEBUG
+    /* in order to look at it all */
+    double  magmax = 0.0;
+    double magsum = 0.0;
+    for(i=0;i < fftsize/2+1;i++){
+        double thisre, thisim;
+        thisre =inbuf[0][i*2];        
+        thisim = inbuf[0][i*2+1]; 
+        re[i] = thisre;
+        im[i] = thisim;
+        mag[i] = sqrt(thisre*thisre + thisim*thisim); 
+        magsum += (mag[i] * mag[i]);
+        if(mag[i] > magmax)
+            magmax = mag[i];
+    }
+    magsum = sqrt(magsum);
+    printf("maxamp of FFT = %.4f\n",magmax);
+    printf("mean level of FFT = %.4f\n",magsum / (fftsize/2+1));
+#endif	
+    
+    delete [] rmsfac;
+    delete [] ampsum;
+
+#ifdef _DEBUG
+    delete [] re;
+    delete [] im;
+    delete [] mag;
+#endif
+    *rms = maxrmsfac;
+	return inframes;
+}
+/* convert from input double buffer (read from mono text impulse file) */
+int genimpframe1(double *insbuf, double*** outbuf, int npoints, double scalefac)
+{
+	int i;	
+    double **inbuf = *outbuf;
+	int fftsize = 1;
+    int insamps;
+#ifdef _DEBUG
+    double *re,*im,*mag;
+#endif
+	insamps = npoints;  // m/c frames
+	if(insamps <= 0)
+        return 0;
+	while(fftsize < insamps*2)			/* convolution style - double-length */
+		fftsize <<= 1;
+	printf("infile size = %d,impulse framesize = %d\n",insamps,fftsize);
+    
+#ifdef _DEBUG
+    re = new double[fftsize/2+1];
+    im = new double[fftsize/2+1];
+    mag = new double[fftsize/2+1];
+#endif
+    	
+	inbuf[0] = new double[fftsize+2];
+	if(inbuf[0]==NULL){
+	    puts("no memory for file buffer!\n");                    
+	    return 0;
+	}
+	memset(inbuf[0],0,(fftsize+2)* sizeof(double)); 
+#ifdef _DEBUG
+    double ampsum = 0.0;
+#endif
+    for(i = 0;i < insamps;i++) {
+#ifdef _DEBUG
+        ampsum += fabs(insbuf[i]);
+#endif        
+        insbuf[i] *= scalefac;
+    }
+#ifdef _DEBUG
+    printf("amplitude sum of impulse file = %f\n",ampsum);
+#endif
+    memcpy(inbuf[0],insbuf,npoints* sizeof(double));
+	//rfftwnd_one_real_to_complex(forward_plan,inbuf[0],NULL);
+	double *anal = inbuf[0];
+	fft_(anal,anal+1,1,fftsize/2,1,-2);
+	reals_(anal,anal+1,fftsize/2,-2);
+    
+#ifdef _DEBUG
+    /* in order to look at it all */
+    for(i=0;i < fftsize/2+1;i++){
+        double thisre, thisim;
+        thisre =inbuf[0][i*2];        
+        thisim = inbuf[0][i*2+1]; 
+        re[i] = thisre;
+        im[i] = thisim;
+        mag[i] = sqrt(thisre*thisre + thisim*thisim);
+    }
+#endif	
+    
+    
+#ifdef _DEBUG
+    delete [] re;
+    delete [] im;
+    delete [] mag;
+#endif
+   
+
+	return insamps;
+}

dev/externals/fastconv/mxfftd.c

@@ -0,0 +1,883 @@
+/*  
+    mxfftd.c: double precision version of mxfft.c
+	 
+	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
+
+*/
+
+/* This program converted from the FORTRAN routines by Singleton in
+ * Section 1.4 of  "Programs for Digital Signal Processing", IEEE Press, 1979.
+ *  Conversion by Trevor Wishart and Keith Henderson, York Univ.
+ */
+
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+
+void fft_(double *, double *, int, int, int, int);
+void fftmx(double *, double *, int, int, int, int, int,
+      int*, double *, double *, double *, double *, int *, int[]);
+void reals_(double *, double *, int, int);
+
+/*
+ *-----------------------------------------------------------------------
+ * subroutine:  fft
+ * multivariate complex fourier transform, computed in place
+ * using mixed-radix fast fourier transform algorithm.
+ *-----------------------------------------------------------------------
+ *
+ *      this is the call from C:
+ *              fft_(anal,banal,&one,&N2,&one,&mtwo);
+ *      CHANGED TO:-
+ *              fft_(anal,banal,one,N2,one,mtwo);
+ */
+
+void
+fft_(double *a, double *b, int nseg, int n, int nspn, int isn)
+  /*    *a,       pointer to array 'anal'  */
+  /*    *b;       pointer to array 'banal' */
+{
+    int nfac[16];               /*  These are one bigger than needed   */
+                                /*  because wish to use Fortran array  */
+                                /* index which runs 1 to n, not 0 to n */
+
+    int         m = 0,
+                nf,
+                k,
+                kt,
+                ntot,
+                j,
+                jj,
+                maxf, maxp=-1;
+
+/* work space pointers */
+    double       *at, *ck, *bt, *sk;
+    int *np;
+
+
+/* reduce the pointers to input arrays - by doing this, FFT uses FORTRAN
+   indexing but retains compatibility with C arrays */
+    a--;        b--;
+
+/*
+ * determine the factors of n
+ */
+    k=nf=abs(n);
+    if (nf==1)
+      return;
+
+    nspn=abs(nf*nspn);
+    ntot=abs(nspn*nseg);
+
+    if (isn*ntot == 0) {
+      printf("\nerror - zero in fft parameters %d %d %d %d",
+              nseg, n, nspn, isn);
+      return;
+    }
+    for (m=0; !(k%16); nfac[++m]=4,k/=16);
+    for (j=3,jj=9; jj<=k; j+=2,jj=j*j)
+      for (; !(k%jj); nfac[++m]=j,k/=jj);
+
+    if (k<=4) {
+      kt = m;
+      nfac[m+1] = k;
+      if (k != 1)
+        m++;
+    }
+    else {
+      if (k%4==0) {
+        nfac[++m]=2;
+        k/=4;
+      }
+
+      kt = m;
+      maxp = (kt+kt+2 > k-1 ? kt+kt+2 : k-1);
+      for (j=2; j<=k; j=1+((j+1)/2)*2)
+        if (k%j==0) {
+          nfac[++m]=j;
+          k/=j;
+        }
+    }
+    if (m <= kt+1)
+      maxp = m + kt + 1;
+    if (m+kt > 15) {
+      printf("\nerror - fft parameter n has more than 15 factors : %d", n);
+      return;
+    }
+    if (kt!=0) {
+      j = kt;
+      while (j)
+        nfac[++m]=nfac[j--];
+    }
+    maxf = nfac[m-kt];
+    if (kt > 0 && maxf <nfac[kt])
+      maxf = nfac[kt];
+
+/*  allocate workspace - assume no errors! */
+    at = (double *) calloc(maxf,sizeof(double));
+    ck = (double *) calloc(maxf,sizeof(double));
+    bt = (double *) calloc(maxf,sizeof(double));
+    sk = (double *) calloc(maxf,sizeof(double));
+    np = (int *) calloc(maxp,sizeof(int));
+
+/* decrement pointers to allow FORTRAN type usage in fftmx */
+    at--;       bt--;   ck--;   sk--;   np--;
+
+/* call fft driver */
+
+    fftmx(a,b,ntot,nf,nspn,isn,m,&kt,at,ck,bt,sk,np,nfac);
+
+/* restore pointers before releasing */
+    at++;       bt++;   ck++;   sk++;   np++;
+
+/* release working storage before returning - assume no problems */
+    free(at);
+    free(sk);
+    free(bt);
+    free(ck);
+    free(np);
+    return;
+}
+
+/*
+ *-----------------------------------------------------------------------
+ * subroutine:  fftmx
+ * called by subroutine 'fft' to compute mixed-radix fourier transform
+ *-----------------------------------------------------------------------
+ */
+void
+fftmx(double *a, double *b, int ntot, int n, int nspan, int isn, int m,
+      int *kt, double *at, double *ck, double *bt, double *sk, int *np, int nfac[])
+{
+    int i,inc,
+      j,jc,jf, jj,
+      k, k1, k2, k3, k4,
+      kk,klim,ks,kspan, kspnn,
+      lim,
+      maxf,mm,
+      nn,nt;
+    double  aa, aj, ajm, ajp, ak, akm, akp,
+      bb, bj, bjm, bjp, bk, bkm, bkp,
+      c1, c2, c3, c72, cd,
+      dr,
+      rad,
+      sd, s1, s2, s3, s72, s120;
+
+    double      xx;     /****** ADDED APRIL 1991 *********/
+    inc=abs(isn);
+    nt = inc*ntot;
+    ks = inc*nspan;
+/******************* REPLACED MARCH 29: ***********************
+                                        rad = atan((double)1.0);
+**************************************************************/
+    rad = 0.785398163397448278900;
+/******************* REPLACED MARCH 29: ***********************
+                                        s72 = rad/0.625;
+                                        c72 = cos(s72);
+                                        s72 = sin(s72);
+**************************************************************/
+    c72 = 0.309016994374947451270;
+    s72 = 0.951056516295153531190;
+/******************* REPLACED MARCH 29: ***********************
+                                        s120 = sqrt((double)0.75);
+**************************************************************/
+    s120 = 0.866025403784438707600;
+
+/* scale by 1/n for isn > 0 ( reverse transform ) */
+
+    if (isn < 0) {
+      s72 = -s72;
+      s120 = -s120;
+      rad = -rad;}
+    else {
+      ak = 1.0/(double)n;
+      for (j=1; j<=nt;j += inc) {
+        a[j] *= (double)ak;
+        b[j] *= (double)ak;
+      }
+    }
+    kspan = ks;
+    nn = nt - inc;
+    jc = ks/n;
+
+/* sin, cos values are re-initialised each lim steps  */
+
+    lim = 32;
+    klim = lim * jc;
+    i = 0;
+    jf = 0;
+    maxf = m - (*kt);
+    maxf = nfac[maxf];
+    if ((*kt) > 0 && maxf < nfac[*kt])
+      maxf = nfac[*kt];
+
+/*
+ * compute fourier transform
+ */
+
+ lbl40:
+    dr = (8.0 * (double)jc)/((double)kspan);
+/*************************** APRIL 1991 POW & POW2 not WORKING.. REPLACE *******
+                    cd = 2.0 * (pow2 ( sin((double)0.5 * dr * rad)) );
+*******************************************************************************/
+    xx =  sin((double)0.5 * dr * rad);
+    cd = 2.0 * xx * xx;
+    sd = sin(dr * rad);
+    kk = 1;
+    if (nfac[++i]!=2) goto lbl110;
+/*
+ * transform for factor of 2 (including rotation factor)
+ */
+    kspan /= 2;
+    k1 = kspan + 2;
+    do {
+      do {
+        k2 = kk + kspan;
+        ak = a[k2];
+        bk = b[k2];
+        a[k2] = (a[kk]) - (double)ak;
+        b[k2] = (b[kk]) - (double)bk;
+        a[kk] = (a[kk]) + (double)ak;
+        b[kk] = (b[kk]) + (double)bk;
+        kk = k2 + kspan;
+      } while (kk <= nn);
+      kk -= nn;
+    } while (kk <= jc);
+    if (kk > kspan) goto lbl350;
+ lbl60:
+    c1 = 1.0 - cd;
+    s1 = sd;
+    mm = (k1/2 < klim ? k1/2 :klim);
+    goto lbl80;
+ lbl70:
+    ak = c1 - ((cd*c1)+(sd*s1));
+    s1 = ((sd*c1)-(cd*s1)) + s1;
+    c1 = ak;
+lbl80:
+    do {
+      do {
+        k2 = kk + kspan;
+        ak = a[kk] - a[k2];
+        bk = b[kk] - b[k2];
+        a[kk] = a[kk] + a[k2];
+        b[kk] = b[kk] + b[k2];
+        a[k2] = (double)((c1 * ak) - (s1 * bk));
+        b[k2] = (double)((s1 * ak) + (c1 * bk));
+        kk = k2 + kspan;
+      } while (kk < nt);
+      k2 = kk - nt;
+      c1 = -c1;
+      kk = k1 - k2;
+    } while (kk > k2);
+    kk += jc;
+    if (kk <= mm) goto lbl70;
+    if (kk < k2)  goto lbl90;
+    k1 += (inc + inc);
+    kk = ((k1-kspan)/2) + jc;
+    if (kk <= (jc+jc)) goto lbl60;
+    goto lbl40;
+ lbl90:
+    s1 = ((double)((kk-1)/jc)) * dr * rad;
+    c1 = cos(s1);
+    s1 = sin(s1);
+    mm = (k1/2 < mm+klim ? k1/2 : mm+klim);
+    goto lbl80;
+/*
+ * transform for factor of 3 (optional code)
+ */
+
+
+ lbl100:
+    k1 = kk + kspan;
+    k2 = k1 + kspan;
+    ak = a[kk];
+    bk = b[kk];
+    aj = a[k1] + a[k2];
+    bj = b[k1] + b[k2];
+    a[kk] = (double)(ak + aj);
+    b[kk] = (double)(bk + bj);
+    ak += (-0.5 * aj);
+    bk += (-0.5 * bj);
+    aj = (a[k1] - a[k2]) * s120;
+    bj = (b[k1] - b[k2]) * s120;
+    a[k1] = (double)(ak - bj);
+    b[k1] = (double)(bk + aj);
+    a[k2] = (double)(ak + bj);
+    b[k2] = (double)(bk - aj);
+    kk = k2 + kspan;
+    if (kk < nn)     goto lbl100;
+    kk -= nn;
+    if (kk <= kspan) goto lbl100;
+    goto lbl290;
+
+/*
+ * transform for factor of 4
+ */
+
+ lbl110:
+    if (nfac[i] != 4) goto lbl230;
+    kspnn = kspan;
+    kspan = kspan/4;
+ lbl120:
+    c1 = 1.0;
+    s1 = 0;
+    mm = (kspan < klim ? kspan : klim);
+    goto lbl150;
+ lbl130:
+    c2 = c1 - ((cd*c1)+(sd*s1));
+    s1 = ((sd*c1)-(cd*s1)) + s1;
+/*
+ * the following three statements compensate for truncation
+ * error.  if rounded arithmetic is used, substitute
+ * c1=c2
+ *
+ * c1 = (0.5/(pow2(c2)+pow2(s1))) + 0.5;
+ * s1 = c1*s1;
+ * c1 = c1*c2;
+ */
+    c1 = c2;
+ lbl140:
+    c2 = (c1 * c1) - (s1 * s1);
+    s2 = c1 * s1 * 2.0;
+    c3 = (c2 * c1) - (s2 * s1);
+    s3 = (c2 * s1) + (s2 * c1);
+ lbl150:
+    k1 = kk + kspan;
+    k2 = k1 + kspan;
+    k3 = k2 + kspan;
+    akp = a[kk] + a[k2];
+    akm = a[kk] - a[k2];
+    ajp = a[k1] + a[k3];
+    ajm = a[k1] - a[k3];
+    a[kk] = (double)(akp + ajp);
+    ajp = akp - ajp;
+    bkp = b[kk] + b[k2];
+    bkm = b[kk] - b[k2];
+    bjp = b[k1] + b[k3];
+    bjm = b[k1] - b[k3];
+    b[kk] = (double)(bkp + bjp);
+    bjp = bkp - bjp;
+    if (isn < 0) goto lbl180;
+    akp = akm - bjm;
+    akm = akm + bjm;
+    bkp = bkm + ajm;
+    bkm = bkm - ajm;
+    if (s1 == 0.0) goto lbl190;
+ lbl160:
+    a[k1] = (double)((akp*c1) - (bkp*s1));
+    b[k1] = (double)((akp*s1) + (bkp*c1));
+    a[k2] = (double)((ajp*c2) - (bjp*s2));
+    b[k2] = (double)((ajp*s2) + (bjp*c2));
+    a[k3] = (double)((akm*c3) - (bkm*s3));
+    b[k3] = (double)((akm*s3) + (bkm*c3));
+    kk = k3 + kspan;
+    if (kk <= nt)   goto lbl150;
+ lbl170:
+    kk -= (nt - jc);
+    if (kk <= mm)   goto lbl130;
+    if (kk < kspan) goto lbl200;
+    kk -= (kspan - inc);
+    if (kk <= jc)   goto lbl120;
+    if (kspan==jc)  goto lbl350;
+    goto lbl40;
+lbl180:
+    akp = akm + bjm;
+    akm = akm - bjm;
+    bkp = bkm - ajm;
+    bkm = bkm + ajm;
+    if (s1 != 0.0)  goto lbl160;
+ lbl190:
+    a[k1] = (double)akp;
+    b[k1] = (double)bkp;
+    a[k2] = (double)ajp;
+    b[k2] = (double)bjp;
+    a[k3] = (double)akm;
+    b[k3] = (double)bkm;
+    kk = k3 + kspan;
+    if (kk <= nt) goto lbl150;
+    goto lbl170;
+ lbl200:
+    s1 = ((double)((kk-1)/jc)) * dr * rad;
+    c1 = cos(s1);
+    s1 = sin(s1);
+    mm = (kspan < mm+klim ? kspan : mm+klim);
+    goto lbl140;
+
+/*
+ * transform for factor of 5 (optional code)
+ */
+
+ lbl210:
+    c2 = (c72*c72) - (s72*s72);
+    s2 = 2.0 * c72 * s72;
+ lbl220:
+    k1 = kk + kspan;
+    k2 = k1 + kspan;
+    k3 = k2 + kspan;
+    k4 = k3 + kspan;
+    akp = a[k1] + a[k4];
+    akm = a[k1] - a[k4];
+    bkp = b[k1] + b[k4];
+    bkm = b[k1] - b[k4];
+    ajp = a[k2] + a[k3];
+    ajm = a[k2] - a[k3];
+    bjp = b[k2] + b[k3];
+    bjm = b[k2] - b[k3];
+    aa = a[kk];
+    bb = b[kk];
+    a[kk] = (double)(aa + akp + ajp);
+    b[kk] = (double)(bb + bkp + bjp);
+    ak = (akp*c72) + (ajp*c2) + aa;
+    bk = (bkp*c72) + (bjp*c2) + bb;
+    aj = (akm*s72) + (ajm*s2);
+    bj = (bkm*s72) + (bjm*s2);
+    a[k1] = (double)(ak - bj);
+    a[k4] = (double)(ak + bj);
+    b[k1] = (double)(bk + aj);
+    b[k4] = (double)(bk - aj);
+    ak = (akp*c2) + (ajp*c72) + aa;
+    bk = (bkp*c2) + (bjp*c72) + bb;
+    aj = (akm*s2) - (ajm*s72);
+    bj = (bkm*s2) - (bjm*s72);
+    a[k2] = (double)(ak - bj);
+    a[k3] = (double)(ak + bj);
+    b[k2] = (double)(bk + aj);
+    b[k3] = (double)(bk - aj);
+    kk = k4 + kspan;
+    if (kk < nn)     goto lbl220;
+    kk -= nn;
+    if (kk <= kspan) goto lbl220;
+    goto lbl290;
+
+/*
+ * transform for odd factors
+ */
+
+ lbl230:
+    k = nfac[i];
+    kspnn = kspan;
+    kspan /= k;
+    if (k==3)   goto lbl100;
+    if (k==5)   goto lbl210;
+    if (k==jf)  goto lbl250;
+    jf = k;
+    s1 = rad/(((double)(k))/8.0);
+    c1 = cos(s1);
+    s1 = sin(s1);
+    ck[jf] = 1.0;
+    sk[jf] = 0.0;
+    for (j=1; j<k ; j++) {
+      ck[j] = (double)((ck[k])*c1 + (sk[k])*s1);
+      sk[j] = (double)((ck[k])*s1 - (sk[k])*c1);
+      k--;
+      ck[k] = ck[j];
+      sk[k] = -(sk[j]);
+    }
+ lbl250:
+    k1 = kk;
+    k2 = kk + kspnn;
+    aa = a[kk];
+    bb = b[kk];
+    ak = aa;
+    bk = bb;
+    j = 1;
+    k1 += kspan;
+    do {
+      k2 -= kspan;
+      j++;
+      at[j] = a[k1] + a[k2];
+      ak = at[j] + ak;
+      bt[j] = b[k1] + b[k2];
+      bk = bt[j] + bk;
+      j++;
+      at[j] = a[k1] - a[k2];
+      bt[j] = b[k1] - b[k2];
+      k1 += kspan;
+    } while (k1 < k2);
+    a[kk] = (double)ak;
+    b[kk] = (double)bk;
+    k1 = kk;
+    k2 = kk + kspnn;
+    j = 1;
+ lbl270:
+    k1 += kspan;
+    k2 -= kspan;
+    jj = j;
+    ak = aa;
+    bk = bb;
+    aj = 0.0;
+    bj = 0.0;
+    k = 1;
+    do {
+      k++;
+      ak = (at[k] * ck[jj]) + ak;
+      bk = (bt[k] * ck[jj]) + bk;
+      k++;
+      aj = (at[k] * sk[jj]) + aj;
+      bj = (bt[k] * sk[jj]) + bj;
+      jj += j;
+      if (jj > jf)
+        jj -= jf;
+    } while (k < jf);
+    k = jf - j;
+    a[k1] = (double)(ak - bj);
+    b[k1] = (double)(bk + aj);
+    a[k2] = (double)(ak + bj);
+    b[k2] = (double)(bk - aj);
+    j++;
+    if (j < k)     goto lbl270;
+    kk += kspnn;
+    if (kk <= nn)  goto lbl250;
+    kk -= nn;
+    if (kk<=kspan) goto lbl250;
+
+/*
+ * multiply by rotation factor (except for factors of 2 and 4)
+ */
+
+ lbl290:
+    if (i==m) goto lbl350;
+    kk = jc + 1;
+ lbl300:
+    c2 = 1.0 - cd;
+    s1 = sd;
+    mm = (kspan < klim ? kspan : klim);
+    goto lbl320;
+ lbl310:
+    c2 = c1 - ((cd*c1) + (sd*s1));
+    s1 = s1 + ((sd*c1) - (cd*s1));
+ lbl320:
+    c1 = c2;
+    s2 = s1;
+    kk += kspan;
+ lbl330:
+    ak = a[kk];
+    a[kk] = (double)((c2*ak) - (s2 * b[kk]));
+    b[kk] = (double)((s2*ak) + (c2 * b[kk]));
+    kk += kspnn;
+    if (kk <= nt) goto lbl330;
+    ak = s1*s2;
+    s2 = (s1*c2) + (c1*s2);
+    c2 = (c1*c2) - ak;
+    kk -= (nt - kspan);
+    if (kk <= kspnn) goto lbl330;
+    kk -= (kspnn - jc);
+    if (kk <= mm)   goto lbl310;
+    if (kk < kspan) goto lbl340;
+    kk -= (kspan - jc - inc);
+    if (kk <= (jc+jc)) goto lbl300;
+    goto lbl40;
+ lbl340:
+    s1 = ((double)((kk-1)/jc)) * dr * rad;
+    c2 = cos(s1);
+    s1 = sin(s1);
+    mm = (kspan < mm+klim ?  kspan :mm+klim);
+    goto lbl320;
+
+/*
+ * permute the results to normal order---done in two stages
+ * permutation for square factors of n
+ */
+
+ lbl350:
+    np[1] = ks;
+    if (!(*kt)) goto lbl440;
+    k = *kt + *kt + 1;
+    if (m < k)
+      k--;
+    np[k+1] = jc;
+    for (j=1; j < k; j++,k--) {
+      np[j+1] = np[j] / nfac[j];
+      np[k] = np[k+1] * nfac[j];
+    }
+    k3 = np[k+1];
+    kspan = np[2];
+    kk = jc + 1;
+    k2 = kspan + 1;
+    j = 1;
+    if (n != ntot) goto lbl400;
+/*
+ * permutation for single-variate transform (optional code)
+ */
+ lbl370:
+    do {
+      ak = a[kk];
+      a[kk] = a[k2];
+      a[k2] = (double)ak;
+      bk = b[kk];
+      b[kk] = b[k2];
+      b[k2] = (double)bk;
+      kk += inc;
+      k2 += kspan;
+    } while (k2 < ks);
+lbl380:
+    do {
+      k2 -= np[j++];
+      k2 += np[j+1];
+    } while (k2 > np[j]);
+    j = 1;
+ lbl390:
+    if (kk < k2) {
+      goto lbl370;
+    }
+    kk += inc;
+    k2 += kspan;
+    if (k2 < ks) goto lbl390;
+    if (kk < ks) goto lbl380;
+    jc = k3;
+    goto lbl440;
+/*
+ * permutation for multivariate transform
+ */
+ lbl400:
+    do {
+      do {
+        k = kk + jc;
+        do {
+          ak = a[kk];
+          a[kk] = a[k2];
+          a[k2] = (double)ak;
+          bk = b[kk];
+          b[kk] = b[k2];
+          b[k2] = (double)bk;
+          kk += inc;
+          k2 += inc;
+        } while (kk < k);
+        kk += (ks - jc);
+        k2 += (ks - jc);
+      } while (kk < nt);
+      k2 -= (nt - kspan);
+      kk -= (nt - jc);
+    } while (k2 < ks);
+ lbl420:
+    do {
+      k2 -= np[j++];
+      k2 += np[j+1];
+    } while (k2 > np[j]);
+    j = 1;
+ lbl430:
+    if (kk < k2)         goto lbl400;
+    kk += jc;
+    k2 += kspan;
+    if (k2 < ks)      goto lbl430;
+    if (kk < ks)      goto lbl420;
+    jc = k3;
+ lbl440:
+    if ((2*(*kt))+1 >= m)
+      return;
+
+    kspnn = *(np + *(kt) + 1);
+    j = m - *kt;
+    nfac[j+1] = 1;
+ lbl450:
+    nfac[j] = nfac[j] * nfac[j+1];
+    j--;
+    if (j != *kt) goto lbl450;
+    *kt = *(kt) + 1;
+    nn = nfac[*kt] - 1;
+    jj = 0;
+    j = 0;
+    goto lbl480;
+ lbl460:
+    jj -= k2;
+    k2 = kk;
+    kk = nfac[++k];
+ lbl470:
+    jj += kk;
+    if (jj >= k2) goto lbl460;
+    np[j] = jj;
+ lbl480:
+    k2 = nfac[*kt];
+    k = *kt + 1;
+    kk = nfac[k];
+    j++;
+    if (j <= nn) goto lbl470;
+/* Determine permutation cycles of length greater than 1 */
+    j = 0;
+    goto lbl500;
+ lbl490:
+    k = kk;
+    kk = np[k];
+    np[k] = -kk;
+    if (kk != j) goto lbl490;
+    k3 = kk;
+ lbl500:
+    kk = np[++j];
+    if (kk < 0)  goto lbl500;
+    if (kk != j) goto lbl490;
+    np[j] = -j;
+    if (j != nn) goto lbl500;
+    maxf *= inc;
+    /* Perform reordering following permutation cycles */
+    goto lbl570;
+ lbl510:
+    j--;
+    if (np[j] < 0) goto lbl510;
+    jj = jc;
+ lbl520:
+    kspan = jj;
+    if (jj > maxf)
+      kspan = maxf;
+    jj -= kspan;
+    k = np[j];
+    kk = (jc*k) + i + jj;
+    k1 = kk + kspan;
+    k2 = 0;
+ lbl530:
+    k2++;
+    at[k2] = a[k1];
+    bt[k2] = b[k1];
+    k1 -= inc;
+    if (k1 != kk) goto lbl530;
+ lbl540:
+    k1 = kk + kspan;
+    k2 = k1 - (jc * (k + np[k]));
+    k = -(np[k]);
+ lbl550:
+    a[k1] = a[k2];
+    b[k1] = b[k2];
+    k1 -= inc;
+    k2 -= inc;
+    if (k1 != kk) goto lbl550;
+    kk = k2;
+    if (k != j)   goto lbl540;
+    k1 = kk + kspan;
+    k2 = 0;
+ lbl560:
+    k2++;
+    a[k1] = at[k2];
+    b[k1] = bt[k2];
+    k1 -= inc;
+    if (k1 != kk) goto lbl560;
+    if (jj)       goto lbl520;
+    if (j  != 1)  goto lbl510;
+lbl570:
+    j = k3 + 1;
+    nt -= kspnn;
+    i = nt - inc + 1;
+    if (nt >= 0)  goto lbl510;
+    return;
+}
+
+
+/*
+ *-----------------------------------------------------------------------
+ * subroutine:
+  reals
+ * used with 'fft' to compute fourier transform or inverse for real data
+ *-----------------------------------------------------------------------
+ *      this is the call from C:
+
+ *              reals_(anal,banal,N2,mtwo);
+ *      which has been changed from CARL call
+ *              reals_(anal,banal,&N2,&mtwo);
+ */
+
+void
+reals_(double *a, double *b, int n, int isn)
+
+  /*    *a,       a refers to an array of floats 'anal'   */
+  /*    *b;       b refers to an array of floats 'banal'  */
+/* See IEEE book for a long comment here on usage */
+
+{
+    int inc,
+      j,
+      k,
+      lim,
+      mm,ml,
+      nf,nk,nh;
+
+    double      aa,ab,
+      ba,bb,
+      cd,cn,
+      dr,
+      em,
+      rad,re,
+      sd,sn;
+    double      xx;     /******* ADDED APRIL 1991 ******/
+    /* adjust  input array pointers (called from C) */
+    a--;        b--;
+    inc=abs(isn);
+    nf=abs(n);
+    if (nf*isn==0) {
+      printf("\nerror - zero in reals parameters : %d : %d ",n,isn);
+      return;
+    }
+    nk = (nf*inc) + 2;
+    nh = nk/2;
+/*****************************
+        rad  = atan((double)1.0);
+******************************/
+    rad = 0.785398163397448278900;
+    dr = -4.0/(double)(nf);
+/********************************** POW2 REMOVED APRIL 1991 *****************
+                                cd = 2.0 * (pow2(sin((double)0.5 * dr * rad)));
+*****************************************************************************/
+    xx = sin((double)0.5 * dr * rad);
+    cd = 2.0 * xx * xx;
+    sd = sin(dr * rad);
+/*
+ * sin,cos values are re-initialised each lim steps
+ */
+    lim = 32;
+    mm = lim;
+    ml = 0;
+    sn = 0.0;
+    if (isn<0) {
+      cn = 1.0;
+      a[nk-1] = a[1];
+      b[nk-1] = b[1]; }
+    else {
+      cn = -1.0;
+      sd = -sd;
+    }
+    for (j=1;j<=nh;j+=inc)      {
+      k = nk - j;
+      aa = a[j] + a[k];
+      ab = a[j] - a[k];
+      ba = b[j] + b[k];
+      bb = b[j] - b[k];
+      re = (cn*ba) + (sn*ab);
+      em = (sn*ba) - (cn*ab);
+      b[k] = (double)((em-bb)*0.5);
+      b[j] = (double)((em+bb)*0.5);
+      a[k] = (double)((aa-re)*0.5);
+      a[j] = (double)((aa+re)*0.5);
+      ml++;
+      if (ml!=mm) {
+        aa = cn - ((cd*cn)+(sd*sn));
+        sn = ((sd*cn) - (cd*sn)) + sn;
+        cn = aa;}
+      else {
+        mm +=lim;
+        sn = ((double)ml) * dr * rad;
+        cn = cos(sn);
+        if (isn>0)
+          cn = -cn;
+        sn = sin(sn);
+      }
+    }
+    return;
+}
+
+

dev/externals/fastconv/readme.txt

@@ -0,0 +1,110 @@
+Preliminary documentation  for fastconv.
+
+
+FASTCONV:  multi-channel fast convolution (using FFTs)
+
+version: 1.0 2010
+
+usage message: 
+fastconv [-aX][-f] infile impulsefile outfile [dry]
+   -aX        : scale output amplitude by X
+   -f         : write output as floats (no clipping)
+  infile      : input soundfile to be processed.
+  impulsefile : soundfile or text file containing impulse response,
+                  e.g. reverb or FIR filter.
+                (text file name must have extension .txt)
+              Supported channel combinations:
+               (a) mono infile, N-channel impulsefile;
+               (b) channels are the same;
+               (c) mono impulsefile, N-channel infile.
+  [dry]       :  set dry/wet mix (e.g. for reverb)
+                 Range: 0.0 - 1.0,  default = 0.0
+                 (uses sin/cos law for constant power mix)
+Note: some recorded reverb impulses effectively include the direct signal.
+In such cases  <dry>  need not be used
+Where impulsefile is a filter response (FIR), optimum length is power-of-two - 1.
+
+
+The primary application of fastconv is convolution reverberation using a 
+sampled impulse response of a building or other responsive space. The term "fast"  
+refers to the use of the Fast Fourier transform (FFT) to perfume the convolution. 
+The program can also be used more experimentally, as the impulse response input can 
+be any mono or multi-channel file (see details of available channel combinations below); 
+a file can also be convolved with itself. 
+The program uses double-precision processing throughout, and files of considerable 
+length can be processed cleanly.
+Note however that the FFT of the impulse response is stored in main memory, 
+so very large files may raise memory demands to critical levels.
+
+
+More on channel options:
+
+  Impulse soundfile can be multi-channel.
+  The .amb format is supported, also Logic Pro SDIR files for reading (change the file extension to .aifc).
+
+  When the infile is multi-channel, impulse file must be mono, or the same number of channels.
+        Where the impulse file is mono, data is duplicated for all input channels.
+		Typical usage: linear-phase filtering. Output should be 100% "wet".
+	The optimum length for a filter impulse response is power-of-two-1, 
+	e.g. 127,255, 511 etc. Most filter creation tools will output a file of this size.
+
+	(More tools to support linear-phase filtering are in preparation!).
+
+
+  When infile is mono, impulse soundfile can be multi-channel, outfile has channel count of
+      impulse response.
+    Typical usage: reverb convolution for spatialization, B-Format convolution
+	It will be usual to supply a non-zero value for "dry", e.g. 0.5. Note however that some recorded 
+       or synthetic impulse responses  may already include a "direct" component. In such cases, 
+	a "dry" value may not be needed.
+
+  The program employs an rms-based gain scaling algorithm which attempts to ensure all outputs are
+ approximately at the same level as the input. In normal use (e.g. a naturally decaying reverb impulse response),
+ the -a flag should not be needed.
+When the -f flag is used, output is forced to floats, with no clipping. The program reports 
+the output level together with a suggested corrective gain factor. This will be of particular relevance to more
+experimental procedures, such as convolving a soundfile with itself or with some other arbitrary source.
+
+
+Note for advanced users - use for FIR linear-phase filtering.
+
+Convolution implements a filtering process - equivalent to multiplication of the spectra of the two inputs. 
+The most common filter in audio work is recursive - it recycles output samples back into the input. The output 
+continues in principle for ever, hence the term "infinite impulse response" (IIR). The advantage of this technique 
+(as employed for example in CDP's "filter' package) is that even a low-order filter (i.e. using a small 
+number of delayed inputs and outputs) can be very powerful in its effects. A disadvantage in some applications 
+can be that an IIR filter changes the phase of components in the input (frequency components are 
+delayed by different amounts). This means among the things that the waveform shape of the input is not 
+preserved. The timbre of the sound (even in regions  not directly boosted 
+or attenuated by the filter) will therefore be changed. In common parlance, such a filter "colours" the signal.
+
+The alternative is a linear phase filter, which preserves all phase relationships in the signal. 
+All frequency components are delayed equally. To achieve this the impulse response must have a symmetrical 
+shape (see illustration). The response decays identically either side of the central peak.  
+
+This requires that there be no recirculated outputs reinjected into the filter. 
+Such a filter has a Finite Impulse Response (FIR). The impulse response data now comprises literally the 
+response of the filter to a single input sample (impulse). An impulse response of 31 samples means 
+that the filter generates and adds together 31 delayed copies of the input, to create each output sample. 
+While IIR filter coefficients may involve only two delayed samples ("a "second-order" filter), FIR responses 
+need to employ many more samples to achieve similar effects. It would not be unusual to use a 511-th order FIR filter. 
+This also means that the overall delay ("latency") of a FIR filter is much longer than that of an IIR filter. 
+
+A FIR filter cannot resonate as an IIR filter can. By computing only delayed inputs, It is unconditionally 
+"stable" - whereas a badly designed IIR filter can "blow up" with output values rapidly exceeding the sample limit.
+
+Fastconv supports the use of FIR coefficient files either in the form of either a short soundfile, or a 
+plain text file containing (in a single column)  the list of coefficients as floating-point numbers within 
+the "normalised" range -1.0 to 1.0. For orthodox filtering purposes a mono soundfile should be used, to process 
+all channels identically. FIR coefficient text files are generated by many engineering-oriented filter design 
+applications.  User may be tempted to  write response files by hand; this can be done, but the results will be 
+virtually impossible to predict or control. 
+
+For maximum efficiency, such files should ideally have a size that is a power-of-two less one: e.g. 255, 511, 1023, etc.
+
+For more information about FIR filters, see :
+
+http://www.labbookpages.co.uk/audio/firWindowing.html
+
+
+August 23 2010