1 /*
2 dnoise.c:
3
4 Copyright (C) 2000 Mark Dolson, John ffitch
5
6 This file is part of Csound.
7
8 The Csound Library is free software; you can redistribute it
9 and/or modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2.1 of the License, or (at your option) any later version.
12
13 Csound is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with Csound; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
21 02110-1301 USA
22 */
23
24 /*
25 * PROGRAM: dnoise - de-noise a recording
26 *
27 * AUTHOR: Mark Dolson
28 *
29 * DATE: August 26, 1989
30 *
31 * COMMENTS: dnoise takes floats from stdin and outputs them
32 * on stdout as a noise-reduced version of the input signal.
33 * dnoise uses the phase vocoder algorithm in which
34 * successsive windows are Fast Fourier Transformed,
35 * noise-gated, and then Inverse Fast Fourier Transformed
36 * and overlap-added back together.
37 *
38 * REVISIONS: John ffitch, September 1999, December 2000
39 * Writes any format using usual Csound functions.
40 *
41 */
42
43 /*
44 This is a noise reduction scheme using frequency-
45 domain noise-gating. This should work best in
46 the case of high signal-to-noise with hiss-type
47 noise. The algorithm is that suggested by
48 Moorer & Berger in "Linear-Phase Bandsplitting:
49 Theory and Applications" presented at the 76th
50 Convention 1984 October 8-11 New York of the Audio
51 Engineering Society (preprint #2132) except that
52 it uses the Weighted Overlap-Add formulation for
53 short-time Fourier analysis-synthesis in place of
54 the recursive formulation suggested by Moorer &
55 Berger. The gain in each frequency bin is computed
56 independently according to
57
58 gain = g0 + (1-g0) * [avg / (avg + th*th*nref)] ^ sh
59
60 where avg and nref are the mean squared signal and
61 noise respectively for the bin in question. (This
62 is slightly different than in Moorer & Berger.) The
63 critical parameters th and g0 are specified in dB
64 and internally converted to decimal values. The nref
65 values are computed at the start of the program on the
66 basis of a noise_soundfile (specified in the command
67 line) which contains noise without signal. The avg
68 values are computed over a rectangular window of m
69 FFT frames looking both ahead and behind the current
70 time. This corresponds to a temporal extent of m*D/R
71 (which is typically (m*N/8)/R). The default settings
72 of N, M, and D should be appropriate for most uses. A
73 higher sample rate than 16KHz might indicate a higher N.
74 */
75
76 #include "std_util.h"
77 #include "soundio.h"
78 #include <math.h>
79 #include <ctype.h>
80 #include <inttypes.h>
81
82
83 #define ERR(x) \
84 { \
85 csound->Message(csound, "%s", x); \
86 return -1; \
87 }
88
89 #define FIND(x) \
90 { \
91 if (*s == '\0') { \
92 if (UNLIKELY(!(--argc) || (((s = *argv++) != NULL) && *s == '-'))) { \
93 csound->Message(csound, "%s\n", Str(x)); \
94 return dnoise_usage(csound, -1); \
95 } \
96 } \
97 }
98
99 static int32_t dnoise_usage(CSOUND *, int32_t);
100 static void hamming(MYFLT *, int32_t, int32_t);
101
102 static int32_t writebuffer(CSOUND *, SNDFILE *, MYFLT *,
103 int32_t, int32_t *, OPARMS *);
104
105 #if 0
106 static void fast(CSOUND *csound, MYFLT *b, int32_t N)
107 {
108 /* The DC term is returned in location b[0] with b[1] set to 0.
109 Thereafter, the i'th harmonic is returned as a complex
110 number stored as b[2*i] + j b[2*i+1]. The N/2 harmonic
111 is returned in b[N] with b[N+1] set to 0. Hence, b must
112 be dimensioned to size N+2. The subroutine is called as
113 fast(b,N) where N=2**M and b is the real array described
114 above.
115 */
116
117 csound->RealFFT(csound, b, N);
118 b[N] = b[1];
119 b[1] = b[N + 1] = FL(0.0);
120 }
121
122
123 static void fsst(CSOUND *csound, MYFLT *b, int32_t N)
124 {
125
126 /* This subroutine synthesizes the real vector b[k] for k=0, 1,
127 ..., N-1 from the fourier coefficients stored in the b
128 array of size N+2. The DC term is in location b[0] with
129 b[1] equal to 0. The i'th harmonic is a complex number
130 stored as b[2*i] + j b[2*i+1]. The N/2 harmonic is in
131 b[N] with b[N+1] equal to 0. The subroutine is called as
132 fsst(b,N) where N=2**M and b is the real array described
133 above.
134 */
135 MYFLT scaleVal;
136 int32_t i;
137
138 scaleVal = csound->GetInverseRealFFTScale(csound, N);
139 b[1] = b[N];
140 b[N] = b[N + 1] = FL(0.0);
141 for (i = 0; i < N; i++)
142 b[i] *= scaleVal;
143 csound->InverseRealFFT(csound, b, N);
144 }
145 #endif
146
fast2(CSOUND * csound,void * setup,MYFLT * b)147 static inline void fast2(CSOUND *csound, void *setup, MYFLT *b)
148 {
149 csound->RealFFT2(csound, setup, b);
150 }
151
fsst2(CSOUND * csound,void * setup,MYFLT * b)152 static inline void fsst2(CSOUND *csound, void *setup, MYFLT *b)
153 {
154 csound->RealFFT2(csound, setup, b);
155 }
156
157
dnoise(CSOUND * csound,int32_t argc,char ** argv)158 static int32_t dnoise(CSOUND *csound, int32_t argc, char **argv)
159 {
160 OPARMS O;
161 MYFLT beg = -FL(1.0), end = -FL(1.0);
162 int64_t Beg = 0, End = 99999999;
163
164 MYFLT
165 *ibuf1, /* pointer to start of input buffer */
166 *ibuf2, /* pointer to start of input buffer */
167 *obuf1, /* pointer to start of output buffer */
168 *obuf2, /* pointer to start of output buffer */
169 *fbuf, /* pointer to start of FFT buffer */
170 *aWin, /* pointer to center of analysis window */
171 *sWin, /* pointer to center of synthesis window */
172 *i0, /* pointer to real channels */
173 *i1, /* pointer to imaginary channels */
174 *j0, /* pointer to real channels */
175 *j1, /* pointer to imaginary channels */
176 *f, /* pointer to FFT buffer */
177 *f0, /* pointer to real channels */
178 *f1, /* pointer to imaginary channels */
179 *w, /* pointer to window */
180 *mbuf, /* m most recent frames of FFT */
181 *nbuf, /* m most recent frames of FFT */
182 *nref, /* noise reference buffer */
183 *rsum, /* running sum of magnitude-squared spectrum */
184 *ssum, /* running sum of magnitude-squared spectrum */
185 *ibp, /* pointer to next input to be read */
186 *ib0, /* pointer to output buffer */
187 *ib1, /* pointer to output buffer */
188 *ib2, /* pointer to output buffer */
189 *obp, /* pointer to next output to be read */
190 *ob0, /* pointer to output buffer */
191 *ob1, /* pointer to output buffer */
192 *ob2; /* pointer to output buffer */
193
194 int32_t
195 N = 0, /* number of phase vocoder channels (bands) */
196 Np2, /* N+2 */
197 M = 0, /* length of aWin impulse response */
198 L = 0, /* length of sWin impulse response */
199 D = 0, /* decimation factor (default will be M/8) */
200 I = 0, /* interpolation factor (default will be I=D)*/
201 W = -1, /* filter overlap factor (determines M, L) */
202 ibuflen, /* size of ibuf */
203 obuflen, /* size of obuf */
204 aLen, /* half-length of analysis window */
205 sLen; /* half-length of synthesis window */
206
207 int64_t
208 oCnt = 0L, /* number of samples written to output */
209 nI, /* current input (analysis) sample */
210 nO, /* current output (synthesis) sample */
211 nImodR, /* current input sample mod R */
212 nMaxOut, /* last output (synthesis) sample */
213 nMin, /* first input (analysis) sample */
214 nMax, /* last input sample (unless EOF) */
215 lnread, /* total input samples read */
216 lj, /* to satisfy lame Microsoft compiler */
217 lk; /* to satisfy lame Microsoft compiler */
218
219 SNDFILE *fp = NULL; /* noise reference file */
220
221 MYFLT
222 Ninv, /* 1. / N */
223 sum, /* scale factor for renormalizing windows */
224 //rIn, /* decimated sampling rate */
225 //rOut, /* pre-interpolated sampling rate */
226 invR, /* 1. / srate */
227 time, /* nI / srate */
228 gain, /* gain of noise gate */
229 g0 = -FL(40.0),/* minimum gain for noise gate */
230 g0m, /* 1. - g0 */
231 th = FL(30.0), /* threshold above noise reference (dB) */
232 avg, /* average square energy */
233 fac, /* factor in gain computation */
234 minv, /* 1 / m */
235 R = -FL(1.0); /* input sampling rate */
236
237 int32_t i,j,k, /* index variables */
238 ibs, /* current starting location in input buffer */
239 ibc, /* current location in input buffer */
240 obs, /* current starting location in output buffer */
241 obc, /* current location in output buffer */
242 m = 5, /* number of frames to save in mbuf */
243 mi = 0, /* frame offset index in mbuf */
244 mj, /* delayed offset index in mbuf */
245 md, /* number of frames of delay in mbuf (m/2) */
246 mp, /* mi * Np2 */
247 sh = 1, /* sharpness control for noise gate gain */
248 nread, /* number of bytes read */
249 N2, /* N/2 */
250 Meven = 0, /* flag for even M */
251 Leven = 0, /* flag for even L */
252 Verbose = 0,/* flag for verbose output to stderr */
253 Chans = -1, /* number of audio input channels (stereo = 2) */
254 chan, /* channel counter */
255 flag = 1, /* end-of-input flag */
256 first = 0; /* first-time-thru flag */
257
258 SOUNDIN *p, *pn;
259 char *infile = NULL, *nfile = NULL;
260 SNDFILE *inf = NULL, *outfd = NULL;
261 char c, *s;
262 int32_t channel = ALLCHNLS;
263 MYFLT beg_time = FL(0.0), input_dur = FL(0.0), sr = FL(0.0);
264 MYFLT beg_ntime = FL(0.0), input_ndur = FL(0.0), srn = FL(0.0);
265 const char *envoutyp = NULL;
266 uint32_t outbufsiz = 0U;
267 int32_t nrecs = 0;
268 csound->GetOParms(csound, &O);
269
270
271 /* audio is now normalised after call to getsndin */
272 /* csound->e0dbfs = csound->dbfs_to_float = FL(1.0); */
273
274 if ((envoutyp = csound->GetEnv(csound, "SFOUTYP")) != NULL) {
275 if (strcmp(envoutyp, "AIFF") == 0)
276 O.filetyp = TYP_AIFF;
277 else if (strcmp(envoutyp, "WAV") == 0)
278 O.filetyp = TYP_WAV;
279 else if (strcmp(envoutyp, "IRCAM") == 0)
280 O.filetyp = TYP_IRCAM;
281 else {
282 csound->Message(csound, Str("%s not a recognised SFOUTYP env setting"),
283 envoutyp);
284 return -1;
285 }
286 }
287 {
288 ++argv;
289 while (--argc>0) {
290 s = *argv++;
291 if (*s++ == '-') { /* read all flags: */
292 while ((c = *s++) != '\0') {
293 switch (c) {
294 case 'o':
295 FIND("no outfilename");
296 O.outfilename = s; /* soundout name */
297 for ( ; *s != '\0'; s++) ;
298 if (UNLIKELY(strcmp(O.outfilename, "stdin") == 0)) {
299 csound->Message(csound, "%s", Str("-o cannot be stdin\n"));
300 return -1;
301 }
302 break;
303 case 'i':
304 FIND("no noisefilename");
305 nfile = s;
306 for ( ; *s != '\0'; s++) ;
307 break;
308 case 'A':
309 if (UNLIKELY(O.filetyp == TYP_WAV))
310 csound->Warning(csound,
311 "%s", Str("-A overriding local default WAV out"));
312 O.filetyp = TYP_AIFF; /* AIFF output request*/
313 break;
314 case 'J':
315 if (UNLIKELY(O.filetyp == TYP_AIFF || O.filetyp == TYP_WAV))
316 csound->Warning(csound, "%s", Str("-J overriding local default "
317 "AIFF/WAV out"));
318 O.filetyp = TYP_IRCAM; /* IRCAM output request */
319 break;
320 case 'W':
321 if (UNLIKELY(O.filetyp == TYP_AIFF))
322 csound->Warning(csound,
323 "%s", Str("-W overriding local default AIFF out"));
324 O.filetyp = TYP_WAV; /* WAV output request */
325 break;
326 case 'h':
327 O.filetyp = TYP_RAW;
328 O.sfheader = 0; /* skip sfheader */
329 break;
330 case 'c':
331 O.outformat = AE_CHAR; /* 8-bit char soundfile */
332 break;
333 case '8':
334 O.outformat = AE_UNCH; /* 8-bit unsigned char file */
335 break;
336 case 'a':
337 O.outformat = AE_ALAW; /* a-law soundfile */
338 break;
339 case 'u':
340 O.outformat = AE_ULAW; /* mu-law soundfile */
341 break;
342 case 's':
343 O.outformat = AE_SHORT; /* short_int soundfile */
344 break;
345 case 'l':
346 O.outformat = AE_LONG; /* long_int soundfile */
347 break;
348 case 'f':
349 O.outformat = AE_FLOAT; /* float soundfile */
350 break;
351 case 'R':
352 O.rewrt_hdr = 1;
353 break;
354 case 'H':
355 if (isdigit(*s)) {
356 int32_t n;
357 sscanf(s, "%d%n", &O.heartbeat, &n);
358 s += n;
359 }
360 else O.heartbeat = 1;
361 break;
362 case 't':
363 FIND(Str("no t argument"));
364 #if defined(USE_DOUBLE)
365 csound->sscanf(s,"%lf",&th);
366 #else
367 csound->sscanf(s,"%f",&th);
368 #endif
369 while (*++s);
370 break;
371 case 'S':
372 FIND("no s arg");
373 sscanf(s,"%d", &sh);
374 while (*++s);
375 break;
376 case 'm':
377 FIND("no m arg");
378 #if defined(USE_DOUBLE)
379 csound->sscanf(s,"%lf",&g0);
380 #else
381 csound->sscanf(s,"%f",&g0);
382 #endif
383 while (*++s);
384 break;
385 case 'n':
386 FIND(Str("no n argument"));
387 sscanf(s,"%d", &m);
388 while (*++s);
389 break;
390 case 'b':
391 FIND(Str("no b argument"));
392 #if defined(USE_DOUBLE)
393 csound->sscanf(s,"%lf",&beg);
394 #else
395 csound->sscanf(s,"%f",&beg);
396 #endif
397 while (*++s);
398 break;
399 case 'B': FIND(Str("no B argument"));
400 sscanf(s,"%" SCNd64, &Beg);
401 while (*++s);
402 break;
403 case 'e': FIND("no e arg");
404 #if defined(USE_DOUBLE)
405 csound->sscanf(s,"%lf",&end);
406 #else
407 csound->sscanf(s,"%f",&end);
408 #endif
409 while (*++s);
410 break;
411 case 'E': FIND(Str("no E argument"));
412 sscanf(s,"%" PRId64, &End);
413 while (*++s);
414 break;
415 case 'N': FIND(Str("no N argument"));
416 sscanf(s,"%d", &N);
417 while (*++s);
418 break;
419 case 'M': FIND(Str("no M argument"));
420 sscanf(s,"%d", &M);
421 while (*++s);
422 break;
423 case 'L': FIND(Str("no L argument"));
424 sscanf(s,"%d", &L);
425 while (*++s);
426 break;
427 case 'w': FIND(Str("no w argument"));
428 sscanf(s,"%d", &W);
429 while (*++s);
430 break;
431 case 'D': FIND(Str("no D argument"));
432 sscanf(s,"%d", &D);
433 while (*++s);
434 break;
435 case 'V':
436 Verbose = 1; break;
437 default:
438 csound->Message(csound, Str("Looking at %c\n"), c);
439 return dnoise_usage(csound, -1); /* this exits with error */
440 }
441 }
442 }
443 else if (infile==NULL) {
444 infile = --s;
445 csound->Message(csound, Str("Infile set to %s\n"), infile);
446 }
447 else {
448 csound->Message(csound, Str("End with %s\n"), s);
449 return dnoise_usage(csound, -1);
450 }
451 }
452 }
453 if (UNLIKELY(infile == NULL)) {
454 csound->Message(csound, "%s", Str("dnoise: no input file\n"));
455 return dnoise_usage(csound, -1);
456 }
457 if (UNLIKELY(nfile == NULL)) {
458 csound->Message(csound, "%s",
459 Str("Must have an example noise file (-i name)\n"));
460 return -1;
461 }
462 if (UNLIKELY((inf = csound->SAsndgetset(csound, infile, &p, &beg_time,
463 &input_dur, &sr, channel)) == NULL)) {
464 csound->Message(csound, Str("error while opening %s"), infile);
465 return -1;
466 }
467 if (O.outformat == 0) O.outformat = p->format;
468 O.sfsampsize = csound->sfsampsize(FORMAT2SF(O.outformat));
469 if (O.filetyp == TYP_RAW) {
470 O.sfheader = 0;
471 O.rewrt_hdr = 0;
472 }
473 else
474 O.sfheader = 1;
475 if (O.outfilename == NULL)
476 O.outfilename = "test";
477 {
478 SF_INFO sfinfo;
479 char *name;
480 memset(&sfinfo, 0, sizeof(SF_INFO));
481 sfinfo.samplerate = (int32_t) p->sr;
482 sfinfo.channels = (int32_t) p->nchanls;
483 sfinfo.format = TYPE2SF(O.filetyp) | FORMAT2SF(O.outformat);
484 if (strcmp(O.outfilename, "stdout") != 0) {
485 name = csound->FindOutputFile(csound, O.outfilename, "SFDIR");
486 if (name == NULL) {
487 csound->Message(csound, Str("cannot open %s.\n"), O.outfilename);
488 return -1;
489 }
490 outfd = sf_open(name, SFM_WRITE, &sfinfo);
491 if (outfd != NULL)
492 csound->NotifyFileOpened(csound, name,
493 csound->type2csfiletype(O.filetyp, O.outformat), 1, 0);
494 csound->Free(csound, name);
495 }
496 else
497 outfd = sf_open_fd(1, SFM_WRITE, &sfinfo, 1);
498 if (UNLIKELY(outfd == NULL)) {
499 csound->Message(csound, Str("cannot open %s."), O.outfilename);
500 return -1;
501 }
502 /* register file to be closed by csoundReset() */
503 (void)csound->CreateFileHandle(csound, &outfd, CSFILE_SND_W,
504 O.outfilename);
505 sf_command(outfd, SFC_SET_CLIPPING, NULL, SF_TRUE);
506 }
507
508 csound->SetUtilSr(csound, (MYFLT)p->sr);
509 csound->SetUtilNchnls(csound, Chans = p->nchanls);
510
511 /* read header info */
512 if (R < FL(0.0))
513 R = (MYFLT)p->sr;
514 if (Chans < 0)
515 Chans = (int32_t) p->nchanls;
516 p->nchanls = Chans;
517
518 if (UNLIKELY(Chans > 2)) {
519 csound->Message(csound, "%s", Str("dnoise: input MUST be mono or stereo\n"));
520 return -1;
521 }
522
523 /* read noise reference file */
524
525 if (UNLIKELY((fp = csound->SAsndgetset(csound, nfile, &pn, &beg_ntime,
526 &input_ndur, &srn, channel)) == NULL)) {
527 csound->Message(csound, "%s",
528 Str("dnoise: cannot open noise reference file\n"));
529 return -1;
530 }
531
532 if (UNLIKELY(sr != srn)) {
533 csound->Message(csound, "%s", Str("Incompatible sample rates\n"));
534 return -1;
535 }
536 /* calculate begin and end times in NOISE file */
537 if (beg >= FL(0.0)) Beg = (int64_t) (beg * R);
538 if (end >= FL(0.0)) End = (int64_t) (end * R);
539 else if (End == 99999999) End = (int64_t) (input_ndur * R);
540
541 nMin = Beg * Chans; /* total number of samples to skip */
542 nMax = End - Beg; /* number of samples per channel to process */
543
544 /* largest valid FFT size is 8192 */
545 if (N == 0)
546 N = 1024;
547 for (i = 1; i < 4096; i *= 2)
548 if (i >= N)
549 break;
550 if (UNLIKELY(i != N))
551 csound->Message(csound,
552 Str("dnoise: warning - N not a valid power of two; "
553 "revised N = %d\n"),i);
554 //FFT setup
555 //printf("NNN %d \n", N);
556 void *fftsetup_fwd = csound->RealFFT2Setup(csound,N,FFT_FWD);
557 void *fftsetup_inv = csound->RealFFT2Setup(csound,N,FFT_INV);
558
559 N = i;
560 N2 = N / 2;
561 Np2 = N + 2;
562 Ninv = FL(1.0) / N;
563
564 if (W != -1) {
565 if (UNLIKELY(M != 0))
566 csound->Message(csound, "%s",
567 Str("dnoise: warning - do not specify both M and W\n"));
568 else if (W == 0)
569 M = 4*N;
570 else if (W == 1)
571 M = 2*N;
572 else if (W == 2)
573 M = N;
574 else if (W == 3)
575 M = N2;
576 else
577 csound->Message(csound, "%s", Str("dnoise: warning - invalid W ignored\n"));
578 }
579
580 if (M == 0)
581 M = N;
582 if ((M%2) == 0)
583 Meven = 1;
584
585 if (L == 0)
586 L = M;
587 if ((L%2) == 0)
588 Leven = 1;
589
590 if (UNLIKELY(M < 7)) {
591 csound->Message(csound, "%s", Str("dnoise: warning - M is too small\n"));
592 exit(~1);
593 }
594 if (D == 0)
595 D = M / 8;
596
597 I = D;
598
599 lj = (int64_t) M + 3 * (int64_t) D;
600 lj *= (int64_t) Chans;
601 if (UNLIKELY(lj > 32767)) {
602 csound->Message(csound, "%s", Str("dnoise: M too large\n"));
603 return -1;
604 }
605 lj = (int64_t) L + 3 * (int64_t) I;
606 lj *= (int64_t) Chans;
607 if (UNLIKELY(lj > 32767)) {
608 csound->Message(csound, "%s", Str("dnoise: L too large\n"));
609 return -1;
610 }
611
612 ibuflen = Chans * (M + 3 * D);
613 obuflen = Chans * (L + 3 * I);
614 outbufsiz = obuflen * sizeof(MYFLT); /* calc outbuf size */
615 #if 0
616 outbuf = csound->Malloc(csound, (size_t) outbufsiz); /* & alloc bufspace */
617 #endif
618 csound->Message(csound, Str("writing %u-byte blks of %s to %s"),
619 outbufsiz, csound->getstrformat(O.outformat),
620 O.outfilename);
621 csound->Message(csound, " (%s)\n", csound->type2string(O.filetyp));
622 /* spoutran = spoutsf; */
623
624 minv = FL(1.0) / (MYFLT)m;
625 md = m / 2;
626 g0 = (MYFLT) pow(10.0,(double)(0.05*(double)g0));
627 g0m = FL(1.0) - g0;
628 th = (MYFLT) pow(10.0,(double)(0.05*(double)th));
629
630 /* set up analysis window: The window is assumed to be symmetric
631 with M total points. After the initial memory allocation,
632 aWin always points to the midpoint of the window (or one
633 half sample to the right, if M is even); aLen is half the
634 true window length (rounded down). If the window duration
635 is longer than the transform (M > N), then the window is
636 multiplied by a sin(x)/x function to meet the condition:
637 aWin[Ni] = 0 for i != 0. In either case, the
638 window is renormalized so that the phase vocoder amplitude
639 estimates are properly scaled. */
640
641 if (UNLIKELY((aWin =
642 (MYFLT*) csound->Calloc(csound,
643 (M+Meven) * sizeof(MYFLT))) == NULL)) {
644 ERR(Str("dnoise: insufficient memory\n"));
645 }
646
647 aLen = M/2;
648 aWin += aLen;
649
650 hamming(aWin, aLen, Meven);
651 for (i = 1; i <= aLen; i++) {
652 aWin[-i] = aWin[i-1];
653 }
654
655 if (M > N) {
656 if (Meven)
657 *aWin *= (MYFLT)N * (MYFLT) sin(HALFPI/(double)N) /( HALFPI_F);
658 for (i = 1; i <= aLen; i++)
659 aWin[i] *= (MYFLT) (N * sin(PI * ((double) i + 0.5 * (double) Meven)
660 / (double) N)
661 / (PI * (i + 0.5 * (double) Meven)));
662 for (i = 1; i <= aLen; i++)
663 aWin[-i] = aWin[i - Meven];
664 }
665
666 sum = FL(0.0);
667 for (i = -aLen; i <= aLen; i++)
668 sum += aWin[i];
669
670 sum = FL(2.0) / sum; /* factor of 2 comes in later in trig identity */
671
672 for (i = -aLen; i <= aLen; i++)
673 aWin[i] *= sum;
674
675 /* set up synthesis window: For the minimal mean-square-error
676 formulation (valid for N >= M), the synthesis window
677 is identical to the analysis window (except for a
678 scale factor), and both are even in length. If N < M,
679 then an interpolating synthesis window is used. */
680
681 if (UNLIKELY((sWin =
682 (MYFLT*) csound->Calloc(csound,
683 (L+Leven) * sizeof(MYFLT))) == NULL)) {
684 ERR(Str("dnoise: insufficient memory\n"));
685 }
686
687 sLen = L/2;
688 sWin += sLen;
689
690 if (M <= N) {
691 hamming(sWin, sLen, Leven);
692 for (i = 1; i <= sLen; i++)
693 sWin[-i] = sWin[i - Leven];
694
695 for (i = -sLen; i <= sLen; i++)
696 sWin[i] *= sum;
697
698 sum = FL(0.0);
699 for (i = -sLen; i <= sLen; i+=I)
700 sum += sWin[i] * sWin[i];
701
702 sum = FL(1.0) / sum;
703
704 for (i = -sLen; i <= sLen; i++)
705 sWin[i] *= sum;
706 }
707 else {
708 hamming(sWin, sLen, Leven);
709 for (i = 1; i <= sLen; i++)
710 sWin[-i] = sWin[i - Leven];
711
712 if (Leven)
713 *sWin *= (MYFLT) (I * sin(HALFPI/(double) I) / (HALFPI));
714 for (i = 1; i <= sLen; i++)
715 sWin[i] *= (MYFLT)(I * sin(PI * ((double) i + 0.5 * (double) Leven)
716 / (double) I)
717 / (PI * ((double) i + 0.5 * (double) Leven)));
718 for (i = 1; i <= sLen; i++)
719 sWin[i] = sWin[i - Leven];
720
721 sum = FL(1.0) / sum;
722
723 for (i = -sLen; i <= sLen; i++)
724 sWin[i] *= sum;
725 }
726
727 /* set up input buffer: nextIn always points to the next empty
728 word in the input buffer (i.e., the sample following
729 sample number (n + aLen)). If the buffer is full,
730 then nextIn jumps back to the beginning, and the old
731 values are written over. */
732
733 if (UNLIKELY((ibuf1 =
734 (MYFLT *) csound->Calloc(csound,
735 ibuflen * sizeof(MYFLT))) == NULL)) {
736 ERR("dnoise: insufficient memory\n");
737 }
738 if (UNLIKELY((ibuf2 =
739 (MYFLT *) csound->Calloc(csound,
740 ibuflen * sizeof(MYFLT))) == NULL)) {
741 ERR(Str("dnoise: insufficient memory\n"));
742 }
743
744 /* set up output buffer: nextOut always points to the next word
745 to be shifted out. The shift is simulated by writing the
746 value to the standard output and then setting that word
747 of the buffer to zero. When nextOut reaches the end of
748 the buffer, it jumps back to the beginning. */
749
750 if (UNLIKELY((obuf1 =
751 (MYFLT*) csound->Calloc(csound,
752 obuflen * sizeof(MYFLT))) == NULL)) {
753 ERR(Str("dnoise: insufficient memory\n"));
754 }
755 if (UNLIKELY((obuf2 =
756 (MYFLT*) csound->Calloc(csound,
757 obuflen * sizeof(MYFLT))) == NULL)) {
758 ERR(Str("dnoise: insufficient memory\n"));
759 }
760
761 /* set up analysis buffer for (N/2 + 1) channels: The input is real,
762 so the other channels are redundant. */
763
764 if (UNLIKELY((fbuf =
765 (MYFLT*) csound->Calloc(csound, Np2 * sizeof(MYFLT))) == NULL)) {
766 ERR(Str("dnoise: insufficient memory\n"));
767 }
768
769 /* noise reduction: calculate noise reference by taking as many
770 consecutive FFT's as possible in noise soundfile, and
771 averaging them all together. Multiply by th*th to
772 establish threshold for noise-gating in each bin. */
773
774 if (UNLIKELY((nref =
775 (MYFLT*) csound->Calloc(csound,
776 (N2 + 1) * sizeof(MYFLT))) == NULL)) {
777 ERR(Str("dnoise: insufficient memory\n"));
778 }
779
780 if (UNLIKELY((mbuf =
781 (MYFLT*) csound->Calloc(csound,
782 (m * Np2) * sizeof(MYFLT))) == NULL)) {
783 ERR(Str("dnoise: insufficient memory\n"));
784 }
785 if (UNLIKELY((nbuf =
786 (MYFLT*) csound->Calloc(csound,
787 (m * Np2) * sizeof(MYFLT))) == NULL)) {
788 ERR(Str("dnoise: insufficient memory\n"));
789 }
790 if (UNLIKELY((rsum =
791 (MYFLT*) csound->Calloc(csound,
792 (N2 + 1) * sizeof(MYFLT))) == NULL)) {
793 ERR(Str("dnoise: insufficient memory\n"));
794 }
795 if (UNLIKELY((ssum =
796 (MYFLT*) csound->Calloc(csound,
797 (N2 + 1) * sizeof(MYFLT))) == NULL)) {
798 ERR(Str("dnoise: insufficient memory\n"));
799 }
800
801 /* skip over nMin samples */
802 while (nMin > (int64_t)ibuflen) {
803 if (UNLIKELY(!csound->CheckEvents(csound)))
804 csound->LongJmp(csound, 1);
805 nread = csound->getsndin(csound, fp, ibuf1, ibuflen, pn);
806 for(i=0; i < nread; i++)
807 ibuf1[i] *= 1.0/csound->Get0dBFS(csound);
808 if (UNLIKELY(nread < ibuflen)) {
809 ERR(Str("dnoise: begin time is greater than EOF of noise file!"));
810 }
811 nMin -= (int64_t) ibuflen;
812 }
813 if (UNLIKELY(!csound->CheckEvents(csound)))
814 csound->LongJmp(csound, 1);
815 i = (int32_t) nMin;
816 nread = csound->getsndin(csound, fp, ibuf1, i, pn);
817 for(i=0; i < nread; i++)
818 ibuf1[i] *= 1.0/csound->Get0dBFS(csound);
819 if (UNLIKELY(nread < i)) {
820 ERR(Str("dnoise: begin time is greater than EOF of noise file!"));
821 }
822 k = 0;
823 lj = Beg; /* single channel only */
824 while (lj < End) {
825 if (UNLIKELY(!csound->CheckEvents(csound)))
826 csound->LongJmp(csound, 1);
827 lj += (int64_t) N;
828 nread = csound->getsndin(csound, fp, fbuf, N, pn);
829 for(i=0; i < nread; i++)
830 fbuf[i] *= 1.0/csound->Get0dBFS(csound);
831 if (nread < N)
832 break;
833
834 fbuf[N] = FL(0.0);
835 fbuf[N + 1] = FL(0.0);
836
837 //fast(csound, fbuf, N);
838 fast2(csound, fftsetup_fwd, fbuf);
839
840 for (i = 0; i <= N+1; i++)
841 fbuf[i] *= Ninv;
842
843 i0 = fbuf;
844 i1 = i0 + 1;
845 for (i = 0; i <= N2; i++, i0 += 2, i1 += 2) {
846 fac = fbuf[2*i] * fbuf[2*i];
847 fac += fbuf[2*i+1] * fbuf[2*i+1];
848 nref[i] += fac;
849 }
850 k++;
851 }
852 if (UNLIKELY(k == 0)) {
853 ERR(Str("dnoise: not enough samples of noise reference\n"));
854 }
855 fac = th * th / k;
856 for (i = 0; i <= N2; i++)
857 nref[i] *= fac; /* nref[i] *= fac; */
858
859 /* initialization: input time starts negative so that the rightmost
860 edge of the analysis filter just catches the first non-zero
861 input samples; output time equals input time. */
862
863 /* zero ibuf1 to start */
864 memset(ibuf1, '\0', ibuflen*sizeof(MYFLT));
865 /* f = ibuf1; */
866 /* for (i = 0; i < ibuflen; i++, f++) */
867 /* *f = FL(0.0); */
868 if (UNLIKELY(!csound->CheckEvents(csound)))
869 csound->LongJmp(csound, 1);
870 /* fill ibuf2 to start */
871 nread = csound->getsndin(csound, inf, ibuf2, ibuflen, p);
872 /* nread = read(inf, ibuf2, ibuflen*sizeof(MYFLT)); */
873 /* nread /= sizeof(MYFLT); */
874 for(i=0; i < nread; i++)
875 ibuf2[i] *= 1.0/csound->Get0dBFS(csound);
876 lnread = nread;
877 memset(ibuf2+nread, '\0', (ibuflen-nread)*sizeof(MYFLT));
878 /* f = ibuf2 + nread; */
879 /* for (i = nread; i < ibuflen; i++, f++) */
880 /* *f = FL(0.0); */
881
882 //rIn = ((MYFLT) R / D);
883 //rOut = ((MYFLT) R / I);
884 invR = FL(1.0) / R;
885 nI = -((int64_t)aLen / D) * D; /* input time (in samples) */
886 nO = nI; /* output time (in samples) */
887 ibs = ibuflen + Chans * (nI - aLen - 1); /* starting position in ib1 */
888 ib1 = ibuf1; /* filled with zeros to start */
889 ib2 = ibuf2; /* first buffer of speech */
890 obs = Chans * (nO - sLen - 1); /* starting position in ob1 */
891 while (obs < 0) {
892 obs += obuflen;
893 first++;
894 }
895 ob1 = obuf1; /* filled with garbage to start */
896 ob2 = obuf2; /* first output buffer */
897 nImodR = nI; /* for reporting progress */
898 mi = 0;
899 mj = m - md;
900 if (mj >= m)
901 mj = 0;
902 mp = mi * Np2;
903
904 nMax = (int64_t)(input_dur * R); /* Do it all */
905 nMaxOut = (int64_t) (nMax * Chans);
906 while (nI < (nMax + aLen)) {
907
908 time = nI * invR;
909
910 for (chan = 0; chan < Chans; chan++) {
911
912 /* prepare for analysis: always begin reading from ib1 */
913 /* always begin writing to ob1 */
914
915 if (ibs >= ibuflen) { /* done reading from ib1 */
916 if (UNLIKELY(!csound->CheckEvents(csound)))
917 csound->LongJmp(csound, 1);
918 /* swap buffers */
919 ib0 = ib1;
920 ib1 = ib2;
921 ib2 = ib0;
922 ibs -= ibuflen;
923 /* fill ib2 */
924 nread = csound->getsndin(csound, inf, ib2, ibuflen, p);
925 for(i=0; i < nread; i++)
926 ib2[i] *= 1.0/csound->Get0dBFS(csound);
927 lnread += nread;
928 memset(ib2+nread, '\0', (ibuflen-nread)*sizeof(MYFLT));
929 /* f = ib2 + nread; */
930 /* for (i = nread; i < ibuflen; i++, f++) */
931 /* *f = FL(0.0); */
932 }
933 ibc = ibs + chan;
934 ibp = ib1 + ibs + chan;
935
936 if (obs >= obuflen) { /* done writing to ob1 */
937 /* dump ob1 (except at beginning) */
938 if (first > 0) {
939 first--;
940 }
941 else {
942 if ((oCnt + obuflen) < nMaxOut) {
943 oCnt += writebuffer(csound, outfd, ob1, obuflen, &nrecs, &O);
944 }
945 else {
946 i = (int32_t) (nMaxOut - oCnt);
947 oCnt += writebuffer(csound, outfd, ob1, i, &nrecs, &O);
948 }
949 }
950 /* zero ob1 */
951 memset(ob1, '\0', ibuflen*sizeof(MYFLT));
952 /* f = ob1; */
953 /* for (i = 0; i < obuflen; i++, f++) */
954 /* *f = FL(0.0); */
955 /* swap buffers */
956 ob0 = ob1;
957 ob1 = ob2;
958 ob2 = ob0;
959 obs -= obuflen;
960 }
961 obc = obs + chan;
962 obp = ob1 + obs + chan;
963
964 /* analysis: The analysis subroutine computes the complex output at
965 time n of (N/2 + 1) of the phase vocoder channels. It operates
966 on input samples (n - aLen) thru (n + aLen).
967 It expects aWin to point to the center of a
968 symmetric window of length (2 * aLen + 1). It is the
969 responsibility of the main program to ensure that these values
970 are correct. The results are returned in fbuf as succesive
971 pairs of real and imaginary values for the lowest (N/2 + 1)
972 channels. The subroutine fast implements an
973 efficient FFT call for a real input sequence. */
974
975 memset(fbuf, '\0', (N+2)*sizeof(MYFLT));
976 /* f = fbuf; */
977 /* for (i = 0; i < N+2; i++, f++) */
978 /* *f = FL(0.0); */
979
980 lk = nI - (int64_t) aLen - 1; /*time shift*/
981 while ((int64_t) lk < 0L)
982 lk += (int64_t) N;
983 k = (int32_t) (lk % (int64_t) N);
984
985 f = fbuf + k;
986 w = aWin - aLen;
987 for (i = -aLen; i <= aLen; i++, k++, f++, w++) {
988 ibp += Chans;
989 ibc += Chans;
990 if (ibc >= ibuflen) {
991 ibc = chan;
992 ibp = ib2 + chan;
993 }
994 if (k >= N) {
995 k = 0;
996 f = fbuf;
997 }
998 *f += *w * *ibp;
999 }
1000
1001 //fast(csound, fbuf, N);
1002 fast2(csound, fftsetup_fwd, fbuf);
1003
1004 /* noise reduction: for each bin, calculate average magnitude-squared
1005 and calculate corresponding gain. Apply this gain to delayed
1006 FFT values in mbuf[mj*Np2 + i?]. */
1007
1008 if (chan == 0) {
1009 f = rsum;
1010 i0 = mbuf + mp;
1011 i1 = i0 + 1;
1012 j0 = mbuf + mj * Np2;
1013 j1 = j0 + 1;
1014 f0 = fbuf;
1015 f1 = f0 + 1;
1016 for (i = 0; i <= N2;
1017 i++, f++, i0+=2, i1+=2, f0+=2, f1+=2, j0+=2, j1+=2) {
1018 /*
1019 * ii0 = 2 * i; // better as in by 2 or shift?
1020 * ii1 = ii0 + 1;
1021 *
1022 * rsum[i] -= mbuf[mp + ii0] * mbuf[mp + ii0];
1023 * rsum[i] -= mbuf[mp + ii1] * mbuf[mp + ii1];
1024 * rsum[i] += fbuf[ii0] * fbuf[ii0];
1025 * rsum[i] += fbuf[ii1] * fbuf[ii1];
1026 */
1027 *f -= *i0 * *i0;
1028 *f -= *i1 * *i1;
1029 *f += *f0 * *f0;
1030 *f += *f1 * *f1;
1031 avg = minv * *f; /* avg = minv * rsum[i]; */
1032 if (avg < FL(0.0))
1033 avg = FL(0.0);
1034 if (avg == FL(0.0))
1035 fac = FL(0.0);
1036 else
1037 fac = avg / (avg + nref[i]);
1038 for (j = 1; j < sh; j++)
1039 fac *= fac;
1040 gain = g0m * fac + g0;
1041 /*
1042 * mbuf[mp + ii0] = fbuf[ii0];
1043 * mbuf[mp + ii1] = fbuf[ii1];
1044 * fbuf[ii0] = gain * mbuf[mj*Np2 + ii0];
1045 * fbuf[ii1] = gain * mbuf[mj*Np2 + ii1];
1046 */
1047 *i0 = *f0;
1048 *i1 = *f1;
1049 *f0 = gain * *j0;
1050 *f1 = gain * *j1;
1051 }
1052 }
1053 else {
1054 f = ssum;
1055 i0 = nbuf + mp;
1056 i1 = i0 + 1;
1057 j0 = nbuf + mj * Np2;
1058 j1 = j0 + 1;
1059 f0 = fbuf;
1060 f1 = f0 + 1;
1061 for (i = 0; i <= N2;
1062 i++, f++, i0+=2, i1+=2, f0+=2, f1+=2, j0+=2, j1+=2) {
1063 /*
1064 * ii0 = 2 * i;
1065 * ii1 = ii0 + 1;
1066 *
1067 * ssum[i] -= nbuf[mp + ii0] * nbuf[mp + ii0];
1068 * ssum[i] -= nbuf[mp + ii1] * nbuf[mp + ii1];
1069 * ssum[i] += fbuf[ii0] * fbuf[ii0];
1070 * ssum[i] += fbuf[ii1] * fbuf[ii1];
1071 */
1072 *f -= *i0 * *i0;
1073 *f -= *i1 * *i1;
1074 *f += *f0 * *f0;
1075 *f += *f1 * *f1;
1076 avg = minv * *f; /* avg = minv * ssum[i]; */
1077 if (avg < FL(0.0))
1078 avg = FL(0.0);
1079 if (avg == FL(0.0))
1080 fac = FL(0.0);
1081 else
1082 fac = avg / (avg + nref[i]);
1083 for (j = 1; j < sh; j++)
1084 fac *= fac;
1085 gain = g0m * fac + g0;
1086 /*
1087 * nbuf[mp + ii0] = fbuf[ii0];
1088 * nbuf[mp + ii1] = fbuf[ii1];
1089 * fbuf[ii0] = gain * nbuf[mj*Np2 + ii0];
1090 * fbuf[ii1] = gain * nbuf[mj*Np2 + ii1];
1091 */
1092 *i0 = *f0;
1093 *i1 = *f1;
1094 *f0 = gain * *j0;
1095 *f1 = gain * *j1;
1096 }
1097 }
1098
1099 if (chan == (Chans - 1)) {
1100 if (++mi >= m)
1101 mi = 0;
1102 if (++mj >= m)
1103 mj = 0;
1104 mp = mi * Np2;
1105 }
1106
1107 /* synthesis: The synthesis subroutine uses the Weighted Overlap-Add
1108 technique to reconstruct the time-domain signal. The (N/2 + 1)
1109 phase vocoder channel outputs at time n are inverse Fourier
1110 transformed, windowed, and added into the output array. */
1111
1112 fsst2(csound, fftsetup_inv, fbuf);
1113 //fsst(csound, fbuf, N);
1114
1115 lk = nO - (int64_t) sLen - 1; /*time shift*/
1116 while (lk < 0)
1117 lk += (int64_t) N;
1118 k = (int32_t) (lk % (int64_t) N);
1119
1120 f = fbuf + k;
1121 w = sWin - sLen;
1122 for (i = -sLen; i <= sLen; i++, k++, f++, w++) {
1123 obp += Chans;
1124 obc += Chans;
1125 if (obc >= obuflen) {
1126 obc = chan;
1127 obp = ob2 + chan;
1128 }
1129 if (k >= N) {
1130 k = 0;
1131 f = fbuf;
1132 }
1133 *obp += *w * *f;
1134 }
1135
1136 if (flag) {
1137 if (nread < ibuflen) { /* EOF detected */
1138 flag = 0;
1139 if ((lnread / Chans) < nMax)
1140 nMax = (lnread / Chans);
1141 }
1142 }
1143
1144 }
1145
1146 ibs += (Chans * D); /* starting point in ibuf */
1147 obs += (Chans * I); /* starting point in obuf */
1148
1149 nI += (int64_t) D; /* increment time */
1150 nO += (int64_t) I;
1151
1152 if (Verbose) {
1153 nImodR += D;
1154 if (nImodR > (int64_t) R) {
1155 nImodR -= (int64_t) R;
1156 csound->Message(csound,
1157 Str("%5.1f seconds of input complete\n"),(time+D*invR));
1158 }
1159 }
1160
1161 }
1162
1163 nMaxOut = (int64_t) (nMax * Chans);
1164 i = (int32_t) (nMaxOut - oCnt);
1165 if (i > obuflen) {
1166 writebuffer(csound, outfd, ob1, obuflen, &nrecs, &O);
1167 i -= obuflen;
1168 ob1 = ob2;
1169 }
1170 if (i > 0)
1171 writebuffer(csound, outfd, ob1, i, &nrecs, &O);
1172
1173 /* csound->rewriteheader(outfd); */
1174 csound->Message(csound, "\n\n");
1175 if (Verbose) {
1176 csound->Message(csound, "%s", Str("processing complete\n"));
1177 csound->Message(csound, "N = %d\n", N);
1178 csound->Message(csound, "M = %d\n", M);
1179 csound->Message(csound, "L = %d\n", L);
1180 csound->Message(csound, "D = %d\n", D);
1181 }
1182 return 0;
1183 }
1184
1185 static const char *usage_txt[] = {
1186 Str_noop("usage: dnoise [flags] input_file"),
1187 "",
1188 Str_noop("flags:"),
1189 Str_noop("i = noise reference soundfile"),
1190 Str_noop("o = output file"),
1191 Str_noop("N = # of bandpass filters (1024)"),
1192 Str_noop("w = filter overlap factor: {0,1,(2),3} DO NOT USE -w AND -M"),
1193 Str_noop("M = analysis window length (N-1 unless -w is specified)"),
1194 Str_noop("L = synthesis window length (M)"),
1195 Str_noop("D = decimation factor (M/8)"),
1196 Str_noop("b = begin time in noise reference soundfile (0)"),
1197 Str_noop("B = starting sample in noise reference soundfile (0)"),
1198 Str_noop("e = end time in noise reference soundfile (end)"),
1199 Str_noop("E = final sample in noise reference soundfile (end)"),
1200 Str_noop("t = threshold above noise reference in dB (30)"),
1201 Str_noop("S = sharpness of noise-gate turnoff (1) (1 to 5)"),
1202 Str_noop("n = number of FFT frames to average over (5)"),
1203 Str_noop("m = minimum gain of noise-gate when off in dB (-40)"),
1204 Str_noop("V : verbose - print status info"),
1205 Str_noop("A : AIFF format output"),
1206 Str_noop("W : WAV format output"),
1207 Str_noop("J : IRCAM format output"),
1208 NULL
1209 };
1210
dnoise_usage(CSOUND * csound,int32_t exitcode)1211 static int32_t dnoise_usage(CSOUND *csound, int32_t exitcode)
1212 {
1213 const char **sp;
1214
1215 for (sp = &(usage_txt[0]); *sp != NULL; sp++)
1216 csound->Message(csound, "%s\n", Str(*sp));
1217
1218 return exitcode;
1219 }
1220
1221 /* report soundfile write(osfd) error */
1222 /* called after chk of write() bytecnt */
1223
sndwrterr(CSOUND * csound,int32_t nret,int32_t nput)1224 static void sndwrterr(CSOUND *csound, int32_t nret, int32_t nput)
1225 {
1226 csound->Message(csound, Str("soundfile write returned sample count of %d, "
1227 "not %d\n"), nret, nput);
1228 csound->Message(csound, "%s", Str("(disk may be full...\n"
1229 " closing the file ...)\n"));
1230 /* FIXME: should clean up */
1231 //csound->Die(csound, "%s", Str("\t... closed\n"));
1232 }
1233
writebuffer(CSOUND * csound,SNDFILE * outfd,MYFLT * outbuf,int32_t nsmps,int32_t * nrecs,OPARMS * O)1234 static int32_t writebuffer(CSOUND *csound, SNDFILE *outfd,
1235 MYFLT *outbuf, int32_t nsmps, int32_t *nrecs, OPARMS *O)
1236 {
1237 int32_t n;
1238
1239 if (UNLIKELY(outfd == NULL)) return 0;
1240 n = sf_write_MYFLT(outfd, outbuf, nsmps);
1241 if (UNLIKELY(n < nsmps)) {
1242 sf_close(outfd);
1243 sndwrterr(csound, n, nsmps);
1244 return -1;
1245 }
1246 if (UNLIKELY(O->rewrt_hdr))
1247 csound->rewriteheader(outfd);
1248
1249 (*nrecs)++; /* JPff fix */
1250 switch (O->heartbeat) {
1251 case 1:
1252 csound->MessageS(csound, CSOUNDMSG_REALTIME, "%c\b", "|/-\\"[*nrecs & 3]);
1253 break;
1254 case 2:
1255 csound->MessageS(csound, CSOUNDMSG_REALTIME, ".");
1256 break;
1257 case 3:
1258 csound->MessageS(csound, CSOUNDMSG_REALTIME, "%d%n", *nrecs, &n);
1259 while (n--) csound->MessageS(csound, CSOUNDMSG_REALTIME, "\b");
1260 break;
1261 case 4:
1262 csound->MessageS(csound, CSOUNDMSG_REALTIME, "\a");
1263 break;
1264 }
1265
1266 return nsmps;
1267 }
1268
hamming(MYFLT * win,int32_t winLen,int32_t even)1269 static void hamming(MYFLT *win, int32_t winLen, int32_t even)
1270 {
1271 double ftmp;
1272 int32_t i;
1273
1274 ftmp = PI / winLen;
1275
1276 if (even) {
1277 for (i = 0; i < winLen; i++)
1278 win[i] = (MYFLT) (0.54 + 0.46 * cos(ftmp * ((double)i + 0.5)));
1279 win[winLen] = FL(0.0);
1280 }
1281 else {
1282 win[0] = FL(1.0);
1283 for (i = 1; i <= winLen; i++)
1284 win[i] = (MYFLT) (0.54 + 0.46 * cos(ftmp * (double)i));
1285 }
1286 }
1287
1288 /* module interface */
1289
dnoise_init_(CSOUND * csound)1290 int32_t dnoise_init_(CSOUND *csound)
1291 {
1292 int32_t retval = csound->AddUtility(csound, "dnoise", dnoise);
1293 if (!retval) {
1294 retval =
1295 csound->SetUtilityDescription(csound, "dnoise",
1296 Str("Removes noise from a sound file"));
1297 }
1298 return retval;
1299 }
1300