Skip site navigation (1)Skip section navigation (2)
Date:      Thu, 14 Nov 1996 02:27:23 -0500 (EST)
From:      "John S. Dyson" <toor@dyson.iquest.net>
To:        multimedia@freebsd.org
Subject:   I have been having fun with a new sw toy!!!
Message-ID:  <199611140727.CAA00238@dyson.iquest.net>

next in thread | raw e-mail | index | archive | help
Since there have been some questions about cdd, I thought that I might
post the latest copy of my audio compressor software.  Note that this
code is *NOT* size compression, it is level compression.  It takes
a 16-bit stereo file and either plays the file onto your sound card
(You need a P133 or better to do this), or will produce a level
compressed copy of the file to stdout.  Sorry that it does not
interpret the wav header or anything like that -- it does what I
want/need for now!!!

I think that you'll notice that this compressor is very very subtile
in its effect.  Also, the control signal is very carefully filtered
to almost totally eliminate audio frequency ripple (holds distortion
at bay.)  This uses true-RMS detection, with fairly long-term averaging.
Delays were added to eliminate most of the negative effects of a
long attack time.

Along with a fairly sophisticated compressor algorithm, this code
has a peak limiter with a few simple tricks to eliminate modulation of
the control signal with audio-freqs also.

This code is work-in-progress, so caveat emptor!!!  Also, the source code
comments itself, so there!!! :-).  I can do things other than VM code :-).

If you want to copy any of your CDs to a cassette tape for the car --
this code will help a little bit...  You might be kind-of suprised at
how this code doesn't totally destroy the sound of the music (at
least by my perception.)  Let me know how this works for you!!!

Have fun!

John
dyson@freebsd.org


/*
 * Copyright (c) 1996, John S. Dyson
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This code (easily) runs realtime on a P5-166 w/EDO, Triton-II on FreeBSD.
 *
 * More info/comments: dyson@freebsd.org
 *
 * This program provides compression of a stereo 16bit audio stream,
 * such as that contained by a 16Bit wav file.  Extreme measures have
 * been taken to make the compression as subtile as possible.  One
 * possible purpose for this code would be to master cassette tapes from
 * CD's for playback in automobiles where dynamic range needs to be
 * restricted.
 *
 * Suitably recoded for an embedded DSP, this would make a killer audio
 * compressor for broadcast or recording.  When writing this code, I
 * ignored the issues of roundoff error or trucation -- Pentiums have
 * really nice FP processors :-).
 */

/*
 * Usage: audcomp [switches] < infile [>outfile]
 *
 * switches:
 *	-p<fpnum>:	percent target is of peak limit.
 *			Essentially specifies how hard the peak limiter
 *			is pushed.  The default of 69% is good.
 *	-t<release>:	release time in seconds.
 *			This compressor is very sophisticated, and actually
 *			the release time is complex.  This is one of the
 *			dominant release time controls, but the actual
 *			release time is dependent on alot of factors regarding
 *			the dynamics of the audio in.
 *	-r<compression-ratio>	This is the compression ratio for the fast
 *			compressor.  Note that I am lazy, and so this is not
 *			really the compression ratio.  -r1.0 is infinity to
 *			one, while the default -r0.50 is 2:1.  Another
 *			really good value is special cased in the code:
 *			-r0.25 is somewhat less than 2:1, and sounds super
 *			smooth.
 *	-R<compression-ratio>	This is the compression ratio for the entire
 *			compressor chain.  The default is -R1.0, and holds
 *			the volume very constant without many nasty
 *			side-effects.  However the dynamics in music are
 *			severely restricted, and a value of -R0.5 might
 *			keep the music more intact.
 *	-g<gain>	This scales input gain.  At low compression ratios,
 *			sources with small signals will not be brought up
 *			to full-scale level.  This param allows you to scale
 *			the input.
 *
 * Many other parameters can be tuned, UTSL!!!
 */

#include <unistd.h>
#include <math.h>
#include <fcntl.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/stat.h>
#include <machine/soundcard.h>


int dspfd;
int filefd;
int ttyfd;

#define DSP "/dev/dsp"

int inptr;
int incnt;
short inbuf[8192];
struct samplebuffer {
	short	left, right;
};


