This is convolvesf.c in view mode; [Download] [Up]
#include <math.h>
#include <stdio.h>
#include <carl/sndio.h>
#include <carl/carl.h>
#include <carl/procom.h>
#define NCMAX 262144
extern int sferror;
extern CSNDFILE *rootsfd;
/*-------------------------------------------------------
convolvesf.c
This program performs fast convolution via the FFT. In
contrast to fastflt, the filter impulse response
is contained in a soundfile as opposed to a filter
file. One potential problem is getting sufficient
memory for the program to run. NCMAX = 262144
gives a maximum impulse response of 5.3 seconds
at 48K, but this requires over 5 Mbytes of memory
and probably won't run! NOTE: This program depends
upon the C versions of the IEEE FFT subroutines
fast and fsst in order to perform FFT's of greater
than 32K.
-------------------------------------------------------*/
main(argc, argv)
int argc;
char **argv;
{
float fsndi();
float *sbuf, /* array for input and FFT */
*tbuf, /* array for overlap-adding output */
*filt; /* array for filter transform */
float a, /* temporary */
b, /* temporary */
temp, /* temporary */
srate, /* sample rate from header */
gain = 1., /* gain */
max = 0.; /* maximum */
long i,
j,
icnt, /* counts input samples */
ocnt = 0, /* counts output samples */
begin = 0, /* first sample of soundfile */
end = 0, /* last sample of soundfile */
N, /* FFT size is twice impulse response */
N2, /* N / 2 */
ip, /* temporary pointer */
ip1; /* temporary pointer */
int nchan; /* number of channels from header */
CSNDFILE *sfd, *opensf();
char ch;
char *cbeg = NULL, *cend = NULL, *cdur = NULL;
char *file = "test";
char cbuf[72]; /* buffer for strings to write to header */
char *dbuf; /* buffer for strings to read from header */
PROP *proplist; /* from header on stdin */
/* Interpret commandline by calling dgl's subroutine crack. */
if (isatty(0)) usage(1);
while ((ch = crack(argc, argv, "b|d|e|g|h", 0)) != NULL) {
switch (ch) {
case 'b': cbeg = arg_option; break;
case 'e': cend = arg_option; break;
case 'd': cdur = arg_option; break;
case 'g': gain = sfexpr(arg_option,1.); break;
case 'h': usage(0);
default: usage(1); /* this exits with error */
}
}
/* Read header from stdin. */
if ((proplist = getheader(stdin)) != NULL) { /* there is a header */
noautocp(); /* suppress hdr copy */
if ((dbuf=getprop(stdin, H_SRATE)) != NULL) srate = atof(dbuf);
if ((dbuf=getprop(stdin, H_NCHANS)) != NULL)nchan = atoi(dbuf);
}
/* Read in soundfile with room response. */
if (arg_index < argc)
file = argv[arg_index];
else fprintf(stderr,"convolvesf: no impulse response\n");
/* open sound file */
if ((sfd = opensf(file, "-r")) == NULL) {
fprintf(stderr,"convolvesf: opensf failed on %s\n", file);
(void) sfallclose();
exit(1);
}
/* calculate times */
if (cbeg != NULL)
begin = sfexpr(cbeg, sfd->sr)*sfd->nc;
if (cend != NULL)
end = sfexpr(cend, sfd->sr)*sfd->nc;
else if (cdur != NULL)
end = sfexpr(cdur, sfd->sr)*sfd->nc + begin;
/* check boundaries */
if (begin < 0 || begin >= sfd->fs) {
fprintf(stderr, "convolvesf: begin time out of range\n");
quit();
}
if (end < 0 || end >= sfd->fs) {
fprintf(stderr, "convolvesf: end time out of range\n");
quit();
}
if (end == 0)
end = sfd->fs;
if ((end - begin) > NCMAX){
end = begin + NCMAX;
fprintf(stderr,"convolvesf: warning - impulse response too long\n");
}
/* Write header to stdout. */
if (srate != sfd->sr)
fprintf(stderr,"convolvesf: warning - incompatible sample rates\n");
if ((nchan > 1) || (sfd->nc > 1))
fprintf(stderr,"convolvesf: warning - mono expected and not found\n");
cpoheader(proplist,stdout);
putheader(stdout);
/* Set up buffers. */
for (N2 = 2; N2 < NCMAX; N2 *= 2)
if (N2 >= (end - begin)) break;
N = 2 * N2;
if ((sbuf = (float *) calloc(N+2, sizeof(float))) == NULL)
malerr("convolvesf: insufficient memory", 1);
if ((tbuf = (float *) calloc(N2, sizeof(float))) == NULL)
malerr("convolvesf: insufficient memory", 1);
if ((filt = (float *) calloc(N+2, sizeof(float))) == NULL)
malerr("convolvesf: insufficient memory", 1);
/* fill filter buffer */
for (i = begin, j = 0; i < end; i++, j++) {
temp = fsndi(sfd, i);
if (sferror) quit();
*(filt+j) = temp;
}
/* close all open files */
(void) sfallclose();
/* Finish initialization: fill buffers, take FFT of filt, normalize */
for (i=j; i<N+2; i++)
*(filt+i) = 0.;
cfast(filt,N);
for (i=0; i<=N2; i++){
a = *(filt + 2*i);
b = *(filt + 2*i + 1);
temp = a*a + b*b;
if (temp > max)
max = temp;
}
if (max != 0.)
max = gain/(sqrt(max));
else {
fprintf(stderr,"convolvesf: impulse response is all zeros\n");
exit(1);
}
for (i=0; i< N+2; i++)
*(filt + i) *= max;
for (i=0; i<N2; i++)
if (getfloat(sbuf+i) <= 0)
break;
icnt = i;
for (; i<N+2; i++)
*(sbuf+i) = 0.;
if (icnt == 0){
fprintf(stderr,"convolvesf: no valid input on stdin\n");
exit(1);
}
/* Main loop: Take N-point FFT of next N/2 input samples and multiply
by FFT of filter impulse response. Inverse FFT and add first N/2
resulting samples to last N/2 samples of previous FFT. */
while(ocnt < icnt){
cfast(sbuf,N);
for (i=0; i<=N2; i++){
ip = 2*i;
ip1 = ip + 1;
a = *(sbuf+ip)* *(filt+ip) - *(sbuf+ip1)* *(filt+ip1);
b = *(sbuf+ip)* *(filt+ip1) + *(sbuf+ip1)* *(filt+ip);
*(sbuf+ip) = a;
*(sbuf+ip1) = b;
}
cfsst(sbuf,N);
for (i=0; i<N2; i++)
*(tbuf+i) += *(sbuf+i);
for (i=0; i<N2; i++)
if (++ocnt <= icnt)
putfloat(tbuf+i);
for (i=0; i<N2; i++)
*(tbuf+i) = *(sbuf+N2+i);
for (i=0; i<N2; i++)
if (getfloat(sbuf+i) <= 0)
break;
icnt += i;
for (; i<N+2; i++)
*(sbuf+i) = 0.;
}
flushfloat();
exit(0);
}
usage(exitcode)
int exitcode;
{
fprintf(stderr,"%s%s%s%s%s%s",
"\nusage: convolvesf [flags] impulse_response_soundfile < floatsams > floatsams\n",
"\nflags:\n",
"\tg = gain factor (1.)\n",
"\tb = begin time in impulse response (first sample to use) (0.)\n",
"\te = end time in impulse response (last sample to use) (end)\n",
"\td = duration of impulse response (end - begin)\n\n");
exit(exitcode);
}
quit() {
fprintf(stderr, "exiting\n");
(void) sfallclose();
exit(1);
}
malerr(str, ex)
char *str;
int ex;
{
fprintf(stderr, "%s\n", str);
fprintf(stderr,"\nconvolvesf: impulse response too long\n");
exit(ex);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.