ref: ca201805c6603ddcb2617976a40566f4b58703cb
parent: 19f64639e51ac6436c0487e8382dfd778a7d6d60
author: cbagwell <cbagwell>
date: Thu Nov 11 15:58:55 EST 1999
Adding new test code.
--- /dev/null
+++ b/test/Makefile
@@ -1,0 +1,26 @@
+# makefile for gcc
+
+# Uncomment the following if you need a specific compiler.
+#CC=gcc
+
+CFLAGS:=-O2 -Wall -L../
+LDFLAGS:=-s
+LIBS:=-lm
+
+all: ding model lding lmodel
+
+ding: ding.c
+ $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $< $(LIBS)
+
+model: model.c
+ $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $< $(LIBS)
+
+lding: ding.c
+ $(CC) $(CFLAGS) -DLSAMPL $(LDFLAGS) -o $@ $< $(LIBS)
+
+lmodel: model.c
+ $(CC) $(CFLAGS) -DLSAMPL $(LDFLAGS) -o $@ $< $(LIBS)
+
+clean:
+ rm -f ding lding model lmodel
+
--- /dev/null
+++ b/test/README
@@ -1,0 +1,275 @@
+Here are some helper programs for checking filter and
+rate conversion accuracy of sox. The programs in the
+test subdirectory which I (Stanley J. Brooks) wrote are
+covered by the GPL. See the Copyright file for details.
+Sox itself is not covered by the GPL, but something similar.
+
+To use this stuff, first cd to this test subdirectory
+and run make, which should produce the ding,lding,
+model, and lmodel executables.
+
+Now to test something...
+
+Say you want to compare the 'band' effect to the 'filter' effect...
+
+./ltest band 800 200 >A
+./ltest filter 600-1000 >B
+gnuplot responseAB
+
+plots the response curves with y being power-gain in dB,
+and x the frequency.
+
+For another, say to compare response and error of resample
+rate conversion from 8000 ->22050 samples per second,
+with linear (default) interpolation versus -qs quadratic
+interpolation with Nuttall window...
+
+./ltest -l resample >A
+./ltest -l resample -qs 0.80 0 >B
+gnuplot plotAB
+
+-------------------------------------------------------------------
+
+About ltest:
+
+-l means use 32-bit signed samples, otherwise 16-bit signed is used.
+
+The other parameters are fed into sox as the 'effect' with parameters.
+
+The rates 8000:22050 are in the script, but you can edit a perl script,
+I hope. just change the ($rate0,$rate1)=(8000,22050) line near the top.
+
+What ltest does:
+
+It uses the 'ding' program to synthesize a sineusoid input file for the
+frequencies 0.00 0.01 ... 0.99 of the Nyquist frequency. This input
+has:
+ 400 samples samples of silence,
+ 4000 samples with smooth envelope rising to volume -v0.5
+ 16000 samples at -v0.5
+ 4000 samples with smooth envelope falling to 0,
+ 400 sample more of silence.
+
+The rising/falling envelopes are shaped like rising, falling portions
+of the (1-cos(x)) function.
+
+Then the filter or rate-change effect is applied to this i0.xx.xx file
+to give an output file.
+
+The output file is examined by the 'model' or 'lmodel' program to analyse
+response level and error level.
+
+model works as follows:
+
+step 1: the entire sample file is read in, and the center-of-gravity of
+ the squared samples is found. This is the time offset which should
+ correspond to the center of the filtered/resampled tone-pulse.
+
+step 2: let N = 16000*(rate1/rate0) be the number of samples at output rate
+ which would correspond to the 16000 sample duration at max volume.
+ We focus on the samples between
+ (center - 0.3*N) and (center + 0.9*N)
+ where the transient effects of attack/release envelope should be
+ small.
+ We do a least-squares fit of a sinusoid at the adjusted frequency
+ to this segment of 0.6*N samples, and print out the component
+ s2max which can be explained, and also the rms level of the 'error'
+ or unexplained part.
+
+That's pretty much it... the ltest perl-script glues it together and outputs
+adjusted data which the gnuplot will like.
+
+Modify the perl and gnuplot scripts to suit your needs.
+
+ SoX: Sound eXchange
+
+
+SoX (also known as Sound eXchange) translates sound samples between different
+file formats, and optionally performs various sound effects.
+
+This release understands:
+
+ o Raw files in various binary formats
+ o Raw textual data
+ o Microsoft .WAV files
+ o PCM, u-law, a-law
+ o MS ADPCM (Read only)
+ o IMA ADPCM (Read only)
+ o MAUD files
+ o Sound Blaster .VOC files
+ o IRCAM SoundFile files
+ o SUN .au files
+ o PCM, u-law, a-law
+ o G7xx ADPCM files (read only)
+ o mutant DEC .au files
+ o Apple/SGI AIFF files
+ o CD-R (music CD format)
+ o Macintosh HCOM files
+ o Sounder files
+ o NeXT .snd files
+ o Soundtool (DOS) files
+ o Psion (palmtop) A-law files
+ o AVR files
+
+The sound effects include:
+
+ o Channel Averaging
+ o Band-pass filter
+ o Chorus effect
+ o Cut out loop samples
+ o Add an echo
+ o Add a sequence of echos
+ o Apply a flanger effect
+ o Apply a high-pass filter
+ o Apply a low-pass filter
+ o Display a list of loops in a file
+ o Add masking noise to a signal
+ o Apply a phaser effect
+ o Convert from stereo to mono
+ o Change sampling rates using several different algorithms.
+ o Apply a reverb effect
+ o Reverse the sound samples (to search for Satanic messages ;-)
+ o Convert from mono to stereo
+ o Swap stereo channels
+ o Display general stats on a sound sample
+ o Add the world-famous Fender Vibro-Champ effect
+
+Big news! Lots of new effects have been added. This includes most the
+popular "Guitar Effects" talked about in the same named FAQ available.
+
+The 'resample' and 'polyphase' effect does high-grade signal rate
+changes using real signal theory. Yes, it's very slow.
+
+History:
+
+This is the 12th release, Patchlevel 17 of the Sound Tools.
+SoX was originally written and maintained by Lance Norskog but
+unfortunetly he has stopped maintaining it since 1995. I, Chris
+Bagwell (cbagwell@sprynet.com), have started maintaining it since
+1996 to the present.
+
+Caveats:
+
+SoX is intended as the Swiss Army knife of sound processing tools. It
+doesn't do anything very well, but sooner or later it comes in very handy.
+SoX is really only usable day-to-day if you hide the wacky options with
+one-line shell scripts.
+
+Installing:
+
+Unless your using a precompiled binary version, you will need to compile
+SoX as described in the INSTALL file. Please read that file for further
+instructions.
+
+Now, read TIPS, CHEAT.eft and CHEAT. These give a background on how
+SoX deals with sound files and how to convert this format
+to that format, and apply various effects.
+
+SoX uses file suffices to determine the nature of a sound sample file.
+If it finds the suffix in its list, it uses the appropriate read
+or write handler to deal with that file. You may override the suffix
+by giving a different type via the '-t type' argument. See the manual
+page for more information.
+
+SoX has an auto-detect feature that attempts to figure out
+the nature of an unmarked sound sample. It works very well.
+This feature is used if you specify '-t auto' for the file type.
+
+I hope to inspire the creation of a common base of sound processing
+tools for computer multimedia work, similar to the PBM toolkit for
+image manipulation.
+
+Sound Tools may be used for any purpose. Source distributions must
+must include the copyright notices, and (lack of) warranty information.
+Binary distributions must include acknowledgements to the creators.
+Files are copyright by their respective authors.
+
+If you have bug fixes/enhancements, please send it to me as I would like
+to coordinate the releases. Please document your changes. I don't
+possess every kind of computer currently sold, and SoX is now beyond
+the phase where I can understand and test most of your contributions.
+
+The majority of SoX features and source code are contributed
+by you the user. Thank you very much for making SoX a success!
+
+ Creator:
+ Lance Norskog thinman@meer.net (inactive currently)
+
+ Mantainer:
+ Chris Bagwell cbagwell@sprynet.com
+
+ Contributors:
+ Juergen Mueller jmueller@uia.ua.ac.be
+ chorus, echo, echos, flanger, phaser, and reverb
+ effects.
+ Guido Van Rossum guido@cwi.nl
+ AU, AIFF, AUTO, HCOM, reverse,
+ many bug fixes
+ Jef Poskanzer jef@well.sf.ca.us
+ original code for u-law and delay line
+ Bill Neisius bill%solaria@hac2arpa.hac.com
+ DOS port, 8SVX, Sounder, Soundtool formats
+ Apollo fixes, stat with auto-picker
+ Rick Richardson rick@digibd.com
+ WAV and SB driver handlers, fixes
+ David Champion dgc3@midway.uchicago.edu
+ Amiga port
+ Pace Willisson pace@blitz.com
+ Fixes for ESIX
+ Leigh Smith leigh@psychokiller.dialix.oz.au
+ SMP and comment movement support.
+ AIFF Loop/MIDI support
+ David Sanderson dws@ssec.wisc.edu
+ AIX3.1 fixes
+ Glenn Lewis glewis@pcocd2.intel.com
+ AIFF chunking fixes
+ Brian Campbell brianc@quantum.qnx.com
+ QNX port and 16-bit fixes
+ Chris Adams gt8741@prism.gatech.edu
+ DOS port fixes
+ John Kohl jtkohl@kolvir.elcr.ca.us
+ BSD386 port, VOC stereo support
+ Ken Kubo ken@hmcvax.claremont.edu
+ VMS port, VOC stereo support
+ Frank Gadegast <phade@cs.tu-berlin.de>
+ Microsoft C 7.0 & C Borland 3.0 ports
+ David Elliot <dce@scmc.sony.com>
+ CD-R format support
+ David Sears <dns@essnj3.essnjay.com>
+ Linux support
+ Tom Littlejohn <tlit@seq1.loc.gov>
+ Raw textual data
+ Boisy G. Pitre boisy@microware.com
+ OS9 port
+ Sun Microsystems, Guido Van Rossum
+ CCITT G.711, G.721, G.723 implementation
+ Graeme Gill graeme@labtam.labtam.oz.au
+ A-LAW format, Good .WAV handling,
+ avg channel expansion
+ Allen Grider grider@hfsi.hfsi.com
+ VOC stereo mode, WAV file handling
+ Michel Fingerhut Michel.Fingerhut@ircam.fr
+ Upgrade 'sf' format to current IRCAM format.
+ Float file support.
+ Chris Knight
+ Achimedes Acorn support
+ Richard Caley R.Caley@ed.ac.uk
+ Psion WVE handler
+ Lutz Vieweg lkv@mania.RoBIN.de
+ MAUD (Amiga) file handler
+ Tim Gardner timg@tpi.com
+ Windows NT port for V7
+ Jimen Ching jiching@wiliki.eng.hawaii.edu
+ Libst porting bugs
+ Lauren Weinstein lauren@vortex.com
+ DOS porting, scripts, professional use
+ Stan Brooks stabro@megsinet.net
+ Rewrite of resample and polyphase code.
+ DSP filter effect. Some test code/scripts.
+ Chris Bagwell cbagwell@sprynet.com
+ OSS and Sun players, bugfixes, ADPCM support,
+ patch collection and maintance.
+ (your name could be here, too)
+ (I've probably lost a few, and several people fixed
+ the same bugs.)
+
--- /dev/null
+++ b/test/ding.c
@@ -1,0 +1,287 @@
+/*
+ ding.c -- program to generate testpatterns for DSP code.
+
+ Copyright (C) 1999 Stanley J. Brooks <stabro@megsinet.net>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program 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 General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+*/
+
+#include <string.h>
+// for open,read,write:
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <errno.h>
+#include <unistd.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define FLOAT double
+
+#ifndef LSAMPL
+# define SAMPL short
+# define MAXSAMPL 0x7fff
+# define MINSAMPL -0x7fff
+#else
+# define SAMPL long
+# define MAXSAMPL 0x7fffffff
+# define MINSAMPL -0x7fffffff
+#endif
+
+static void Usage(void)__attribute__((noreturn));
+
+static void Usage(void)
+{+ fprintf(stderr, "Usage: ./ding [options] [<in-file>] <out-file>\n");
+ fprintf(stderr, " Options:\n");
+ fprintf(stderr, " [-v <vol>] float, volume, 1.00 is max\n");
+ fprintf(stderr, " [-f <freq>] float, frequency = freq*nyquist_rate\n");
+ fprintf(stderr, " [-d <decay>] float, per-sample decay factor\n");
+ fprintf(stderr, " [-o <n>] int, zero-pad samples before tone\n");
+ fprintf(stderr, " [-e <n>] int, length in samples of tone\n");
+ fprintf(stderr, " [-p <n>] int, zero-pad samples after tone\n");
+ exit(-1);
+}
+
+static int ReadN(int fd, SAMPL *v, int n)
+{+ int r,m;
+ static SAMPL frag=0;
+ static int fraglen=0;
+ char *q;
+
+ q=(char*)v;
+ if (fraglen)
+ memcpy(q, (char*)&frag, fraglen);
+
+ m = n*sizeof(SAMPL) - fraglen;
+ q += fraglen;
+ if (fd>=0) {+ do {+ r = read(fd, q, m);
+ }while(r==-1 && errno==EINTR);
+ if (r==-1) {+ perror("Error reading fd1"); exit(-1);+ }
+ }else{+ bzero(q,m);
+ r = m;
+ }
+ r += fraglen;
+ fraglen = r%sizeof(SAMPL);
+ if (fraglen)
+ memcpy((char*)&frag, (char*)v, fraglen);
+
+ return (r/sizeof(SAMPL));
+}
+
+#define BSIZ 0x10000
+
+int main(int argct, char **argv)
+{+ int optc;
+ int fd1,fd2;
+ char *fnam1,*fnam2;
+ int len, st;
+ SAMPL *ibuff,max,min;
+ int poflo,moflo;
+ FLOAT Vol0=1;
+ FLOAT Freq=0; /* of nyquist */
+ int Offset=0;
+ int Pad=0;
+ struct _Env {+ int r; /* rise */
+ int c; /* center */
+ int f; /* fall */
+ FLOAT v; /* volume */
+ FLOAT d; /* decay */
+ } Env = {0,0,0,0.5,1.0};+
+ double x=1,y=0;
+ double thx=1,thy=0;
+
+ static inline void a_init(double frq)
+ {+ if (frq != 0) {+ x = 1; y = 0;
+ } else {+ x = 0; y = 1;
+ }
+ thx = cos(frq*M_PI);
+ thy = sin(frq*M_PI);
+ }
+
+ static inline const double a(int k, int L)
+ {+ double a, u;
+ u = (double)k/L; /* in [0,1] */
+ a = 0.5*(1 + cos(M_PI * u));
+ //printf("a(%d,%d) = %.5f, l=%.1f, u=%.6f\n",k,L,a,l,u); + return a;
+ }
+
+ static inline void a_post(int k)
+ {+ double x1;
+ x1 = x*thx - y*thy;
+ y = x*thy + y*thx;
+ x = x1;
+ /* update (x,y) for next tic */
+ if (!(k&0x0f)) { /* norm correction each 16 samples */+ x1 = 1/sqrt(x*x+y*y);
+ x *= x1;
+ y *= x1;
+ }
+ }
+
+ /* Parse the options */
+ while ((optc = getopt(argct, argv, "d:o:e:t:p:f:v:h")) != -1) {+ char *p;
+ switch(optc) {+ case 'd':
+ Env.d = strtod(optarg,&p);
+ if (p==optarg || *p) {+ fprintf(stderr,"option -%c expects float value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'f':
+ Freq = strtod(optarg,&p);
+ if (p==optarg || *p) {+ fprintf(stderr,"option -%c expects float value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'e':
+ p = optarg;
+ Env.c = strtol(p,&p,10);
+ if (*p++ == ':') {+ Env.r = Env.f = Env.c;
+ Env.c = strtol(p,&p,10);
+ }
+ if (*p++ == ':') {+ Env.f = strtol(p,&p,10);
+ }
+ if (*p || Env.c<=0 || Env.r<0 || Env.f<0) {+ fprintf(stderr,"option -%c not valid (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'o':
+ Offset = strtol(optarg,&p,10);
+ if (p==optarg || *p || Offset<0) {+ fprintf(stderr,"option -%c expects int value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'p':
+ Pad = strtol(optarg,&p,10);
+ if (p==optarg || *p || Pad<0) {+ fprintf(stderr,"option -%c expects int value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'v':
+ Env.v = strtod(optarg,&p);
+ if (p==optarg || *p) {+ fprintf(stderr,"option -%c expects float value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'h':
+ default:
+ Usage();
+ }
+ }
+
+ Env.v *= MAXSAMPL;
+ if (Freq==0.0) Env.v *= sqrt(0.5);
+ //fprintf(stderr,"Vol0 %8.3f\n", Vol0);
+
+ len = Offset+Env.r+Env.c+Env.f+Pad;
+ fnam1=NULL; fd1=-1;
+ //fprintf(stderr,"optind=%d argct=%d\n",optind,argct);
+ if (optind <= argct-2) {+ fnam1=argv[optind++];
+ fd1=open(fnam1,O_RDONLY);
+ if (fd1<0) {+ fprintf(stderr,"Open: %s %s\n",fnam1,strerror(errno)); return(1);
+ }
+ len=lseek(fd1,0,SEEK_END)/2;
+ lseek(fd1,0,SEEK_SET);
+ }
+
+ if (optind != argct-1) Usage();
+
+ fd2 = 1; /* stdout */
+ fnam2=argv[optind++];
+ if (strcmp(fnam2,"-")) {+ fd2=open(fnam2,O_WRONLY|O_CREAT|O_TRUNC,0644);
+ if (fd2<0) {+ fprintf(stderr,"Open: %s %s\n",fnam2,strerror(errno)); return(1);
+ }
+ }
+ //fprintf(stderr, "Files: %s %s\n",fnam1,fnam2);
+
+ a_init(Freq);
+
+ ibuff=(SAMPL*)malloc(BSIZ*sizeof(SAMPL));
+ poflo=moflo=0;
+ max =MINSAMPL; min=MAXSAMPL;
+ for(st=0; st<len; ){+ int ct;
+ SAMPL *ibp;
+
+ ct = len - st; if (ct>BSIZ) ct=BSIZ;
+ ReadN(fd1,ibuff,ct);
+ for (ibp=ibuff; ibp<ibuff+ct; ibp++,st++) {+ int k;
+ double v = *ibp;
+
+ if (max<*ibp) max=*ibp;
+ else if (min>*ibp) min=*ibp;
+ v *= Vol0;
+
+ k = st-Offset;
+ if (k>=0 && k<Env.r) {+ v += y*a(Env.r-k,Env.r)*Env.v;
+ a_post(st);
+ }else if (k-=Env.r, k>=0 && k<Env.c) {+ v += y*Env.v;
+ a_post(st);
+ }else if (k-=Env.c, k>=0 && k<Env.f) {+ v += y*a(k,Env.f)*Env.v;
+ a_post(st);
+ }
+
+ if (v>MAXSAMPL) {+ poflo++;
+ v=MAXSAMPL;
+ } else if (v<MINSAMPL) {+ moflo++;
+ v=MINSAMPL;
+ }
+ *ibp = rint(v);
+ }
+
+ write(fd2,(char*)ibuff,ct*sizeof(SAMPL));
+ }
+
+ fprintf(stderr,"input range: [%ld,%ld] pos/neg oflos: %d/%d\n",min,max,poflo,moflo);
+ return 0;
+
+}
--- /dev/null
+++ b/test/ltest.pl
@@ -1,0 +1,106 @@
+#!/usr/bin/perl -w
+
+# ltest.pl -- program to help generate response/error graphs for DSP code.
+#
+# Copyright (C) 1999 Stanley J. Brooks <stabro@megsinet.net>
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program 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 General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+use strict;
+$|=1;
+
+# set default values
+my $sox='../sox';
+my $effect='resample';
+#my $effect='resample -qs';
+#my $effect='resample -q';
+#my $effect='resample -ql';
+#my $effect='polyphase -cutoff 0.90';
+#my $effect='filter 400-2000';
+#my $effect='filter 400-2000 1024';
+my ($rate0,$rate1)=(8000,22050); # sample rates
+my $p=400; # silence before/after tonepulse
+my $env="-e4000:16000:4000"; # attack, duration, drop
+
+# parse commandline arguments
+my $updown = 0; # set to 1 for up/down rate-conversion test
+my ($ding,$model,$t,$rms,$lim)=("./ding","./model","sw",11585.2, 0.5);+while ($_ = shift @ARGV) {+ if ($_ eq "-l") {+ ($ding,$model,$t,$rms,$lim)=("./lding","./lmodel","sl",759250125.0, 50.0);+ } elsif ($_ eq "-ud") {+ $updown=1;
+ } else {+ unshift @ARGV,$_;
+ last;
+ }
+}
+if ($#ARGV >= 0) {+ $effect="@ARGV";
+}
+
+my $ratechange=0;
+if ($effect =~ m{^(resample|polyphase|rate)}) {+ $ratechange=1;
+}
+
+# output a nice header explaining the data
+if ($ratechange==0) {+ print("# Testing $sox -r$rate0 d/i0.xx.$t d/j0.xx.$t $effect\n");+} else {+ print("# Testing $sox -r$rate0 d/i*.$t -r$rate1 d/u*.$t $effect\n");+ if ($updown==1) {+ print("# then back down to $rate0\n");+ }
+}
+print("# with tone pulses from 0.00 to 0.99 percent of Nyquist\n");+print("# produced by $ding -f0.xx -v0.5 -o$p $env -p$p d/i0.xx.$t\n");+
+# generate the test data
+mkdir("d",0775);+my $f;
+my %q;
+for ($f=0.00; $f<1.0; $f+=0.01) {+ my @mod;
+ my $s=sprintf("%4.2f",$f);+ #print "$ding -v0.5 -o$p $env -p$p -f$s -d1.0 d/i$s.$t\n";
+ qx{$ding -v0.5 -o$p $env -p$p -f$s -d1.0 d/i$s.$t &>/dev/null};+ if ($ratechange==0) {+ qx{$sox -r$rate0 d/i$s.$t -r$rate0 d/j$s.$t $effect 2>/dev/null};+ @mod = grep {/v2max/} qx{$model -f$s $env $rate0 d/j$s.$t 2>&1};+ } else {+ qx{$sox -r$rate0 d/i$s.$t -r$rate1 d/u$s.$t $effect 2>/dev/null};+ if ($updown) {+ qx{$sox -r$rate1 d/u$s.$t -r$rate0 d/o$s.$t $effect 2>/dev/null};+ @mod = grep {/v2max/} qx{$model -f$s $env $rate0:$rate0 d/o$s.$t 2>&1};+ }else{+ @mod = grep {/v2max/} qx{$model -f$s $env $rate0:$rate1 d/u$s.$t 2>&1};+ }
+ }
+ print STDERR "$s: ",@mod;
+ $_=shift(@mod);
+ if (m{s2max *([0-9.]+), *v2max *([0-9.]+), *rmserr *(-?[0-9.]+)}) {+ #print("$s $1\n");+ #print("$s $1 $3\n");+ my $v = ($1 > $lim)? $1 : $lim;
+ my $r = ($3 > $lim)? $3 : $lim;
+ my $dbv = 20*log($v/$rms)/log(10);
+ my $dbr = 20*log($r/$rms)/log(10);
+ printf("%s %.3f %.3f\n",$s,$dbv,$dbr);+ }
+ #unlink "d/i$s.$t","d/u$s.$t","d/o$s.$t";
+}
+
+exit 0;
--- /dev/null
+++ b/test/model.c
@@ -1,0 +1,319 @@
+/*
+ model.c -- program to help test DSP filter and rate-conversion code.
+
+ Copyright (C) 1999 Stanley J. Brooks <stabro@megsinet.net>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program 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 General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+*/
+
+#include <string.h>
+// for open,read,write:
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <errno.h>
+#include <unistd.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define FLOAT double
+
+#ifndef LSAMPL
+# define SAMPL short
+# define MAXSAMPL 0x7fff
+# define MINSAMPL -0x8000
+#else
+# define SAMPL long
+# define MAXSAMPL 0x7fffff
+# define MINSAMPL -0x800000
+#endif
+
+static struct _Env {+ int r; /* rise */
+ int c; /* center */
+ int f; /* fall */
+ double v; /* volume */
+ double d; /* decay */
+} Env = {0,0,0,0.5,1.0};+
+static void Usage(void)__attribute__((noreturn));
+
+static void Usage(void)
+{+ fprintf(stderr,"Usage: ./model rate0,rate1 [-l len] [-f freq] <in-file>\n"); exit(-1);
+}
+
+static int ReadN(int fd, SAMPL *v, int n)
+{+ int r;
+
+ do {+ r = read(fd, (char*)(v), n*sizeof(SAMPL));
+ }while(r==-1 && errno==EINTR);
+ if (r==-1 || r%sizeof(SAMPL)) {+ perror("Error reading input"); exit(-1);+ }
+ return r/sizeof(SAMPL);
+}
+
+#define BSIZ 400000/*0x10000*/
+
+static double bestfit(double sx, double sy,
+ double h11, double h12, double h22,
+ double *cx, double *cy)
+{+ double a,su,uu,cu,w2;
+
+ a=0; *cy=0;
+ if (h22>1e-20*h11) {+ a = h12/h22; /* u[] = x[] - a*y[] is orthogonal to y[] */
+ *cy = sy/h22; /* *cy*y[] is projection of s[] onto y[] */
+ }
+ /* <s,x-ay> = sx -a*sy */
+ su = sx - a*sy; /* su is iprod of s[] and u[] */
+ /* <u,u> = <x-ay,x-ay> = h11 - 2a*h11 + a*a*h22 */
+ uu = h11 - 2*a*h12 + a*a*h22;
+ cu = 0;
+ if (uu>1e-20*h11)
+ cu = su/uu; /* s[] = *cy * y[] + cu * u[] is orthogonal decomposition */
+ w2 = *cy * *cy * h22;
+ w2 += cu * cu * uu;
+ *cx = cu;
+ *cy -= a*cu;
+ return w2;
+}
+
+static double eCenter(const SAMPL *ibuff, int len)
+{+ double v0,v1,v2;
+ int n;
+
+ v0 = v1 = v2 = 0;
+ for (n=0; n<len; n++) {+ double w = ibuff[n];
+ w = w*w;
+ v0 += w;
+ v1 += n*w;
+ v2 += (n*n)*w;
+ }
+ v1 /= v0;
+ v2 /= v0;
+ v2 -= v1*v1;
+ fprintf(stderr,"eCenter %.2f, STD %.2f\n",v1,sqrt(v2));
+ return v1;
+}
+
+static void
+bigcalc(double Factor, double Freq1, const SAMPL *ibuff, int len)
+{+ double c,del;
+ double x,y;
+ double thx,thy;
+ double Len1;
+ int k1,n;
+ long double h11,h12,h22;
+ long double sx,sy,ss;
+ double cn,cx,cy;
+ double s2=0, v2=0;
+ const SAMPL *ip;
+
+ static inline void a_init(double frq)
+ {+ x = 1;
+ y = 0;
+ thx = cos(frq*M_PI);
+ thy = sin(frq*M_PI);
+ }
+
+ static inline const double a(double k, double L)
+ {+ double a, l, u;
+ l = L/2.0;
+ u = k/l - 1; /* in (-1,+1) */
+ a = 0.5*(1 + cos(M_PI * u));
+ return a;
+ }
+
+ static inline void a_post(int k)
+ {+ double x1;
+ x1 = x*thx - y*thy;
+ y = x*thy + y*thx;
+ x = x1;
+ /* update (x,y) for next tic */
+ if ((k&0x0f)==0x0f) { /* norm correction each 16 samples */+ x1 = 1/sqrt(x*x+y*y);
+ x *= x1;
+ y *= x1;
+ }
+ }
+
+ c = eCenter(ibuff,len);
+ Len1 = Env.c*Factor*0.6;
+ c += Env.c*Factor*0.15;
+ {+ double a;
+ int b;
+ a = c-Len1/2;
+ b = ceil(a-0.001);
+ ip = ibuff + b;
+ del = b-a;
+ }
+ fprintf(stderr,"del %.3f\n", del);
+ k1 = Len1-del;
+
+ a_init(Freq1);
+ h11 = h12 = h22 = 0;
+ sx = sy = ss = 0;
+ for(n=0; n<=k1; n++) {+ double f,s;
+ double u,v;
+
+ s = ip[n]; /* sigval at n */
+
+ //if (fabs(n-Len1/2)<Len1/3) {+ f = 1; //a(n+del,Len1);
+ u = f*x; v = f*y;
+ sx += s*u;
+ sy += s*v;
+ ss += s*s;
+ h11 += u*u;
+ h12 += u*v;
+ h22 += v*v;
+ //}
+ a_post(n);
+ }
+ //fprintf(stderr,"h12 %.10f\n", (double)h12/(sqrt(h11)*sqrt(h22)));
+
+ v2 = ss/(k1+1);
+ s2 = bestfit(sx,sy,h11,h12,h22,&cx,&cy)/(k1+1);
+ cn = sqrt(cx*cx+cy*cy);
+ fprintf(stderr,"amp %.1f, cx %.10f, cy %.10f\n", cn, cx/cn, cy/cn);
+
+ fprintf(stderr,"[%.1f,%.1f] s2max %.2f, v2max %.2f, rmserr %.2f\n",
+ c/Factor, c, sqrt(s2), sqrt(v2), (s2<=v2)? sqrt(v2-s2):-sqrt(s2-v2));
+
+}
+
+int main(int argct, char **argv)
+{+ int optc;
+ int fd1;
+ int len;
+ SAMPL *ibuff;
+ int rate0, rate1;
+ double Freq0=0, Freq1; /* of nyquist */
+ double Factor;
+ char *fnam1;
+
+ /* Parse the options */
+ while ((optc = getopt(argct, argv, "d:o:e:p:f:h")) != -1) {+ char *p;
+ switch(optc) {+ case 'd':
+ Env.d = strtod(optarg,&p);
+ if (p==optarg || *p) {+ fprintf(stderr,"option -%c expects float value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'f':
+ Freq0 = strtod(optarg,&p);
+ if (p==optarg || *p) {+ fprintf(stderr,"option -%c expects float value (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'e':
+ p = optarg;
+ Env.c = strtol(p,&p,10);
+ if (*p++ == ':') {+ Env.r = Env.f = Env.c;
+ Env.c = strtol(p,&p,10);
+ }
+ if (*p++ == ':') {+ Env.f = strtol(p,&p,10);
+ }
+ if (*p || Env.c<=0 || Env.r<0 || Env.f<0) {+ fprintf(stderr,"option -%c not valid (%s)\n", optc, optarg);
+ Usage();
+ }
+ break;
+ case 'h':
+ default:
+ Usage();
+ }
+ }
+
+ //fprintf(stderr,"optind=%d argct=%d\n",optind,argct);
+
+ if (optind != argct-2) Usage();
+
+ {+ char *p0, *p;
+ p0 = argv[optind];
+ rate0 = rate1 = strtol(p0,&p,10);
+ if (*p) {+ if (*p != ':') {+ fprintf(stderr,"Invalid rate parameter (%s)\n", p0);
+ Usage();
+ }
+ p0 = p+1;
+ rate1 = strtol(p0,&p,10);
+ if (*p) {+ fprintf(stderr,"Invalid 2nd rate parameter (%s)\n", p0);
+ Usage();
+ }
+ }
+ if (rate0<=0 || rate1<=0) {+ fprintf(stderr,"Invalid rate parameter (%s)\n", argv[optind]);
+ Usage();
+ }
+ optind++;
+ }
+
+ Factor = (double)rate1 / (double)rate0;
+ Freq1 = Freq0/Factor;
+
+ //if (optind != argct-1) Usage();
+
+ fnam1=NULL; fd1=-1;
+ fnam1=argv[optind++];
+ fd1=open(fnam1,O_RDONLY);
+ if (fd1<0) {+ fprintf(stderr,"Open: %s %s\n",fnam1,strerror(errno)); return(1);
+ }
+
+ //fprintf(stderr, "Files: %s %s\n",fnam1,fnam2);
+
+ ibuff=(SAMPL*)malloc(BSIZ*sizeof(SAMPL));
+ {+ int n,r;
+
+ for (n=0; n<BSIZ; n+=r) {+ r = ReadN(fd1,ibuff+n,BSIZ-n);
+ if (r==0) break;
+ }
+ len = n;
+ }
+ fprintf(stderr,"Read %d input samples\n",len);
+
+ bigcalc(Factor, Freq1, ibuff, len);
+
+ return 0;
+
+}
--- /dev/null
+++ b/test/plotA
@@ -1,0 +1,16 @@
+#
+set autoscale
+#set logscale x
+set mxtics 5
+set mytics 2
+set grid xtics ytics
+set grid mxtics mytics
+show grid
+plot [0:4000] 'A' using (4000*$1):2 title "response" with lines,\
+ 'A' using (4000*$1):3 title "error" with lines
+pause -1 "Hit return to continue"
+quit
+
+#set key outside below
+#set key left box
+#reset
--- /dev/null
+++ b/test/plotAB
@@ -1,0 +1,18 @@
+#
+set autoscale
+#set logscale x
+set mxtics 5
+set mytics 2
+set grid xtics ytics
+set grid mxtics mytics
+show grid
+plot [0:4000] 'A' using (4000*$1):2 title "A response" with lines,\
+ 'A' using (4000*$1):3 title "A error" with lines,\
+ 'B' using (4000*$1):2 title "B response" with lines,\
+ 'B' using (4000*$1):3 title "B error" with lines
+pause -1 "Hit return to continue"
+quit
+
+#set key outside below
+#set key left box
+#reset
--- /dev/null
+++ b/test/responseAB
@@ -1,0 +1,16 @@
+#
+set autoscale
+#set logscale x
+set mxtics 5
+set mytics 2
+set grid xtics ytics
+set grid mxtics mytics
+show grid
+plot [0:4000] 'A' using (4000*$1):2 title "A response" with lines,\
+ 'B' using (4000*$1):3 title "B response" with lines
+pause -1 "Hit return to continue"
+quit
+
+#set key outside below
+#set key left box
+#reset