__inline int
getsample(short *right, short *left) {
	if (incnt == 0) {
		incnt = read(0, inbuf, sizeof inbuf);
		if (incnt == 0 || incnt < 0)
			return 0;
		incnt /= sizeof (short);
		inptr = 0;
	}

	*left = inbuf[inptr++];
	*right = inbuf[inptr++];
	incnt -= 2;
	return 1;
}

#define OBSIZE 4096
int outptr;
short outbuf[OBSIZE * sizeof(short)];

__inline void
putsample(short right, short left) {
	outbuf[outptr++] = left;
	outbuf[outptr++] = right;

	if (outptr >= OBSIZE) {
		write(dspfd, outbuf, outptr * sizeof(short));
		outptr = 0;
	}
}

/*
 * Normally, elements with number suffices are progressively
 * slower.
 */
#define NFILT 12
#define NEFILT 17
struct compstruct {
/* Simple level running average */
	double rlevelsq0, rlevelsq1;
	double rlevelsq0filter, rlevelsq1filter;
	double rlevelsqn[NFILT];
	double rlevelsqefilter;
	double rlevelsqe[NEFILT];
	double rlevelsq0ffilter;
	int ndelay; /* delay for rlevelsq0ffilter delay */
	int ndelayptr; /* ptr for the input */
	double *rightdelay;
	double *leftdelay;
/* Simple gain running average */
	double rgain;
	double rgainfilter;
	double lastrgain;
/* Max fast agc gain, slow agc gain */
	double maxfastgain, maxslowgain;
/* Fast gain compression ratio */
/*	Note that .5 is 2:1, 1.0 is infinity (hard) */
	double fastgaincompressionratio;
	double compressionratio;
/* Max level, target level, floor level */
	double maxlevel, targetlevel, floorlevel;
/* Gainriding gain */
	double rmastergain0filter;
	double rmastergain0;
/* Peak limit gain */
	double rpeakgain0, rpeakgain1, rpeakgainfilter;
	int peaklimitdelay, rpeaklimitdelay;
/* Running total gain */
	double totalgain;
/* Idle gain */
	double npeakgain;
/* Compress enabled */
	int compress;
};

void
compinit(struct compstruct *cs) {
	int i;
	cs->rgain = cs->rmastergain0 = 1.0;
	cs->rlevelsq0 = cs->rlevelsq1 = 0;
	cs->compress = 1;
	cs->ndelay = 1.0/cs->rlevelsq0ffilter;
	cs->rightdelay = malloc(sizeof(double) * cs->ndelay);
	cs->leftdelay = malloc(sizeof(double) * cs->ndelay);
	cs->rpeakgain0 = 1.0;
	cs->rpeakgain1 = 1.0;
	cs->rpeaklimitdelay = 0;
	cs->ndelayptr = 0;
	for(i=0;i<NFILT;i++)
		cs->rlevelsqn[i] = 0;
	for(i=0;i<NEFILT;i++)
		cs->rlevelsqe[i] = 0;
}

__inline double
hardlimit(double value, double knee, double limit)
{
	double lrange = (limit - knee);
	double ab = fabs(value);
/*
	if (ab > knee) {
		double abslimit = (limit * 1.1);
		if (ab < abslimit) 
			value = knee + lrange * sin( ((value - knee)/abslimit) * (3.14 / (4*1.1)));
	}
*/
	if (ab >= limit)
		value = value > 0 ? limit : -limit;
	return value;
}

void
compressor(double *righta, double *lefta, struct compstruct *cs) {
	
	double level, levelsq0, levelsq1, levelsqe;
	double gain, qgain, tgain;
	double newright, newleft;
	double efilt;
	double fastgain, slowgain, tslowgain;
	double leveldelta;
	double right, left, rightd, leftd;
	int delayed;
	double nrgain, nlgain, ngain, ngsq;
	double sqrtrpeakgain;
	int i;
	int skipmode;

	right = *righta;
	left = *lefta;

	cs->rightdelay[cs->ndelayptr] = right;
	cs->leftdelay[cs->ndelayptr] = left;
	cs->ndelayptr++;
	if (cs->ndelayptr >= cs->ndelay)
		cs->ndelayptr = 0;
/* enable/disable compression */

	skipmode = 0;
	if (cs->compress == 0) {
		skipmode = 1;
		goto skipagc;
	}
	levelsq0 = (right) * (right) + (left) * (left);
	
	if (levelsq0 > cs->rlevelsq0) {
		cs->rlevelsq0 = (levelsq0 * cs->rlevelsq0ffilter) +
			cs->rlevelsq0 * (1 - cs->rlevelsq0ffilter);
	} else {
		cs->rlevelsq0 = (levelsq0 * cs->rlevelsq0filter) +
			cs->rlevelsq0 * (1 - cs->rlevelsq0filter);
	}

	if (cs->rlevelsq0 <= cs->floorlevel * cs->floorlevel)
		goto skipagc;

	if (cs->rlevelsq0 > cs->rlevelsq1) {
		cs->rlevelsq1 = cs->rlevelsq0;
	} else {
		cs->rlevelsq1 = cs->rlevelsq0 * cs->rlevelsq1filter +
			cs->rlevelsq1 * (1 - cs->rlevelsq1filter);
	}

	cs->rlevelsqn[0] = cs->rlevelsq1;
	for(i=0;i<NFILT-1;i++) {
		if (cs->rlevelsqn[i] > cs->rlevelsqn[i+1])
			cs->rlevelsqn[i+1] = cs->rlevelsqn[i];
		else
			cs->rlevelsqn[i+1] = cs->rlevelsqn[i] * cs->rlevelsq1filter +
				cs->rlevelsqn[i+1] * (1 - cs->rlevelsq1filter);
	}
	

	efilt = cs->rlevelsqefilter;
	levelsqe = cs->rlevelsqe[0] = cs->rlevelsqn[NFILT-1];
	for(i=0;i<NEFILT-1;i++) {
		cs->rlevelsqe[i+1] = cs->rlevelsqe[i] * efilt +
			cs->rlevelsqe[i+1] * (1.0 - efilt);
		if (cs->rlevelsqe[i+1] > levelsqe)
			levelsqe = cs->rlevelsqe[i+1];
		efilt *= 1.0/1.5;
	}

	gain = cs->targetlevel / sqrt(levelsqe);
	if (cs->compressionratio < 0.99) {
		if (cs->compressionratio == 0.50)
			gain = sqrt(gain);
		else
			gain = exp(log(gain) * cs->compressionratio);
	}

	if (gain < cs->rgain)
		cs->rgain = gain * cs->rlevelsqefilter/2 +
				cs->rgain * (1 - cs->rlevelsqefilter/2);
	else
		cs->rgain = gain * cs->rgainfilter +
				cs->rgain * (1 - cs->rgainfilter);

	cs->lastrgain = cs->rgain;
	if ( gain < cs->lastrgain)
		cs->lastrgain = gain;

skipagc:;

	tgain = cs->lastrgain;

	leftd = cs->leftdelay[cs->ndelayptr];
	rightd = cs->rightdelay[cs->ndelayptr];

	fastgain = tgain;
	if (fastgain > cs->maxfastgain)
		fastgain = cs->maxfastgain;

	if (fastgain < 0.0001)
		fastgain = 0.0001;

	if (cs->fastgaincompressionratio == 0.25) {
		qgain = sqrt(sqrt(fastgain));
	} else if (cs->fastgaincompressionratio == 0.5) {
		qgain = sqrt(fastgain);
	} else if (cs->fastgaincompressionratio == 1.0) {
		qgain = fastgain;
	} else {
		qgain = exp(log(fastgain) * cs->fastgaincompressionratio);
	}

	tslowgain = tgain / qgain;
	if (tslowgain > cs->maxslowgain)
		tslowgain = cs->maxslowgain;
	if (tslowgain < cs->rmastergain0)
		cs->rmastergain0 = tslowgain;
	else
		cs->rmastergain0 = tslowgain * cs->rmastergain0filter +
				(1 - cs->rmastergain0filter) * cs->rmastergain0;

	slowgain = cs->rmastergain0;
	if (skipmode == 0)
		cs->npeakgain = slowgain * qgain;

/**/
	newright = rightd * cs->npeakgain;
	if (fabs(newright) >= cs->maxlevel)
		nrgain = cs->maxlevel / fabs(newright);
	else
		nrgain = 1.0;

	newleft = leftd * cs->npeakgain;
	if (fabs(newleft) >= cs->maxlevel)
		nlgain = cs->maxlevel / fabs(newleft);
	else
		nlgain = 1.0;

	ngain = nrgain;
	if (nlgain < ngain)
		ngain = nlgain;

	ngsq = ngain * ngain;
	if (ngsq <= cs->rpeakgain0) {
#if 0
		if (ngsq < cs->rpeakgain0) {
			fprintf(stderr,"*");
		}
#endif
		cs->rpeakgain0 = ngsq /* * 0.50 + cs->rpeakgain0 * 0.50 */;
		cs->rpeaklimitdelay = cs->peaklimitdelay;
	} else if (cs->rpeaklimitdelay == 0) {
		double tnrgain;
		if (nrgain > 1.0)
			tnrgain = 1.0;
		else
			tnrgain = nrgain;
		cs->rpeakgain0 = tnrgain * cs->rpeakgainfilter +
			(1.0 - cs->rpeakgainfilter) * cs->rpeakgain0;
	}

	if (cs->rpeakgain0 <= cs->rpeakgain1) {
		cs->rpeakgain1 = cs->rpeakgain0;
		cs->rpeaklimitdelay = cs->peaklimitdelay;
	} else if (cs->rpeaklimitdelay == 0) {
		cs->rpeakgain1 = cs->rpeakgainfilter * cs->rpeakgain0 +
			(1.0 - cs->rpeakgainfilter) * cs->rpeakgain1;
	} else {
		--cs->rpeaklimitdelay;
	}

	sqrtrpeakgain = sqrt(cs->rpeakgain1);
	cs->totalgain = cs->npeakgain * sqrtrpeakgain;

	right = newright * sqrtrpeakgain;
	*righta = hardlimit(right, 32200, 32767);

	left = newleft * sqrtrpeakgain;
	*lefta = hardlimit(left, 32200, 32767);

	if (right != *righta || left != *lefta) {
		fprintf(stderr,"!");
	}
}

main(int argc, char *argv[]) {
	int i;
	int writecount;
	int toread;
	int parm;
	char *readb;
	short dright, dleft;
	double fpdright, fpdleft;
	double maxgain, mingain, maxlevel;
	int sample;
	int compress = 1;
	char ttyc;
	struct compstruct cs;
	int ttyflag;
	int ch;
	double fratio;
	int fratioflag;
	double ratio;
	int  ratioflag;
	double releasetime;
	int releasetimeflag;
	double peakpercent;
	int peakpercentflag;
	double gain = 1.0;

	peakpercentflag = 0;
	releasetimeflag = 0;
	fratioflag = 0;
	ratioflag = 0;

	while ((ch = getopt(argc, argv, "p:t:g:R:r:")) !=( -1)) {
		switch (ch) {
		case 'p':
			peakpercent = atof(optarg);
			peakpercentflag = 1;
			break;
		case 't':
			releasetime = atof(optarg);
			releasetimeflag = 1;
			break;
			
		case 'r':
			fratio = atof(optarg);
			fratioflag = 1;
			break;
		case 'R':
			ratio = atof(optarg);
			ratioflag = 1;
			break;
		case 'g':
			gain = atof(optarg);
			break;

		case '?':
		default:
			fprintf(stderr,"argument error -- %s\n", ch);
			exit(1);
		}
	}
	argc -= optind;
	argv += optind;
			
	ttyflag = isatty(1);

#ifdef DISPLAY
	if (ttyflag) {
		ttyfd = open("/dev/tty", O_RDWR);
		if (ttyfd < 0) {
			perror("ttyopen");
			exit(1);
		};

		fcntl(ttyfd, F_SETFL, O_NONBLOCK);
	}
#endif

	if (ttyflag) {
		dspfd = open(DSP, O_RDWR);
		if (dspfd < 0) {
			perror("dspopen");
			exit(1);
		}

		parm = 16;
		if (ioctl(dspfd, SOUND_PCM_READ_BITS, &parm)) {
			perror("ioctlrd16");
			exit(1);
		}

		parm = 16;
		if (ioctl(dspfd, SOUND_PCM_WRITE_BITS, &parm)) {
			perror("ioctlwr16");
			exit(1);
		}

		parm = 2;
		if (ioctl(dspfd, SOUND_PCM_READ_CHANNELS, &parm)) {
			perror("ioctlrd2");
			exit(1);
		}

		parm = 2;
		if (ioctl(dspfd, SOUND_PCM_WRITE_CHANNELS, &parm)) {
			perror("ioctlwr2");
			exit(1);
		}

		parm = 44100;
		if (ioctl(dspfd, SOUND_PCM_READ_RATE, &parm)) {
			perror("ioctlsp44100");
			exit(1);
		}

		parm = 44100;
		if (ioctl(dspfd, SOUND_PCM_WRITE_RATE, &parm)) {
			perror("ioctlsp44100");
			exit(1);
		}
	} else {
		dspfd = 1;
	}

	mingain = 10000;
	maxgain = 0;

	/* These filters should filter at least the lowest audio freq */
	cs.rlevelsq0filter = .002;
	cs.rlevelsq1filter = .010;
	/* These are the attack time for the rms measurement */
	cs.rlevelsq0ffilter = .002;

	cs.rlevelsqefilter = .002;

	/*
	 * Linear gain filters as opposed to the level measurement filters
	 * above
	 */
	if (releasetimeflag) {
		if (releasetime < 0.05)
			releasetime = 0.05;
		cs.rgainfilter = 1.0/(releasetime * 44100);
	} else {
		cs.rgainfilter = 0.002;
	}

	/*
	 * compression ratio for fast gain.  This will determine how
	 * much the audio is made more dense.  .5 is equiv to 2:1
	 * compression.  1.0 is equiv to inf:1 compression.
	 */
	if (fratioflag)
		cs.fastgaincompressionratio = fratio;
	else
		cs.fastgaincompressionratio = 0.50;

	if (ratioflag)
		cs.compressionratio = ratio;
	else
		cs.compressionratio = 1.0;

	/*
	 * Limiter level
	 */
	cs.maxlevel = 32000;
	/*
	 * maximum gain for fast compressor
	 */
	cs.maxfastgain = 3;
	/*
	 * maximum gain for slow compressor
	 */
	cs.maxslowgain = 9;
	/*
	 * target level for compression
	 */
	if (peakpercentflag)
		cs.targetlevel = cs.maxlevel * peakpercent / 100.0;
	else
		cs.targetlevel = cs.maxlevel * 0.69;

	/*
	 * Level below which gain tracking shuts off
	 */
	cs.floorlevel = 2000;

	/*
	 * Slow compressor time constants
	 */
	cs.rmastergain0filter = .000003;

	cs.rpeakgainfilter = .001;
	cs.rpeaklimitdelay = 2500;

	compinit(&cs);
	while (getsample(&dright, &dleft)) {

		fpdright = dright;
		fpdright *= gain;

		fpdleft = dleft;
		fpdleft *= gain;

		if (compress)
			compressor(&fpdright, &fpdleft, &cs);

		dright = fpdright;
		dleft = fpdleft;

		putsample(dright, dleft);
		sample++;
		if (cs.totalgain > maxgain)
			maxgain = cs.totalgain;
		if (cs.totalgain < mingain)
			mingain = cs.totalgain;
		if (fpdright > maxlevel)
			maxlevel = fpdright;
		if (fpdleft > maxlevel)
			maxlevel = fpdleft;

#ifdef DISPLAY
		if ((sample % (4410*5)) == 0) {
			double dbmin, dbmax, dbmaster, pctlvl;
			if (mingain < .001)
				mingain = .001;
			if (maxgain < .001)
				maxgain = .001;
			dbmin = 20*log10(mingain);
			dbmax = 20*log10(maxgain);
			dbmaster = 20*log10(cs.rmastergain0);
			printf("min: %5.2f, max: %5.2f, master: %5.2f, max: %6.0f\n",
				dbmin, dbmax, dbmaster, (maxlevel/cs.maxlevel) * 100.0);
			maxgain = 0.0;
			mingain = 10000.0;
			maxlevel = 0.0;
			if (read(ttyfd, &ttyc, 1)) {
				if (ttyc == 'y')
					cs.compress = 1;
				else if (ttyc == 'n')
					cs.compress = 0;
				else if (ttyc == '+') {
					int i;
					for(i=0;i<44100*60;i++)
						getsample(&dright, &dleft);
				}
			}
		}
#endif
	}
	exit(0);
}



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?199611140727.CAA00238>