1 /*
2 atsa.c:
3
4 ATS analysis utility
5 Copyright (C) 2002-2004 Oscar Pablo Di Liscia, Pete Moss, Juan Pampin
6 Ported to Csound by Istvan Varga, original version is available at
7 http://sourceforge.net/projects/atsa/
8
9 This file is part of Csound.
10
11 The Csound Library is free software; you can redistribute it
12 and/or modify it under the terms of the GNU Lesser General Public
13 License as published by the Free Software Foundation; either
14 version 2.1 of the License, or (at your option) any later version.
15
16 Csound is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 GNU Lesser General Public License for more details.
20
21 You should have received a copy of the GNU Lesser General Public
22 License along with Csound; if not, write to the Free Software
23 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
24 02110-1301 USA
25 */
26
27 #define _FILE_OFFSET_BITS 64
28
29 #include "std_util.h"
30
31 #include <math.h>
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <errno.h>
36
37 #include <sndfile.h>
38
39 #if defined(__GNUC__) && defined(__STRICT_ANSI__)
40 # ifndef inline
41 # define inline __inline__
42 # endif
43 #endif
44
45 typedef float mus_sample_t;
46
47 /* window types */
48 #define BLACKMAN 0
49 #define BLACKMAN_H 1
50 #define HAMMING 2
51 #define VONHANN 3
52
53 /********** ANALYSIS PARAMETERS ***********/
54 /* start time */
55 #define ATSA_START 0.0f
56 /* duration */
57 #define ATSA_DUR 0.0f
58 /* lowest frequency (hertz) */
59 #define ATSA_LFREQ 20.0f
60 /* highest frequency (hertz) */
61 #define ATSA_HFREQ 20000.0f
62 /* frequency deviation (ratio) */
63 #define ATSA_FREQDEV 0.1f
64 /* number of f0 cycles in window */
65 #define ATSA_WCYCLES 4
66 /* window type */
67 #define ATSA_WTYPE BLACKMAN_H
68 /* window size */
69 #define ATSA_WSIZE 1024
70 /* hop size proportional to window size (ratio) */
71 #define ATSA_HSIZE 0.25f
72 /* lowest magnitude for peaks (amp) */
73 #define ATSA_LMAG -60.0f
74 /* length of analysis tracks (frames) */
75 #define ATSA_TRKLEN 3
76 /* minimum short partial length (frames) */
77 #define ATSA_MSEGLEN 3
78 /* minimum short partial SMR avreage (dB SPL) */
79 #define ATSA_MSEGSMR 60.0f
80 /* minimum gap length (frames) */
81 #define ATSA_MGAPLEN 3
82 /* threshold for partial SMR average (dB SPL) */
83 #define ATSA_SMRTHRES 30.0f
84 /* last peak contribution for tracking (ratio) */
85 #define ATSA_LPKCONT 0.0f
86 /* SMR contribution for tracking (ratio) */
87 #define ATSA_SMRCONT 0.5f
88 /* minimum number of frames for analysis (frames) */
89 #define ATSA_MFRAMES 4
90 /* default analysis file type
91 * 1 =only amp. and freq.
92 * 2 =amp., freq. and phase
93 * 3 =amp., freq. and noise
94 * 4 =amp., freq., phase, and noise
95 */
96 #define ATSA_TYPE 4
97 /* default residual file */
98 #if defined(LINUX) || defined(MACOSX)
99 # define ATSA_RES_FILE "/tmp/atsa_res.wav"
100 #else
101 # define ATSA_RES_FILE "/atsa_res.wav"
102 #endif
103
104 /* constants and macros */
105 #define NIL (-1)
106 #define ATSA_MAX_DB_SPL (100.0)
107 #define ATSA_NOISE_THRESHOLD (-120)
108 #define ATSA_CRITICAL_BANDS (25)
109 #define ATSA_NOISE_VARIANCE (0.04)
110 /* array of critical band frequency edges base on data from:
111 * Zwicker, Fastl (1990) "Psychoacoustics Facts and Models",
112 * Berlin ; New York : Springer-Verlag
113 */
114 #define ATSA_CRITICAL_BAND_EDGES {0.0, 100.0, 200.0, 300.0, 400.0, 510.0, \
115 630.0, 770.0, 920.0, 1080.0, 1270.0, \
116 1480.0, 1720.0, 2000.0, 2320.0, 2700.0, \
117 3150.0, 3700.0, 4400.0, 5300.0, 6400.0, \
118 7700.0, 9500.0, 12000.0, 15500.0, 20000.0}
119
120 //#define AMP_DB(amp) ((amp) != 0.0 ? (float) log10((amp) * 20.0) : -32767.0f)
121 //#define DB_AMP(db) ((float) pow(10.0, (db) / 20.0))
122
123 /* data structures */
124
125 /* ANARGS
126 * ======
127 * analysis parameters
128 */
129 typedef struct {
130 /* args[0] is infile, args[1] is outfile */
131 char *args[2];
132 float start;
133 float duration;
134 float lowest_freq;
135 float highest_freq;
136 float freq_dev;
137 int win_cycles;
138 int win_type;
139 int win_size;
140 float hop_size;
141 float lowest_mag;
142 int track_len;
143 int min_seg_len;
144 int min_gap_len;
145 float last_peak_cont;
146 float SMR_cont;
147 float SMR_thres;
148 float min_seg_SMR;
149 /* parameters computed from command line */
150 int first_smp;
151 int cycle_smp;
152 int hop_smp;
153 int total_samps;
154 int srate;
155 int fft_size;
156 float fft_mag;
157 int lowest_bin;
158 int highest_bin;
159 int frames;
160 int type;
161 } ANARGS;
162
163 /* ATS_FFT
164 * fft data
165 */
166
167 typedef struct {
168 int size;
169 int rate;
170 MYFLT *data;
171 } ATS_FFT;
172
173 /* ATS_PEAK
174 * ========
175 * spectral peak data
176 */
177 typedef struct {
178 double amp;
179 double frq;
180 double pha;
181 double smr;
182 int track;
183 } ATS_PEAK;
184
185 /* ATS_FRAME
186 * =========
187 * analysis frame data
188 */
189 typedef struct {
190 ATS_PEAK *peaks;
191 int n_peaks;
192 double time;
193 } ATS_FRAME;
194
195 /* ATS_HEADER
196 * ==========
197 * ats file header data
198 */
199 typedef struct {
200 /* Magic Number for ID of file, must be 123.00 */
201 double mag;
202 /* sampling rate */
203 double sr;
204 /* Frame size (samples) */
205 double fs;
206 /* Window size (samples) */
207 double ws;
208 /* number of partials per frame */
209 double par;
210 /* number of frames present */
211 double fra;
212 /* max. amplitude */
213 double ma;
214 /* max. frequency */
215 double mf;
216 /* duration (secs) */
217 double dur;
218 /* type (1,2 3 or 4)
219 * 1 =only amp. and freq.
220 * 2 =amp., freq. and phase
221 * 3 =amp., freq. and noise
222 * 4 =amp., freq., phase, and noise
223 */
224 double typ;
225 } ATS_HEADER;
226
227 /* ATS_SOUND
228 * =========
229 * ats analysis data
230 */
231 typedef struct {
232 /* global sound info */
233 int srate;
234 int frame_size;
235 int window_size;
236 int partials;
237 int frames;
238 double dur;
239 /* info deduced from analysis */
240 /* we use optimised to keep the
241 * # of partials killed by optimisation
242 */
243 int optimized;
244 double ampmax;
245 double frqmax;
246 ATS_PEAK *av;
247 /* sinusoidal data */
248 /* all of these ** are accessed as [partial][frame] */
249 double **time;
250 double **frq;
251 double **amp;
252 double **pha;
253 double **smr;
254 /* noise data */
255 int *bands;
256 double **res;
257 double **band_energy;
258 } ATS_SOUND;
259
260 /* Interface:
261 * ==========
262 * grouped by file in alphabetical order
263 */
264
265 /* atsa.c */
266
267 /* main_anal
268 * =========
269 * main analysis function
270 * soundfile: path to input file
271 * out_file: path to output ats file
272 * anargs: pointer to analysis parameters
273 * resfile: path to residual file
274 * returns error status
275 */
276 static int main_anal(CSOUND *csound, char *soundfile, char *ats_outfile,
277 ANARGS *anargs, char *resfile);
278
279 /* critical-bands.c */
280
281 /* evaluate_smr
282 * ============
283 * evalues the masking curves of an analysis frame
284 * peaks: pointer to an array of peaks
285 * peaks_size: number of peaks
286 */
287 static void evaluate_smr(ATS_PEAK *peaks, int peaks_size);
288
289 /* other-utils.c */
290
291 /* window_norm
292 * ===========
293 * computes the norm of a window
294 * returns the norm value
295 * win: pointer to a window
296 * size: window size
297 */
298 static float window_norm(float *win, int size);
299
300 /* make_window
301 * ===========
302 * makes an analysis window, returns a pointer to it.
303 * win_type: window type, available types are:
304 * BLACKMAN, BLACKMAN_H, HAMMING and VONHANN
305 * win_size: window size
306 */
307 static float *make_window(CSOUND *csound, int win_type, int win_size);
308
309 /* push_peak
310 * =========
311 * pushes a peak into an array of peaks
312 * re-allocating memory and updating its size
313 * returns a pointer to the array of peaks.
314 * new_peak: pointer to new peak to push into the array
315 * peaks_list: list of peaks
316 * peaks_size: pointer to the current size of the array.
317 */
318 static ATS_PEAK *push_peak(CSOUND *csound, ATS_PEAK *new_peak,
319 ATS_PEAK *peaks, int *peaks_size);
320
321 /* peak_frq_inc
322 * ============
323 * function used by qsort to sort an array of peaks
324 * in increasing frequency order.
325 */
326 static int peak_frq_inc(void const *a, void const *b);
327
328 /* peak_amp_inc
329 * ============
330 * function used by qsort to sort an array of peaks
331 * in increasing amplitude order.
332 */
333 static int peak_amp_inc(void const *a, void const *b);
334
335 #if 0
336 /* peak_smr_dec
337 * ============
338 * function used by qsort to sort an array of peaks
339 * in decreasing SMR order.
340 */
341 static int peak_smr_dec(void const *a, void const *b);
342 #endif
343
344 /* peak-detection.c */
345
346 /* peak_detection
347 * ==============
348 * detects peaks in a ATS_FFT block
349 * returns an array of detected peaks.
350 * ats_fft: pointer to ATS_FFT structure
351 * lowest_bin: lowest fft bin to start detection
352 * highest_bin: highest fft bin to end detection
353 * lowest_mag: lowest magnitude to detect peaks
354 * norm: analysis window norm
355 * peaks_size: pointer to size of the returned peaks array
356 */
357 static ATS_PEAK *peak_detection(CSOUND *csound, ATS_FFT *ats_fft,
358 int lowest_bin, int highest_bin,
359 float lowest_mag, double norm,
360 int *peaks_size);
361
362 /* peak-tracking.c */
363
364 /* peak_tracking
365 * =============
366 * connects peaks from one analysis frame to tracks
367 * returns a pointer to the analysis frame.
368 * tracks: pointer to the tracks
369 * tracks_size: numeber of tracks
370 * peaks: peaks to connect
371 * peaks_size: number of peaks
372 * frq_dev: frequency deviation from tracks
373 * SMR_cont: contribution of SMR to tracking
374 * n_partials: pointer to the number of partials before tracking
375 */
376 static ATS_FRAME *peak_tracking(CSOUND *csound, ATS_PEAK *tracks,
377 int *tracks_size, ATS_PEAK *peaks,
378 int *peaks_size, float frq_dev,
379 float SMR_cont, int *n_partials);
380
381 /* update_tracks
382 * =============
383 * updates analysis tracks
384 * returns a pointer to the tracks.
385 * tracks: pointer to the tracks
386 * tracks_size: numeber of tracks
387 * track_len: length of tracks
388 * frame_n: analysis frame number
389 * ana_frames: pointer to previous analysis frames
390 * last_peak_cont: contribution of last peak to the track
391 */
392 static ATS_PEAK *update_tracks(CSOUND *csound, ATS_PEAK *tracks,
393 int *tracks_size, int track_len, int frame_n,
394 ATS_FRAME *ana_frames, float last_peak_cont);
395
396 /* save-load-sound.c */
397
398 /* ats_save
399 * ========
400 * saves an ATS_SOUND to disk.
401 * sound: pointer to ATS_SOUND structure
402 * outfile: pointer to output ats file
403 * SMR_thres: partials with and avreage SMR
404 * below this value are considered masked
405 * and not written out to the ats file
406 * type: file type
407 * NOTE: sound MUST be optimised using optimize_sound
408 * before calling this function
409 */
410 static void ats_save(CSOUND *csound, ATS_SOUND *sound, FILE *outfile,
411 float SMR_thres, int type);
412
413 /* tracker.c */
414
415 /* tracker
416 * =======
417 * partial tracking function
418 * returns an ATS_SOUND with data issued from analysis
419 * anargs: pointer to analysis parameters
420 * soundfile: path to input file
421 * resfile: path to residual file
422 */
423 static ATS_SOUND *tracker(CSOUND *csound, ANARGS *anargs, char *soundfile,
424 char *resfile);
425
426 /* utilities.c */
427
428 /* ppp2
429 * ====
430 * returns the closest power of two
431 * greater than num
432 */
433 static inline unsigned int ppp2(int num);
434
435 /* various conversion functions
436 * to deal with dB and dB SPL
437 * they take and return double floats
438 */
439 static inline double amp2db(double amp);
440 static inline double db2amp(double db);
441 static inline double amp2db_spl(double amp);
442 // static inline double db2amp_spl(double db_spl);
443
444 /* optimize_sound
445 * ==============
446 * optimises an ATS_SOUND in memory before saving
447 * anargs: pointer to analysis parameters
448 * sound: pointer to ATS_SOUND structure
449 */
450 static void optimize_sound(CSOUND *csound, ANARGS *anargs, ATS_SOUND *sound);
451
452 /* residual.c */
453
454 /* compute_residual
455 * ================
456 * Computes the difference between the synthesis and the original sound.
457 * the <win-samps> array contains the sample numbers in the input file
458 * corresponding to each frame
459 * fil: pointer to analysed data
460 * fil_len: length of data in samples
461 * output_file: output file path
462 * sound: pointer to ATS_SOUND
463 * win_samps: pointer to array of analysis windows center times
464 * file_sampling_rate: sampling rate of analysis file
465 */
466 static void compute_residual(CSOUND *csound, mus_sample_t **fil,
467 int fil_len, char *output_file,
468 ATS_SOUND *sound, int *win_samps,
469 int file_sampling_rate);
470
471 /* residual-analysis.c */
472
473 /* residual_analysis
474 * =================
475 * performs the critical-band analysis of the residual file
476 * file: name of the sound file containing the residual
477 * sound: sound to store the residual data
478 */
479 static void residual_analysis(CSOUND *csound, char *file, ATS_SOUND *sound);
480
481 #if 0
482 /* band_energy_to_res
483 * ==================
484 * transfers residual engergy from bands to partials
485 * sound: sound structure containing data
486 * frame: frame number
487 */
488 static void band_energy_to_res(CSOUND *csound, ATS_SOUND *sound, int frame);
489 #endif
490
491 #if 0
492 /* res_to_band_energy
493 * ==================
494 * transfers residual engergy from partials to bands
495 * sound: sound structure containing data
496 * frame: frame number
497 */
498 static void res_to_band_energy(ATS_SOUND *sound, int frame);
499 #endif
500
501 /* init_sound
502 * ==========
503 * initialises a new sound allocating memory
504 */
505 static void init_sound(CSOUND *csound, ATS_SOUND *sound, int sampling_rate,
506 int frame_size, int window_size, int frames,
507 double duration, int partials, int use_noise);
508
509 /* free_sound
510 * ==========
511 * frees sound's memory
512 */
513 static void free_sound(CSOUND *csound, ATS_SOUND *sound);
514
515 /* ------------------------------------------------------------------------ */
516
517 /* main_anal
518 * =========
519 * main analysis function
520 * soundfile: path to input file
521 * out_file: path to output ats file
522 * anargs: pointer to analysis parameters
523 * returns error status
524 */
main_anal(CSOUND * csound,char * soundfile,char * ats_outfile,ANARGS * anargs,char * resfile)525 static int main_anal(CSOUND *csound, char *soundfile, char *ats_outfile,
526 ANARGS *anargs, char *resfile)
527 {
528 /* create pointers and structures */
529 ATS_SOUND *sound = NULL;
530 FILE *outfile;
531 void *fd;
532
533 /* open output file */
534 fd = csound->FileOpen2(csound, &outfile, CSFILE_STD, ats_outfile, "wb",
535 NULL, CSFTYPE_ATS, 0);
536 if (UNLIKELY(fd == NULL)) {
537 csound->Die(csound, Str("\nCould not open %s for writing, %s\nbye...\n"),
538 ats_outfile, Str(sf_strerror(NULL)));
539 }
540 /* call tracker */
541 sound = tracker(csound, anargs, soundfile, resfile);
542 /* save sound */
543 if (LIKELY(sound != NULL)) {
544 csound->Message(csound, "%s", Str("saving ATS data..."));
545 ats_save(csound, sound, outfile, anargs->SMR_thres, anargs->type);
546 csound->Message(csound, "%s", Str("done!\n"));
547 }
548 else {
549 /* file I/O error */
550 return -2;
551 }
552 /* close output file */
553 csound->FileClose(csound, fd);
554 /* free ATS_SOUND memory */
555 free_sound(csound, sound);
556 return 0;
557 }
558
559 /* ------------------------------------------------------------------------ */
560
usage(CSOUND * csound)561 static CS_NOINLINE CS_NORETURN void usage(CSOUND *csound)
562 {
563 csound->Message(csound, "ATSA 1.0\n");
564 csound->Message(csound, "%s", Str("atsa soundfile atsfile [flags]\n"));
565 csound->Message(csound, "%s", Str("Flags:\n"));
566 csound->Message(csound, Str("\t -b start (%f seconds)\n"), ATSA_START);
567 csound->Message(csound, Str("\t -e duration (%f seconds or end)\n"),
568 ATSA_DUR);
569 csound->Message(csound, Str("\t -l lowest frequency (%f Hertz)\n"),
570 ATSA_LFREQ);
571 csound->Message(csound, Str("\t -H highest frequency (%f Hertz)\n"),
572 ATSA_HFREQ);
573 csound->Message(csound,
574 Str("\t -d frequency deviation (%f of partial freq.)\n"),
575 ATSA_FREQDEV);
576 csound->Message(csound, Str("\t -c window cycles (%d cycles)\n"),
577 ATSA_WCYCLES);
578 csound->Message(csound, Str("\t -w window type (type: %d)\n"), ATSA_WTYPE);
579 csound->Message(csound, "%s", Str("\t\t(Options: 0=BLACKMAN, 1=BLACKMAN_H, "
580 "2=HAMMING, 3=VONHANN)\n"));
581 csound->Message(csound, Str("\t -h hop size (%f of window size)\n"),
582 ATSA_HSIZE);
583 csound->Message(csound, Str("\t -m lowest magnitude (%f)\n"), ATSA_LMAG);
584 csound->Message(csound, Str("\t -t track length (%d frames)\n"),
585 ATSA_TRKLEN);
586 csound->Message(csound, Str("\t -s min. segment length (%d frames)\n"),
587 ATSA_MSEGLEN);
588 csound->Message(csound, Str("\t -g min. gap length (%d frames)\n"),
589 ATSA_MGAPLEN);
590 csound->Message(csound, Str("\t -T SMR threshold (%f dB SPL)\n"),
591 ATSA_SMRTHRES);
592 csound->Message(csound, Str("\t -S min. segment SMR (%f dB SPL)\n"),
593 ATSA_MSEGSMR);
594 csound->Message(csound, Str("\t -P last peak contribution "
595 "(%f of last peak's parameters)\n"),
596 ATSA_LPKCONT);
597 csound->Message(csound, Str("\t -M SMR contribution (%f)\n"), ATSA_SMRCONT);
598 csound->Message(csound, Str("\t -F File Type (type: %d)\n"), ATSA_TYPE);
599 csound->Message(csound, "%s", Str("\t\t(Options: 1=amp.and freq. only, "
600 "2=amp.,freq. and phase, "
601 "3=amp.,freq. and residual, "
602 "4=amp.,freq.,phase, and residual)\n\n"));
603 csound->LongJmp(csound, 1);
604 }
605
atsa_main(CSOUND * csound,int argc,char ** argv)606 static int atsa_main(CSOUND *csound, int argc, char **argv)
607 {
608 int i, val, end_of_flags = 0;
609 ANARGS *anargs;
610 char *soundfile = (char *) NULL, *ats_outfile = (char *) NULL;
611 char *s = (char *) NULL;
612 char cur_opt = '\0';
613
614 anargs = (ANARGS *) csound->Calloc(csound, sizeof(ANARGS));
615
616 /* default values for analysis args */
617 anargs->start = ATSA_START;
618 anargs->duration = ATSA_DUR;
619 anargs->lowest_freq = ATSA_LFREQ;
620 anargs->highest_freq = ATSA_HFREQ;
621 anargs->freq_dev = ATSA_FREQDEV;
622 anargs->win_cycles = ATSA_WCYCLES;
623 anargs->win_type = ATSA_WTYPE;
624 anargs->hop_size = ATSA_HSIZE;
625 anargs->lowest_mag = ATSA_LMAG;
626 anargs->track_len = ATSA_TRKLEN;
627 anargs->min_seg_len = ATSA_MSEGLEN;
628 anargs->min_gap_len = ATSA_MGAPLEN;
629 anargs->SMR_thres = ATSA_SMRTHRES;
630 anargs->min_seg_SMR = ATSA_MSEGSMR;
631 anargs->last_peak_cont = ATSA_LPKCONT;
632 anargs->SMR_cont = ATSA_SMRCONT;
633 anargs->type = ATSA_TYPE;
634
635 for (i = 1; i < argc; ++i) {
636 if (cur_opt == '\0') {
637 if (argv[i][0] != '-' || end_of_flags) {
638 if (soundfile == NULL)
639 soundfile = argv[i];
640 else if (ats_outfile == NULL)
641 ats_outfile = argv[i];
642 else
643 usage(csound);
644 continue;
645 }
646 else if (argv[i][1] == '-' && argv[i][2] == '\0') {
647 end_of_flags = 1;
648 continue;
649 }
650 else if (argv[i][1] == '\0')
651 usage(csound);
652 else {
653 cur_opt = argv[i][1];
654 s = &(argv[i][2]);
655 if (*s == '\0')
656 continue;
657 }
658 }
659 else
660 s = argv[i];
661 if (*s == '\0')
662 usage(csound);
663 switch (cur_opt) {
664 case 'b':
665 anargs->start = (float) atof(s);
666 break;
667 case 'e':
668 anargs->duration = (float) atof(s);
669 break;
670 case 'l':
671 anargs->lowest_freq = (float) atof(s);
672 break;
673 case 'H':
674 anargs->highest_freq = (float) atof(s);
675 break;
676 case 'd':
677 anargs->freq_dev = (float) atof(s);
678 break;
679 case 'c':
680 anargs->win_cycles = (int) atoi(s);
681 break;
682 case 'w':
683 anargs->win_type = (int) atoi(s);
684 break;
685 case 'h':
686 anargs->hop_size = (float) atof(s);
687 break;
688 case 'm':
689 anargs->lowest_mag = (float) atof(s);
690 break;
691 case 't':
692 anargs->track_len = (int) atoi(s);
693 break;
694 case 's':
695 anargs->min_seg_len = (int) atoi(s);
696 break;
697 case 'g':
698 anargs->min_gap_len = (int) atoi(s);
699 break;
700 case 'T':
701 anargs->SMR_thres = (float) atof(s);
702 break;
703 case 'S':
704 anargs->min_seg_SMR = (float) atof(s);
705 break;
706 case 'P':
707 anargs->last_peak_cont = (float) atof(s);
708 break;
709 case 'M':
710 anargs->SMR_cont = (float) atof(s);
711 break;
712 case 'F':
713 anargs->type = (int) atoi(s);
714 break;
715 default:
716 usage(csound);
717 }
718 cur_opt = '\0';
719 }
720 if (cur_opt != '\0' ||
721 soundfile == NULL || soundfile[0] == '\0' ||
722 ats_outfile == NULL || ats_outfile[0] == '\0')
723 usage(csound);
724 #ifdef WIN32
725 {
726 char buffer[160];
727 char * tmp = getenv("TEMP");
728 strNcpy(buffer, tmp, 160);
729 // MKG 2014 Jan 29: No linkage for strlcat with MinGW here.
730 // but wrong; corrected
731 //strlcat(buffer, ATSA_RES_FILE, 160);
732 strncat(buffer, ATSA_RES_FILE, 160-strlen(buffer)); buffer[159] = '\0';
733 val = main_anal(csound, soundfile, ats_outfile, anargs, buffer);
734 }
735 #else
736 val = main_anal(csound, soundfile, ats_outfile, anargs, ATSA_RES_FILE);
737 #endif
738 csound->Free(csound, anargs);
739 return (val);
740 }
741
742 /* ------------------------------------------------------------------------ */
743
744 /* private function prototypes */
745 static void clear_mask(ATS_PEAK *peaks, int peaks_size);
746 static double compute_slope_r(double val);
747 static double frq2bark(double frq, double *edges);
748 static int find_band(double frq, double *edges);
749
750 /* frq2bark
751 * ========
752 * frequency to bark scale conversion
753 */
frq2bark(double frq,double * edges)754 static double frq2bark(double frq, double *edges)
755 {
756 double lo_frq, hi_frq;
757 int band;
758
759 if (frq <= 400.0)
760 return (frq * 0.01);
761 if (UNLIKELY(frq >= 20000.0))
762 return (NIL);
763
764 band = find_band(frq, edges);
765 lo_frq = edges[band];
766 hi_frq = edges[band + 1];
767 return (1.0 + band + fabs(log10(frq / lo_frq) / log10(lo_frq / hi_frq)));
768 }
769
770 /* find_band
771 * =========
772 * returns the critical band number
773 * corresponding to frq
774 */
find_band(double frq,double * edges)775 static int find_band(double frq, double *edges)
776 {
777 int i = 0;
778
779 while (frq > edges[i++]);
780 return (i - 2);
781 }
782
783 /* compute_slope_r
784 * ===============
785 * computes masking curve's right slope from val
786 */
compute_slope_r(double val)787 static double compute_slope_r(double val)
788 {
789 double i = val - 40.0;
790
791 return (((i > 0.0) ? i : 0.0) * 0.37 - 27.0);
792 }
793
794 /* clear_mask
795 * ==========
796 * clears masking curves
797 * peaks: array of peaks representing the masking curve
798 * peaks_size: number of peaks in curve
799 */
clear_mask(ATS_PEAK * peaks,int peaks_size)800 static void clear_mask(ATS_PEAK *peaks, int peaks_size)
801 {
802 while (peaks_size--)
803 peaks[peaks_size].smr = 0.0;
804 }
805
806 /* evaluate_smr
807 * ============
808 * evalues the masking curves of an analysis frame
809 * setting the peaks smr slot.
810 * peaks: pointer to an array of peaks
811 * peaks_size: number of peaks
812 */
evaluate_smr(ATS_PEAK * peaks,int peaks_size)813 static void evaluate_smr(ATS_PEAK *peaks, int peaks_size)
814 {
815 double slope_l = -27.0, slope_r, delta_dB = -50.0;
816 double frq_masker, amp_masker, frq_maskee, amp_maskee, mask_term;
817 int i, j;
818 ATS_PEAK *maskee;
819 double edges[ATSA_CRITICAL_BANDS + 1] = ATSA_CRITICAL_BAND_EDGES;
820
821 clear_mask(peaks, peaks_size);
822 if (peaks_size == 1)
823 peaks[0].smr = amp2db_spl(peaks[0].amp);
824 else
825 for (i = 0; i < peaks_size; i++) {
826 maskee = &peaks[i];
827 frq_maskee = frq2bark(maskee->frq, edges);
828 amp_maskee = amp2db_spl(maskee->amp);
829 for (j = 0; j < peaks_size; j++)
830 if (i != j) {
831 frq_masker = frq2bark(peaks[j].frq, edges);
832 amp_masker = amp2db_spl(peaks[j].amp);
833 slope_r = compute_slope_r(amp_masker);
834 mask_term = (frq_masker < frq_maskee) ?
835 (amp_masker + delta_dB +
836 (slope_r * (frq_maskee - frq_masker))) : (amp_masker +
837 delta_dB +
838 (slope_l *
839 (frq_masker -
840 frq_maskee)));
841 if (mask_term > maskee->smr)
842 maskee->smr = mask_term;
843 }
844 maskee->smr = amp_maskee - maskee->smr;
845 }
846 }
847
848 /* ------------------------------------------------------------------------ */
849
850 /* make_window
851 * ===========
852 * makes an analysis window, returns a pointer to it.
853 * win_type: window type, available types are:
854 * BLACKMAN, BLACKMAN_H, HAMMING and VONHANN
855 * win_size: window size
856 */
make_window(CSOUND * csound,int win_type,int win_size)857 static float *make_window(CSOUND *csound, int win_type, int win_size)
858 {
859 float *buffer;
860 int i;
861 float arg = TWOPI / (float) (win_size - 1);
862
863 buffer = (float *) csound->Malloc(csound, win_size * sizeof(float));
864
865 for (i = 0; i < win_size; i++) {
866 switch (win_type) {
867 case BLACKMAN: /* Blackman (3 term) */
868 buffer[i] = (float)(0.42 - 0.5 * cos(arg * i) + 0.08 * cos(arg * (i+i)));
869 break;
870 case BLACKMAN_H: /* Blackman-Harris (4 term) */
871 buffer[i] =(float)(
872 0.35875 - 0.48829 * cos(arg * i) + 0.14128 * cos(arg * (i+i)) -
873 0.01168 * cos(arg * (i+i+i)));
874 break;
875 case HAMMING: /* Hamming */
876 buffer[i] = (float)(0.54 - 0.46 * cos(arg * i));
877 break;
878 case VONHANN: /* Von Hann ("hanning") */
879 buffer[i] = (float)(0.5 - 0.5 * cos(arg * i));
880 break;
881 }
882 }
883 return (buffer);
884 }
885
886 /* window_norm
887 * ===========
888 * computes the norm of a window
889 * returns the norm value
890 * win: pointer to a window
891 * size: window size
892 */
window_norm(float * win,int size)893 static float window_norm(float *win, int size)
894 {
895 float acc = 0.0f;
896 int i;
897
898 for (i = 0; i < size; i++) {
899 acc += win[i];
900 }
901 return (2.0f / acc);
902 }
903
904 /* push_peak
905 * =========
906 * pushes a peak into an array of peaks
907 * re-allocating memory and updating its size
908 * returns a pointer to the array of peaks.
909 * new_peak: pointer to new peak to push into the array
910 * peaks_list: list of peaks
911 * peaks_size: pointer to the current size of the array.
912 */
push_peak(CSOUND * csound,ATS_PEAK * new_peak,ATS_PEAK * peaks_list,int * peaks_size)913 static ATS_PEAK *push_peak(CSOUND *csound, ATS_PEAK *new_peak,
914 ATS_PEAK *peaks_list, int *peaks_size)
915 {
916 peaks_list =
917 (ATS_PEAK *) csound->ReAlloc(csound, peaks_list,
918 sizeof(ATS_PEAK) * ++*peaks_size);
919 peaks_list[*peaks_size - 1] = *new_peak;
920 return (peaks_list);
921 }
922
923 /* peak_frq_inc
924 * ============
925 * function used by qsort to sort an array of peaks
926 * in increasing frequency order.
927 */
peak_frq_inc(void const * a,void const * b)928 static int peak_frq_inc(void const *a, void const *b)
929 {
930 return (int)(1000.0 * (((ATS_PEAK *) a)->frq - ((ATS_PEAK *) b)->frq));
931 }
932
933 /* peak_amp_inc
934 * ============
935 * function used by qsort to sort an array of peaks
936 * in increasing amplitude order.
937 */
peak_amp_inc(void const * a,void const * b)938 static int peak_amp_inc(void const *a, void const *b)
939 {
940 return (int)(1000.0 * (((ATS_PEAK *) a)->amp - ((ATS_PEAK *) b)->amp));
941 }
942
943 #if 0
944 /* peak_smr_dec
945 * ============
946 * function used by qsort to sort an array of peaks
947 * in decreasing SMR order.
948 */
949 static int peak_smr_dec(void const *a, void const *b)
950 {
951 return (int)(1000.0 * (((ATS_PEAK *) b)->smr - ((ATS_PEAK *) a)->smr));
952 }
953 #endif
954
atsa_sound_read_noninterleaved(SNDFILE * sf,mus_sample_t ** bufs,int nChannels,int nFrames)955 static CS_NOINLINE void atsa_sound_read_noninterleaved(SNDFILE *sf,
956 mus_sample_t **bufs,
957 int nChannels,
958 int nFrames)
959 {
960 mus_sample_t tmpBuf[128];
961 int i, j, k, m, n;
962
963 m = 128 / nChannels;
964 k = m * nChannels; /* samples in tmpBuf[] */
965 j = k; /* position in tmpBuf[] */
966 for (i = 0; i < nFrames; i++) {
967 if (j >= k) {
968 if ((nFrames - i) < m) {
969 m = (nFrames - i);
970 k = m * nChannels;
971 }
972 if (sizeof(mus_sample_t) == sizeof(float))
973 n = (int) sf_readf_float(sf, (void *) &(tmpBuf[0]), (sf_count_t) m);
974 else
975 n = (int) sf_readf_double(sf, (void *) &(tmpBuf[0]), (sf_count_t) m);
976 if (n < 0)
977 n = 0;
978 n *= nChannels;
979 for (; n < k; n++)
980 tmpBuf[n] = (mus_sample_t) 0;
981 j = 0;
982 }
983 for (n = 0; n < nChannels; n++)
984 bufs[n][i] = tmpBuf[j++];
985 }
986 }
987
atsa_sound_write_noninterleaved(SNDFILE * sf,mus_sample_t ** bufs,int nChannels,int nFrames)988 static CS_NOINLINE void atsa_sound_write_noninterleaved(SNDFILE *sf,
989 mus_sample_t **bufs,
990 int nChannels,
991 int nFrames)
992 {
993 mus_sample_t tmpBuf[128];
994 int i, j, k, m, n;
995
996 m = 128 / nChannels;
997 k = m * nChannels; /* samples in tmpBuf[] */
998 j = 0; /* position in tmpBuf[] */
999 for (i = 0; i < nFrames; i++) {
1000 for (n = 0; n < nChannels; n++)
1001 tmpBuf[j++] = bufs[n][i];
1002 if (j >= k || i == (nFrames - 1)) {
1003 n = j / nChannels;
1004 if (sizeof(mus_sample_t) == sizeof(float))
1005 sf_writef_float(sf, (void *) &(tmpBuf[0]), (sf_count_t) n);
1006 else
1007 sf_writef_double(sf, (void *) &(tmpBuf[0]), (sf_count_t) n);
1008 j = 0;
1009 }
1010 }
1011 }
1012
1013 /* ------------------------------------------------------------------------ */
1014
1015 /* private function prototypes */
1016 static void parabolic_interp(double alpha, double beta, double gamma,
1017 double *offset, double *height);
1018 static double phase_interp(double PeakPhase, double OtherPhase, double offset);
1019 static void to_polar(ATS_FFT *ats_fft, double *mags, double *phase, int N,
1020 double norm);
1021
1022 /* peak_detection
1023 * ==============
1024 * detects peaks in a ATS_FFT block
1025 * returns pointer to an array of detected peaks.
1026 * ats_fft: pointer to ATS_FFT structure
1027 * lowest_bin: lowest fft bin to start detection
1028 * highest_bin: highest fft bin to end detection
1029 * lowest_mag: lowest magnitude to detect peaks
1030 * norm: analysis window norm
1031 * peaks_size: pointer to size of the returned peaks array
1032 */
peak_detection(CSOUND * csound,ATS_FFT * ats_fft,int lowest_bin,int highest_bin,float lowest_mag,double norm,int * peaks_size)1033 static ATS_PEAK *peak_detection(CSOUND *csound, ATS_FFT *ats_fft,
1034 int lowest_bin, int highest_bin,
1035 float lowest_mag, double norm,
1036 int *peaks_size)
1037 {
1038 int k, N = (highest_bin ? highest_bin : ats_fft->size / 2);
1039 int first_bin = (lowest_bin ? ((lowest_bin > 2) ? lowest_bin : 2) : 2);
1040 double fft_mag =
1041 ((double) ats_fft->rate / ats_fft->size), *fftmags, *fftphase;
1042 double right_bin, left_bin, central_bin, offset;
1043 ATS_PEAK ats_peak, *peaks = NULL;
1044
1045 lowest_mag = (float) db2amp(lowest_mag);
1046
1047 /* init peak */
1048 ats_peak.amp = 0.0;
1049 ats_peak.frq = 0.0;
1050 ats_peak.pha = 0.0;
1051 ats_peak.smr = 0.0;
1052
1053 fftmags = (double *) csound->Malloc(csound, N * sizeof(double));
1054 fftphase = (double *) csound->Malloc(csound, N * sizeof(double));
1055 /* convert spectrum to polar coordinates */
1056 to_polar(ats_fft, fftmags, fftphase, N, norm);
1057 central_bin = fftmags[first_bin - 2];
1058 right_bin = fftmags[first_bin - 1];
1059 /* peak detection */
1060 for (k = first_bin; k < N; k++) {
1061 left_bin = central_bin;
1062 central_bin = right_bin;
1063 right_bin = fftmags[k];
1064 if ((central_bin > (double) lowest_mag) && (central_bin > right_bin) &&
1065 (central_bin > left_bin)) {
1066 parabolic_interp(left_bin, central_bin, right_bin, &offset,
1067 &ats_peak.amp);
1068 ats_peak.frq = fft_mag * ((k - 1) + offset);
1069 ats_peak.pha =
1070 (offset < 0.0) ? phase_interp(fftphase[k - 2], fftphase[k - 1],
1071 fabs(offset))
1072 : phase_interp(fftphase[k - 1], fftphase[k],
1073 offset);
1074 ats_peak.track = -1;
1075 /* push peak into peaks list */
1076 peaks = push_peak(csound, &ats_peak, peaks, peaks_size);
1077 }
1078 }
1079 /* free up fftmags and fftphase */
1080 csound->Free(csound, fftmags);
1081 csound->Free(csound, fftphase);
1082 return (peaks);
1083 }
1084
1085 /* to_polar
1086 * ========
1087 * rectangular to polar conversion
1088 * values are also scaled by window norm
1089 * and stored into separate arrays of
1090 * magnitudes and phases.
1091 * ats_fft: pointer to ATS_FFT structure
1092 * mags: pointer to array of magnitudes
1093 * phase: pointer to array of phases
1094 * N: highest bin in fft data array
1095 * norm: window norm used to scale magnitudes
1096 */
to_polar(ATS_FFT * ats_fft,double * mags,double * phase,int N,double norm)1097 static void to_polar(ATS_FFT *ats_fft, double *mags, double *phase, int N,
1098 double norm)
1099 {
1100 int k;
1101 double x, y;
1102
1103 for (k = 0; k < N; k++) {
1104 x = (double) ats_fft->data[k << 1];
1105 y = (double) ats_fft->data[(k << 1) + 1];
1106 mags[k] = norm * hypot(x, y);
1107 phase[k] = ((x == 0.0 && y == 0.0) ? 0.0 : atan2(y, x));
1108 }
1109 }
1110
1111 /* parabolic_interp
1112 * ================
1113 * parabolic peak interpolation
1114 */
parabolic_interp(double alpha,double beta,double gamma,double * offset,double * height)1115 static void parabolic_interp(double alpha, double beta, double gamma,
1116 double *offset, double *height)
1117 {
1118 double dbAlpha = amp2db(alpha), dbBeta = amp2db(beta), dbGamma = amp2db(gamma);
1119 *offset = 0.5 * ((dbAlpha - dbGamma) / (dbAlpha - 2.0 * dbBeta + dbGamma));
1120 *height = db2amp(dbBeta - ((dbAlpha - dbGamma) * 0.25 * *offset));
1121 }
1122
1123 /* phase_interp
1124 * ============
1125 * phase interpolation
1126 */
phase_interp(double PeakPhase,double RightPhase,double offset)1127 static double phase_interp(double PeakPhase, double RightPhase, double offset)
1128 {
1129 if ((PeakPhase - RightPhase) > PI * 1.5)
1130 RightPhase += TWOPI;
1131 else if ((RightPhase - PeakPhase) > PI * 1.5)
1132 PeakPhase += TWOPI;
1133 return ((RightPhase - PeakPhase) * offset + PeakPhase);
1134 }
1135
1136 /* ------------------------------------------------------------------------ */
1137
1138 /* private types */
1139 typedef struct {
1140 int size;
1141 ATS_PEAK *cands;
1142 } ATS_CANDS;
1143
1144 /* private function prototypes */
1145 static ATS_PEAK *find_candidates(CSOUND *csound, ATS_PEAK *peaks,
1146 int peaks_size, double lo, double hi,
1147 int *cand_size);
1148 static void sort_candidates(ATS_CANDS *cands, ATS_PEAK peak, float SMR_cont);
1149
1150 /* peak_tracking
1151 * =============
1152 * connects peaks from one analysis frame to tracks
1153 * returns a pointer to two frames of orphaned peaks.
1154 * tracks: pointer to the tracks
1155 * tracks_size: numeber of tracks
1156 * peaks: peaks to connect
1157 * peaks_size: number of peaks
1158 * frq_dev: frequency deviation from tracks
1159 * SMR_cont: contribution of SMR to tracking
1160 * n_partials: pointer to the number of partials before tracking
1161 */
peak_tracking(CSOUND * csound,ATS_PEAK * tracks,int * tracks_size,ATS_PEAK * peaks,int * peaks_size,float frq_dev,float SMR_cont,int * n_partials)1162 static ATS_FRAME *peak_tracking(CSOUND *csound, ATS_PEAK *tracks,
1163 int *tracks_size, ATS_PEAK *peaks,
1164 int *peaks_size, float frq_dev,
1165 float SMR_cont, int *n_partials)
1166 {
1167 ATS_CANDS *track_candidates =
1168 (ATS_CANDS *) csound->Malloc(csound, *peaks_size * sizeof(ATS_CANDS));
1169 double lo, hi;
1170 int k, j, used, goback;
1171 ATS_FRAME *returned_peaks =
1172 (ATS_FRAME *) csound->Malloc(csound, 2 * sizeof(ATS_FRAME));
1173
1174 returned_peaks[0].peaks = returned_peaks[1].peaks = NULL;
1175 returned_peaks[0].n_peaks = returned_peaks[1].n_peaks = 0;
1176
1177 /* sort data to prepare for matching */
1178 qsort(tracks, *tracks_size, sizeof(ATS_PEAK), peak_frq_inc);
1179 qsort(peaks, *peaks_size, sizeof(ATS_PEAK), peak_frq_inc);
1180
1181 /* find candidates for each peak and set each peak to best candidate */
1182 for (k = 0; k < *peaks_size; k++) {
1183 /* find frq limits for candidates */
1184 lo = peaks[k].frq - (0.5 * peaks[k].frq * frq_dev);
1185 hi = peaks[k].frq + (0.5 * peaks[k].frq * frq_dev);
1186 /* get possible candidates */
1187 track_candidates[k].size = 0;
1188 track_candidates[k].cands =
1189 find_candidates(csound, tracks, *tracks_size, lo, hi,
1190 &track_candidates[k].size);
1191 if (track_candidates[k].size) {
1192 sort_candidates(&track_candidates[k], peaks[k], SMR_cont);
1193 peaks[k].track = track_candidates[k].cands[0].track;
1194 }
1195 }
1196
1197 /* compare adjacent peaks track numbers to insure unique track numbers */
1198 do {
1199 goback = 0;
1200 for (j = 0; j < (*peaks_size - 1); j++)
1201 if ((peaks[j].track == peaks[j + 1].track) && (peaks[j].track > -1)) {
1202 if (track_candidates[j].cands[0].amp >
1203 track_candidates[j + 1].cands[0].amp) {
1204 track_candidates[j].cands[0].amp = ATSA_HFREQ;
1205 qsort(track_candidates[j].cands, track_candidates[j].size,
1206 sizeof(ATS_PEAK), peak_amp_inc);
1207 if (track_candidates[j].cands[0].amp < ATSA_HFREQ) {
1208 peaks[j].track = track_candidates[j].cands[0].track;
1209 goback = 1;
1210 }
1211 else
1212 peaks[j].track = -1;
1213 }
1214 else {
1215 track_candidates[j + 1].cands[0].amp = ATSA_HFREQ;
1216 qsort(track_candidates[j + 1].cands, track_candidates[j + 1].size,
1217 sizeof(ATS_PEAK), peak_amp_inc);
1218 if (track_candidates[j + 1].cands[0].amp < ATSA_HFREQ)
1219 peaks[j + 1].track = track_candidates[j + 1].cands[0].track;
1220 else
1221 peaks[j + 1].track = -1;
1222 }
1223 }
1224 } while (goback);
1225
1226 /* by this point, all peaks will either have a unique track number, or -1
1227 now we need to take care of those left behind */
1228 for (k = 0; k < *peaks_size; k++)
1229 if (peaks[k].track == -1) {
1230 peaks[k].track = (*n_partials)++;
1231 returned_peaks[1].peaks =
1232 push_peak(csound, &peaks[k], returned_peaks[1].peaks,
1233 &returned_peaks[1].n_peaks);
1234 }
1235
1236 /* check for tracks that didnt get assigned */
1237 for (k = 0; k < *tracks_size; k++) {
1238 used = 0;
1239 for (j = 0; j < *peaks_size; j++)
1240 if (tracks[k].track == peaks[j].track) {
1241 used = 1;
1242 break;
1243 }
1244 if (!used)
1245 returned_peaks[0].peaks =
1246 push_peak(csound, &tracks[k], returned_peaks[0].peaks,
1247 &returned_peaks[0].n_peaks);
1248 }
1249
1250 for (k = 0; k < *peaks_size; k++)
1251 csound->Free(csound, track_candidates[k].cands);
1252 csound->Free(csound, track_candidates);
1253 return (returned_peaks);
1254 }
1255
1256 /* find_candidates
1257 * ===============
1258 * find candidates to continue a track form an array of peaks
1259 * returns a pointer to an array of candidates
1260 * peaks: pointer to array of peaks
1261 * peaks_size: number of peaks
1262 * lo: lowest frequency to consider candidates
1263 * hi: highest frequency to consider candidates
1264 * cand_size: pointer to the number of candidates returned
1265 */
find_candidates(CSOUND * csound,ATS_PEAK * peaks,int peaks_size,double lo,double hi,int * cand_size)1266 static ATS_PEAK *find_candidates(CSOUND *csound, ATS_PEAK *peaks,
1267 int peaks_size, double lo, double hi,
1268 int *cand_size)
1269 {
1270 int i;
1271 ATS_PEAK *cand_list = NULL;
1272
1273 for (i = 0; i < peaks_size; i++)
1274 if ((lo <= peaks[i].frq) && (peaks[i].frq <= hi))
1275 cand_list = push_peak(csound, &peaks[i], cand_list, cand_size);
1276
1277 return (cand_list);
1278 }
1279
1280 /* sort_candidates
1281 * ===================
1282 * sorts candidates from best to worst according to frequency and SMR
1283 * peak_candidates: pointer to an array of candidate peaks
1284 * peak: the peak we are matching
1285 * SMR_cont: contribution of SMR to the matching
1286 */
sort_candidates(ATS_CANDS * cands,ATS_PEAK peak,float SMR_cont)1287 static void sort_candidates(ATS_CANDS *cands, ATS_PEAK peak, float SMR_cont)
1288 {
1289 int i;
1290
1291 /* compute delta values and store them in cands.amp
1292 (dont worry, the candidate amps are useless otherwise!) */
1293 for (i = 0; i < cands->size; i++)
1294 cands->cands[i].amp =
1295 (fabs(cands->cands[i].frq - peak.frq) +
1296 (SMR_cont * fabs(cands->cands[i].smr - peak.smr))) / (SMR_cont + 1);
1297
1298 /* sort list by amp (increasing) */
1299 qsort(cands->cands, cands->size, sizeof(ATS_PEAK), peak_amp_inc);
1300 }
1301
1302 /* update_tracks
1303 * =============
1304 * updates analysis tracks
1305 * returns a pointer to the tracks.
1306 * tracks: pointer to the tracks
1307 * tracks_size: numeber of tracks
1308 * track_len: length of tracks
1309 * frame_n: analysis frame number
1310 * ana_frames: pointer to previous analysis frames
1311 * last_peak_cont: contribution of last peak to the track
1312 */
update_tracks(CSOUND * csound,ATS_PEAK * tracks,int * tracks_size,int track_len,int frame_n,ATS_FRAME * ana_frames,float last_peak_cont)1313 static ATS_PEAK *update_tracks(CSOUND *csound, ATS_PEAK *tracks,
1314 int *tracks_size, int track_len, int frame_n,
1315 ATS_FRAME *ana_frames, float last_peak_cont)
1316 {
1317 int frames, first_frame, track, g, i, k;
1318 double frq_acc, last_frq, amp_acc, last_amp, smr_acc, last_smr;
1319 int f, a, s;
1320 ATS_PEAK *l_peaks, *peak;
1321
1322 if (tracks != NULL) {
1323 frames = (frame_n < track_len) ? frame_n : track_len;
1324 first_frame = frame_n - frames;
1325 for (g = 0; g < *tracks_size; g++) {
1326 track = tracks[g].track;
1327 frq_acc = last_frq = amp_acc = last_amp = smr_acc = last_smr = 0.0;
1328 f = a = s = 0;
1329 for (i = first_frame; i < frame_n; i++) {
1330 l_peaks = ana_frames[i].peaks;
1331 peak = NULL;
1332 for (k = 0; k < ana_frames[i].n_peaks; k++)
1333 if (l_peaks[k].track == track) {
1334 peak = &l_peaks[k];
1335 break;
1336 }
1337 if (peak != NULL) {
1338 if (peak->frq > 0.0) {
1339 last_frq = peak->frq;
1340 frq_acc += peak->frq;
1341 f++;
1342 }
1343 if (peak->amp > 0.0) {
1344 last_amp = peak->amp;
1345 amp_acc += peak->amp;
1346 a++;
1347 }
1348 if (peak->smr > 0.0) {
1349 last_smr = peak->smr;
1350 smr_acc += peak->smr;
1351 s++;
1352 }
1353 }
1354 }
1355 if (f)
1356 tracks[g].frq =
1357 (last_peak_cont * last_frq) +
1358 ((1 - last_peak_cont) * (frq_acc / f));
1359 if (a)
1360 tracks[g].amp =
1361 (last_peak_cont * last_amp) +
1362 ((1 - last_peak_cont) * (amp_acc / a));
1363 if (s)
1364 tracks[g].smr =
1365 (last_peak_cont * last_smr) +
1366 ((1 - last_peak_cont) * (smr_acc / s));
1367 }
1368 }
1369 else
1370 for (g = 0; g < ana_frames[frame_n - 1].n_peaks; g++)
1371 tracks =
1372 push_peak(csound, &ana_frames[frame_n - 1].peaks[g], tracks,
1373 tracks_size);
1374
1375 return (tracks);
1376 }
1377
1378 /* ------------------------------------------------------------------------ */
1379
1380 #define ATSA_RES_MIN_FFT_SIZE 4096
1381 #define ATSA_RES_PAD_FACTOR 2
1382 #define MAG_SQUARED(re, im, norm) (norm * (re*re+im*im))
1383
1384 /* private function prototypes */
1385 static int residual_get_N(int M, int min_fft_size, int factor);
1386 static void residual_get_bands(double fft_mag, double *true_bands,
1387 int *limits, int bands);
1388 //static double residual_compute_time_domain_energy(ATS_FFT *fft_struct);
1389 static double residual_get_band_energy(int lo, int hi, ATS_FFT *fft_struct,
1390 double norm);
1391 static void residual_compute_band_energy(ATS_FFT *fft_struct,
1392 int *band_limits, int bands,
1393 double *band_energy, double norm);
1394
residual_get_N(int M,int min_fft_size,int factor)1395 static int residual_get_N(int M, int min_fft_size, int factor)
1396 {
1397 int def_size = factor * M;
1398
1399 while (def_size < min_fft_size)
1400 def_size = ppp2(def_size + 1);
1401 return (def_size);
1402 }
1403
residual_get_bands(double fft_mag,double * true_bands,int * limits,int bands)1404 static void residual_get_bands(double fft_mag, double *true_bands,
1405 int *limits, int bands)
1406 {
1407 int k;
1408
1409 for (k = 0; k < bands; k++)
1410 limits[k] = (int)floor(true_bands[k] / fft_mag);
1411 }
1412
1413 /* static double residual_compute_time_domain_energy(ATS_FFT *fft) */
1414 /* { */
1415 /* /\* Parseval's Theorem states: */
1416
1417 /* N-1 N-1 */
1418 /* sum(|x(n)^2|) = 1/N* sum (|X(k)|^2) */
1419 /* n=0 k=0 */
1420
1421 /* then we multiply the time domain energy by 1/2 */
1422 /* because we only compute frequency energy between */
1423 /* 0 Hz and Nyquist only (0 -> N/2) */
1424 /* *\/ */
1425 /* int n; */
1426 /* double sum = 0.0; */
1427
1428 /* for (n = 0; n < fft->size; n++) */
1429 /* sum += fabs((double) fft->data[n] * (double) fft->data[n]); */
1430 /* return (sum); */
1431 /* } */
1432
residual_get_band_energy(int lo,int hi,ATS_FFT * fft,double norm)1433 static double residual_get_band_energy(int lo, int hi, ATS_FFT *fft,
1434 double norm)
1435 {
1436 /* does 1/N * sum(re^2+im^2) within a band around <center>
1437 from <lo> lower bin to <hi> upper bin in <fft-struct> */
1438 int k;
1439 double sum = 0.0;
1440
1441 if (lo < 0)
1442 lo = 0;
1443 if (hi > fft->size / 2)
1444 hi = fft->size / 2; /* was (int)floor(fft->size * 0.5) */
1445 for (k = lo; k < hi; k++) {
1446 double re = (double) fft->data[k << 1];
1447 double im = (double) fft->data[(k << 1) + 1];
1448
1449 sum += MAG_SQUARED(re, im, norm);
1450 }
1451 return (sum / (double) fft->size);
1452 }
1453
residual_compute_band_energy(ATS_FFT * fft,int * band_limits,int bands,double * band_energy,double norm)1454 static void residual_compute_band_energy(ATS_FFT *fft, int *band_limits,
1455 int bands, double *band_energy,
1456 double norm)
1457 {
1458 /* loop through bands and evaluate energy
1459 we compute energy of one band as:
1460 (N-1)/2
1461 1/N * sum(|X(k)|^2)
1462 k=0
1463 N=fft size, K=bins in band */
1464 int b;
1465
1466 for (b = 0; b < bands - 1; b++)
1467 band_energy[b] =
1468 residual_get_band_energy(band_limits[b], band_limits[b + 1], fft,
1469 norm);
1470 }
1471
1472 /* residual_analysis
1473 * =================
1474 * performs the critical-band analysis of the residual file
1475 * file: name of the sound file containing the residual
1476 * sound: sound to store the residual data
1477 */
residual_analysis(CSOUND * csound,char * file,ATS_SOUND * sound)1478 static void residual_analysis(CSOUND *csound, char *file, ATS_SOUND *sound)
1479 {
1480 int file_sampling_rate, sflen, hop, M, N, frames, *band_limits;
1481 int M_2, st_pt, filptr, i, frame_n, k;
1482 double norm = 1.0, threshold, fft_mag, **band_arr = NULL, *band_energy;
1483 //double time_domain_energy = 0.0, freq_domain_energy = 0.0, sum = 0.0;
1484 double edges[ATSA_CRITICAL_BANDS + 1] = ATSA_CRITICAL_BAND_EDGES;
1485 ATS_FFT fft;
1486 SF_INFO sfinfo;
1487 mus_sample_t **bufs;
1488 SNDFILE *sf;
1489 void *fd;
1490
1491 memset(&sfinfo, 0, sizeof(SF_INFO));
1492 fd = csound->FileOpen2(csound, &sf, CSFILE_SND_R, file, &sfinfo, "SFDIR;SSDIR",
1493 CSFTYPE_UNKNOWN_AUDIO, 0);
1494 if (UNLIKELY(fd == NULL)) {
1495 csound->Die(csound, Str("atsa: error opening residual file '%s'"), file);
1496 }
1497 if (UNLIKELY(sfinfo.channels != 2)) {
1498 csound->Die(csound,
1499 Str("atsa: residual file has %d channels, must be stereo !"),
1500 (int) sfinfo.channels);
1501 }
1502 file_sampling_rate = sfinfo.samplerate;
1503 sflen = (int) sfinfo.frames;
1504 hop = sound->frame_size;
1505 M = sound->window_size;
1506 N = residual_get_N(M, ATSA_RES_MIN_FFT_SIZE, ATSA_RES_PAD_FACTOR);
1507 bufs = (mus_sample_t **) csound->Malloc(csound, 2 * sizeof(mus_sample_t *));
1508 bufs[0] =
1509 (mus_sample_t *) csound->Malloc(csound, sflen * sizeof(mus_sample_t));
1510 bufs[1] =
1511 (mus_sample_t *) csound->Malloc(csound, sflen * sizeof(mus_sample_t));
1512 fft.size = N;
1513 fft.rate = file_sampling_rate;
1514 fft.data = (MYFLT *) csound->Malloc(csound, (N + 2) * sizeof(MYFLT));
1515 threshold = /*AMP_DB*/(ATSA_NOISE_THRESHOLD);
1516 frames = sound->frames;
1517 fft_mag = (double) file_sampling_rate / (double) N;
1518 band_limits =
1519 (int *) csound->Malloc(csound, sizeof(int) * (ATSA_CRITICAL_BANDS + 1));
1520 residual_get_bands(fft_mag, edges, band_limits, ATSA_CRITICAL_BANDS + 1);
1521 band_arr = sound->band_energy;
1522 band_energy =
1523 (double *) csound->Malloc(csound, ATSA_CRITICAL_BANDS * sizeof(double));
1524
1525 M_2 = (int)floor(((double) M - 1.0) * 0.5);
1526 st_pt = N - M_2;
1527 filptr = M_2 * -1;
1528 /* read sound into memory */
1529 atsa_sound_read_noninterleaved(sf, bufs, 2, sflen);
1530
1531 for (frame_n = 0; frame_n < frames; frame_n++) {
1532 for (i = 0; i < (N + 2); i++) {
1533 fft.data[i] = (MYFLT) 0;
1534 }
1535 for (k = 0; k < M; k++) {
1536 if (filptr >= 0 && filptr < sflen)
1537 fft.data[(k + st_pt) % N] = (MYFLT) bufs[0][filptr];
1538 filptr++;
1539 }
1540 //smp = filptr - M_2 - 1;
1541 //time_domain_energy = residual_compute_time_domain_energy(&fft);
1542 /* take the fft */
1543 csound->RealFFTnp2(csound, fft.data, N);
1544 residual_compute_band_energy(&fft, band_limits, ATSA_CRITICAL_BANDS + 1,
1545 band_energy, norm);
1546 //sum = 0.0;
1547 //for (k = 0; k < ATSA_CRITICAL_BANDS; k++) {
1548 // sum += band_energy[k];
1549 //}
1550 //freq_domain_energy = 2.0 * sum;
1551 for (k = 0; k < ATSA_CRITICAL_BANDS; k++) {
1552 if (band_energy[k] < threshold) {
1553 band_arr[k][frame_n] = 0.0;
1554 }
1555 else {
1556 band_arr[k][frame_n] = band_energy[k];
1557 }
1558 }
1559 filptr = filptr - M + hop;
1560 }
1561 /* save data in sound */
1562 sound->band_energy = band_arr;
1563 csound->Free(csound, fft.data);
1564 csound->Free(csound, band_energy);
1565 csound->Free(csound, band_limits);
1566 csound->Free(csound, bufs[0]);
1567 csound->Free(csound, bufs[1]);
1568 csound->Free(csound, bufs);
1569 }
1570
1571 #if 0
1572 /* band_energy_to_res
1573 * ==================
1574 * transfers residual engergy from bands to partials
1575 * sound: sound structure containing data
1576 * frame: frame number
1577 */
1578 static void band_energy_to_res(CSOUND *csound, ATS_SOUND *sound, int frame)
1579 {
1580 int i, j;
1581 double edges[] = ATSA_CRITICAL_BAND_EDGES;
1582 double bandsum[ATSA_CRITICAL_BANDS];
1583 double partialfreq, partialamp;
1584 double *partialbandamp; /* amplitude of the band that the partial is in */
1585 int *bandnum; /* the band number that the partial is in */
1586
1587 partialbandamp = csound->Malloc(csound, sizeof(double) * sound->partials);
1588 bandnum = csound->Malloc(csound, sizeof(int) * sound->partials);
1589 /* initialise the sum per band */
1590 for (i = 0; i < ATSA_CRITICAL_BANDS; i++)
1591 bandsum[i] = 0;
1592
1593 /* find find which band each partial is in */
1594 for (i = 0; i < sound->partials; i++) {
1595 partialfreq = sound->frq[i][frame];
1596 partialamp = sound->amp[i][frame];
1597 for (j = 0; j < 25; j++) {
1598 if ((partialfreq < edges[j + 1]) && (partialfreq >= edges[j])) {
1599 bandsum[j] += partialamp;
1600 bandnum[i] = j;
1601 partialbandamp[i] = sound->band_energy[j][frame];
1602 break;
1603 }
1604 }
1605 }
1606 /* compute energy per partial */
1607 for (i = 0; i < sound->partials; i++) {
1608 if (bandsum[bandnum[i]] > 0.0)
1609 sound->res[i][frame] =
1610 sound->amp[i][frame] * partialbandamp[i] / bandsum[bandnum[i]];
1611 else
1612 sound->res[i][frame] = 0.0;
1613 }
1614 csound->Free(csound, partialbandamp);
1615 csound->Free(csound, bandnum);
1616 }
1617 #endif
1618
1619 /* res_to_band_energy
1620 * ==================
1621 * transfers residual engergy from partials to bands
1622 * sound: sound structure containing data
1623 * frame: frame number
1624 */
1625 #if 0
1626 static void res_to_band_energy(ATS_SOUND *sound, int frame)
1627 {
1628 int j, par;
1629 double sum;
1630 double edges[ATSA_CRITICAL_BANDS + 1] = ATSA_CRITICAL_BAND_EDGES;
1631
1632 par = 0;
1633 for (j = 0; j < ATSA_CRITICAL_BANDS; j++) {
1634 sum = 0.0;
1635 while (sound->frq[par][frame] >= edges[j] &&
1636 sound->frq[par][frame] < edges[j + 1]) {
1637 sum += sound->res[par][frame];
1638 par++;
1639 }
1640 sound->band_energy[j][frame] = sum;
1641 }
1642 }
1643 #endif
1644
1645 /* ------------------------------------------------------------------------ */
1646
1647 /* private function prototypes */
1648 static int compute_m(double pha_1, double frq_1, double pha, double frq,
1649 int buffer_size);
1650 static double compute_aux(double pha_1, double pha, double frq_1,
1651 int buffer_size, int M);
1652 static double compute_alpha(double aux, double frq_1, double frq,
1653 int buffer_size);
1654 static double compute_beta(double aux, double frq_1, double frq,
1655 int buffer_size);
1656 static double interp_phase(double pha_1, double frq_1, double alpha,
1657 double beta, int i);
1658 static void read_frame(mus_sample_t **fil, int fil_len, int samp_1,
1659 int samp_2, double *in_buffer);
1660 static void synth_buffer(double a1, double a2, double f1, double f2,
1661 double p1, double p2, double *buffer,
1662 int frame_samps);
1663
1664 /* Functions for phase interpolation
1665 * All this comes from JOS/XJS article on PARSHL.
1666 * Original phase interpolation eqns. by Qualtieri/McAulay.
1667 */
1668
compute_m(double pha_1,double frq_1,double pha,double frq,int buffer_size)1669 static int compute_m(double pha_1, double frq_1, double pha, double frq,
1670 int buffer_size)
1671 {
1672 /* int val = (int) ((((pha_1 + (frq_1 * (double) buffer_size) - pha)
1673 + ((frq - frq_1) * 0.5 * (double) buffer_size)) / TWOPI)
1674 + 0.5); */
1675 return ((int)
1676 ((((pha_1 + (frq_1 * (double) buffer_size) - pha) +
1677 ((frq - frq_1) * 0.5 * (double) buffer_size)) / TWOPI) + 0.5));
1678 }
1679
compute_aux(double pha_1,double pha,double frq_1,int buffer_size,int M)1680 static double compute_aux(double pha_1, double pha, double frq_1,
1681 int buffer_size, int M)
1682 {
1683 /* double val = (double) ((pha + (TWOPI * (double) M))
1684 - (pha_1 + (frq_1 * (double) buffer_size))); */
1685 return ((double)
1686 ((pha + (TWOPI * (double) M)) -
1687 (pha_1 + (frq_1 * (double) buffer_size))));
1688 }
1689
compute_alpha(double aux,double frq_1,double frq,int buffer_size)1690 static double compute_alpha(double aux, double frq_1, double frq,
1691 int buffer_size)
1692 {
1693 /* double val = (double) (((3.0 / (double) (buffer_size * buffer_size)) * aux)
1694 - ((frq - frq_1) / (double) buffer_size)); */
1695 return ((double)
1696 (((3.0 / (double) (buffer_size * buffer_size)) * aux) -
1697 ((frq - frq_1) / (double) buffer_size)));
1698 }
1699
compute_beta(double aux,double frq_1,double frq,int buffer_size)1700 static double compute_beta(double aux, double frq_1, double frq,
1701 int buffer_size)
1702 {
1703 /* double val = (double) (((-2.0 / (double) (buffer_size * buffer_size
1704 * buffer_size)) * aux)
1705 + ((frq - frq_1)
1706 / (double) (buffer_size * buffer_size))); */
1707 return ((double)
1708 (((-2.0 / (double) (buffer_size * buffer_size * buffer_size)) *
1709 aux) + ((frq - frq_1) / (double) (buffer_size * buffer_size))));
1710 }
1711
interp_phase(double pha_1,double frq_1,double alpha,double beta,int i)1712 static double interp_phase(double pha_1, double frq_1, double alpha,
1713 double beta, int i)
1714 {
1715 /* double val = (double) ((beta * (double) (i * i * i))
1716 + (alpha * (double) (i * i))
1717 + (frq_1 * (double) i) + pha_1); */
1718 return ((double)
1719 ((beta * (double) (i * i * i)) + (alpha * (double) (i * i)) +
1720 (frq_1 * (double) i) + pha_1));
1721 }
1722
1723 /* read_frame
1724 * ==========
1725 * reads a frame from the input file
1726 * fil: pointer to an array with sound data
1727 * fil_len: length of datas in samples
1728 * samp_1: first sample number in frame
1729 * samp_2: last sample number in frame
1730 * in_buffer: pointer to input buffer
1731 * which is filled out by the function
1732 * NOTE: caller should allocate memory for buffer
1733 */
read_frame(mus_sample_t ** fil,int fil_len,int samp_1,int samp_2,double * in_buffer)1734 static void read_frame(mus_sample_t **fil, int fil_len, int samp_1,
1735 int samp_2, double *in_buffer)
1736 {
1737 int i, index, samps = samp_2 - samp_1;
1738 mus_sample_t tmp;
1739
1740 /* samps = samp_2 - samp_1; */
1741 for (i = 0; i < samps; i++) {
1742 index = samp_1 + i;
1743 if (index < fil_len)
1744 tmp = fil[0][index];
1745 else
1746 tmp = (mus_sample_t) 0.0;
1747 in_buffer[i] = (double) tmp;
1748 }
1749 }
1750
1751 /* synth_buffer
1752 * ============
1753 * synthesizes a buffer of sound using
1754 * amplitude linear interpolation and
1755 * phase cubic interpolation
1756 * a1: strating amplitude
1757 * a2: ending amplitude
1758 * f1: starting frequency in radians per sample
1759 * f2: ending frequency in radians per sample
1760 * p1: starting phase in radians
1761 * p2: ending phase in radians
1762 * buffer: pointer to synthsis buffer
1763 * which is filled out by the function
1764 * NOTE: caller should allocate memory for buffer
1765 * frame_samps: number of samples in frame (buffer)
1766 */
synth_buffer(double a1,double a2,double f1,double f2,double p1,double p2,double * buffer,int frame_samps)1767 static void synth_buffer(double a1, double a2, double f1, double f2,
1768 double p1, double p2, double *buffer,
1769 int frame_samps)
1770 {
1771 int k, M;
1772 double aux, alpha, beta, amp, amp_inc, int_pha;
1773
1774 M = compute_m(p1, f1, p2, f2, frame_samps);
1775 aux = compute_aux(p1, p2, f1, frame_samps, M);
1776 alpha = compute_alpha(aux, f1, f2, frame_samps);
1777 beta = compute_beta(aux, f1, f2, frame_samps);
1778 amp = a1;
1779 amp_inc = (a2 - a1) / (double) frame_samps;
1780 for (k = 0; k < frame_samps; k++) {
1781 int_pha = interp_phase(p1, f1, alpha, beta, k);
1782 buffer[k] += amp * cos(int_pha);
1783 amp += amp_inc;
1784 }
1785 }
1786
1787 /* compute_residual
1788 * ================
1789 * Computes the difference between the synthesis and the original sound.
1790 * the <win-samps> array contains the sample numbers in the input file
1791 * corresponding to each frame
1792 * fil: pointer to analysed data
1793 * fil_len: length of data in samples
1794 * output_file: output file path
1795 * sound: pointer to ATS_SOUND
1796 * win_samps: pointer to array of analysis windows center times
1797 * file_sampling_rate: sampling rate of analysis file
1798 */
compute_residual(CSOUND * csound,mus_sample_t ** fil,int fil_len,char * output_file,ATS_SOUND * sound,int * win_samps,int file_sampling_rate)1799 static void compute_residual(CSOUND *csound, mus_sample_t **fil,
1800 int fil_len, char *output_file,
1801 ATS_SOUND *sound, int *win_samps,
1802 int file_sampling_rate)
1803 {
1804 int i, frm, frm_1, frm_2, par, frames, partials, frm_samps, out_smp = 0;
1805 double *in_buff, *synth_buff, mag, a1, a2, f1, f2, p1, p2, diff, synth;
1806 mus_sample_t **obuf;
1807 SF_INFO sfinfo;
1808 SNDFILE *sf;
1809 void *fd;
1810
1811 frames = sound->frames;
1812 partials = sound->partials;
1813 frm_samps = sound->frame_size;
1814 mag = TWOPI / (double) file_sampling_rate;
1815 in_buff = (double *) csound->Malloc(csound, frm_samps * sizeof(double));
1816 synth_buff = (double *) csound->Malloc(csound, frm_samps * sizeof(double));
1817 /* open output file */
1818 memset(&sfinfo, 0, sizeof(SF_INFO));
1819 //sfinfo.frames = (sf_count_t)0; /* was -1 */
1820 sfinfo.samplerate = file_sampling_rate;
1821 sfinfo.channels = 2;
1822 sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
1823 fd = csound->FileOpen2(csound, &sf, CSFILE_SND_W, output_file, &sfinfo,
1824 NULL, CSFTYPE_WAVE, 0);
1825 if (UNLIKELY(fd == NULL)) {
1826 csound->Die(csound, Str("\nERROR: cannot open file %s for writing\n"),
1827 output_file);
1828 }
1829 sf_set_string(sf, SF_STR_SOFTWARE, "created by ATSA");
1830 /* allocate memory */
1831 obuf = (mus_sample_t **) csound->Malloc(csound, 2 * sizeof(mus_sample_t *));
1832 obuf[0] =
1833 (mus_sample_t *) csound->Calloc(csound,
1834 frm_samps * sizeof(mus_sample_t));
1835 obuf[1] =
1836 (mus_sample_t *) csound->Calloc(csound,
1837 frm_samps * sizeof(mus_sample_t));
1838 /* compute residual frame by frame */
1839 for (frm = 1; frm < frames; frm++) {
1840 /* clean buffers up */
1841 memset(in_buff, '\0', frm_samps * sizeof(double));
1842 memset(synth_buff, '\0', frm_samps * sizeof(double));
1843 /* for (i = 0; i < frm_samps; i++) */
1844 /* in_buff[i] = synth_buff[i] = 0.0; */
1845 frm_1 = frm - 1;
1846 frm_2 = frm;
1847 /* read frame from input */
1848 read_frame(fil, fil_len, win_samps[frm_1], win_samps[frm_2], in_buff);
1849 /* compute one synthesis frame */
1850 for (par = 0; par < partials; par++) {
1851 a1 = sound->amp[par][frm_1];
1852 a2 = sound->amp[par][frm_2];
1853 /* have to convert the frequency into radians per sample!!! */
1854 f1 = sound->frq[par][frm_1] * mag;
1855 f2 = sound->frq[par][frm_2] * mag;
1856 /* f1 *= mag; */
1857 /* f2 *= mag; */
1858 p1 = sound->pha[par][frm_1];
1859 p2 = sound->pha[par][frm_2];
1860 if (!(a1 <= 0.0 && a2 <= 0.0)) {
1861 /* check amp 0 in frame 1 */
1862 if (a1 <= 0.0) {
1863 f1 = f2;
1864 p1 = p2 - (f2 * frm_samps);
1865 while (p1 > PI)
1866 p1 -= TWOPI;
1867 while (p1 < (PI * -1))
1868 p1 += TWOPI;
1869 }
1870 /* check amp 0 in frame 2 */
1871 if (a2 <= 0.0) {
1872 f2 = f1;
1873 p2 = p1 + (f1 * frm_samps);
1874 while (p2 > PI)
1875 p2 -= TWOPI;
1876 while (p2 < (PI * -1))
1877 p2 += TWOPI;
1878 }
1879 synth_buffer(a1, a2, f1, f2, p1, p2, synth_buff, frm_samps);
1880 }
1881 }
1882 /* write output: chan 0 residual chan 1 synthesis */
1883 for (i = 0; i < frm_samps; i++) {
1884 synth = synth_buff[i];
1885 diff = in_buff[i] - synth;
1886 obuf[0][i] = (mus_sample_t) diff;
1887 obuf[1][i] = (mus_sample_t) synth;
1888 out_smp++;
1889 }
1890 atsa_sound_write_noninterleaved(sf, obuf, 2, frm_samps);
1891 }
1892 csound->Free(csound, in_buff);
1893 csound->Free(csound, synth_buff);
1894 /* update header and close output file */
1895 csound->FileClose(csound, fd);
1896 csound->Free(csound, obuf[0]);
1897 csound->Free(csound, obuf[1]);
1898 csound->Free(csound, obuf);
1899 }
1900
1901 /* ------------------------------------------------------------------------ */
1902
1903 /* ats_save
1904 * ========
1905 * saves an ATS_SOUND to disk.
1906 * sound: pointer to ATS_SOUND structure
1907 * outfile: pointer to output ats file
1908 * SMR_thres: partials with and avreage SMR
1909 * below this value are considered masked
1910 * and not written out to the ats file
1911 * type: file type
1912 * NOTE: sound MUST be optimised using optimize_sound
1913 * before calling this function
1914 */
ats_save(CSOUND * csound,ATS_SOUND * sound,FILE * outfile,float SMR_thres,int type)1915 static void ats_save(CSOUND *csound, ATS_SOUND *sound, FILE *outfile,
1916 float SMR_thres, int type)
1917 {
1918 int frm, i, par, dead = 0;
1919 double daux;
1920 ATS_HEADER header;
1921
1922 if (UNLIKELY(sound->optimized == NIL)) {
1923 csound->Die(csound, "%s", Str("Error: sound not optimised !"));
1924 }
1925 /* count how many partials are dead
1926 * unfortunately we have to do this first to
1927 * write the number of partials in the header
1928 */
1929 for (i = 0; i < sound->partials; i++) {
1930 /* see if partial is dead */
1931 if (!(sound->av[i].frq > 0.0) || !(sound->av[i].smr >= SMR_thres)) {
1932 dead++;
1933 }
1934 }
1935 /* sort partials by increasing frequency */
1936 qsort(sound->av, sound->partials, sizeof(ATS_PEAK), peak_frq_inc);
1937 /* fill header up */
1938 header.mag = 123.0;
1939 header.sr = (double) sound->srate;
1940 header.fs = (double) sound->frame_size;
1941 header.ws = (double) sound->window_size;
1942 header.par = (double) (sound->partials - dead);
1943 header.fra = (double) sound->frames;
1944 header.ma = sound->ampmax;
1945 header.mf = sound->frqmax;
1946 header.dur = sound->dur;
1947 header.typ = (double) type;
1948 /* write header */
1949 fseek(outfile, 0, SEEK_SET);
1950 if (UNLIKELY(1!=fwrite(&header, sizeof(ATS_HEADER), 1, outfile)))
1951 fprintf(stderr, "%s", Str("Write failure\n"));
1952 /* write frame data */
1953 for (frm = 0; frm < sound->frames; frm++) {
1954 daux = sound->time[0][frm];
1955 if (UNLIKELY(1!=fwrite(&daux, sizeof(double), 1, outfile)))
1956 fprintf(stderr, "%s", Str("Write failure\n"));
1957 for (i = 0; i < sound->partials; i++) {
1958 /* we ouput data in increasing frequency order
1959 * and we check for dead partials
1960 */
1961 if ((sound->av[i].frq > 0.0) && (sound->av[i].smr >= SMR_thres)) {
1962 /* get partial number from sound */
1963 par = sound->av[i].track;
1964 /* output data to file */
1965 daux = sound->amp[par][frm];
1966 if (UNLIKELY(1!=fwrite(&daux, sizeof(double), 1, outfile)))
1967 fprintf(stderr, "%s", Str("Write failure\n"));
1968 daux = sound->frq[par][frm];
1969 if (UNLIKELY(1!=fwrite(&daux, sizeof(double), 1, outfile)))
1970 fprintf(stderr, "%s", Str("Write failure\n"));
1971 if (type == 2 || type == 4) {
1972 daux = sound->pha[par][frm];
1973 if (UNLIKELY(1!=fwrite(&daux, sizeof(double), 1, outfile)))
1974 fprintf(stderr, "%s", Str("Write failure\n"));
1975 }
1976 }
1977 }
1978 /* write noise data */
1979 if (type == 3 || type == 4) {
1980 for (i = 0; i < ATSA_CRITICAL_BANDS; i++) {
1981 daux = sound->band_energy[i][frm];
1982 if (UNLIKELY(1!=fwrite(&daux, sizeof(double), 1, outfile)))
1983 fprintf(stderr, "%s", Str("Write failure\n"));
1984 }
1985 }
1986 }
1987 }
1988
1989 /* ------------------------------------------------------------------------ */
1990
1991 /* private function prototypes */
1992 static int compute_frames(ANARGS *anargs);
1993
1994 /* ATS_SOUND *tracker (ANARGS *anargs, char *soundfile)
1995 * partial tracking function
1996 * anargs: pointer to analysis parameters
1997 * soundfile: path to input file
1998 * returns an ATS_SOUND with data issued from analysis
1999 */
tracker(CSOUND * csound,ANARGS * anargs,char * soundfile,char * resfile)2000 static ATS_SOUND *tracker(CSOUND *csound, ANARGS *anargs, char *soundfile,
2001 char *resfile)
2002 {
2003 int M_2, first_point, filptr, n_partials = 0;
2004 int frame_n, k, sflen, *win_samps, peaks_size, tracks_size = 0;
2005 int i, frame, i_tmp;
2006 float *window, norm, sfdur, f_tmp;
2007
2008 /* declare structures and buffers */
2009 ATS_SOUND *sound = NULL;
2010 ATS_PEAK *peaks, *tracks = NULL, cpy_peak;
2011 ATS_FRAME *ana_frames = NULL, *unmatched_peaks = NULL;
2012 mus_sample_t **bufs;
2013 ATS_FFT fft;
2014 SF_INFO sfinfo;
2015 SNDFILE *sf;
2016 void *fd;
2017
2018 /* open input file
2019 we get srate and total_samps in file in anargs */
2020 memset(&sfinfo, 0, sizeof(SF_INFO));
2021 fd = csound->FileOpen2(csound, &sf, CSFILE_SND_R, soundfile, &sfinfo,
2022 "SFDIR;SSDIR", CSFTYPE_UNKNOWN_AUDIO, 0);
2023 if (UNLIKELY(fd == NULL)) {
2024 csound->ErrorMsg(csound, Str("atsa: cannot open input file '%s': %s"),
2025 soundfile, Str(sf_strerror(NULL)));
2026 return NULL;
2027 }
2028 /* warn about multi-channel sound files */
2029 if (UNLIKELY(sfinfo.channels != 1)) {
2030 csound->ErrorMsg(csound,
2031 Str("atsa: file has %d channels, must be mono !"),
2032 (int) sfinfo.channels);
2033 return NULL;
2034 }
2035
2036 csound->Message(csound, "%s", Str("tracking...\n"));
2037
2038 /* get sample rate and # of frames from file header */
2039 anargs->srate = sfinfo.samplerate;
2040 sflen = (int) sfinfo.frames;
2041 sfdur = (float) sflen / anargs->srate;
2042 /* check analysis parameters */
2043 /* check start time */
2044 if (UNLIKELY(!(anargs->start >= 0.0 && anargs->start < sfdur))) {
2045 csound->Warning(csound, Str("start %f out of bounds, corrected to 0.0"),
2046 anargs->start);
2047 anargs->start = 0.0f;
2048 }
2049 /* check duration */
2050 if (anargs->duration == ATSA_DUR) {
2051 anargs->duration = sfdur - anargs->start;
2052 }
2053 f_tmp = anargs->duration + anargs->start;
2054 if (UNLIKELY(!(anargs->duration > 0.0 && f_tmp <= sfdur))) {
2055 csound->Warning(csound, Str("duration %f out of bounds, "
2056 "limited to file duration"),
2057 anargs->duration);
2058 anargs->duration = sfdur - anargs->start;
2059 }
2060 /* print time bounds */
2061 csound->Message(csound, Str("start: %f duration: %f file dur: %f\n"),
2062 anargs->start, anargs->duration, sfdur);
2063 /* check lowest frequency */
2064 if (UNLIKELY(!
2065 (anargs->lowest_freq > 0.0 &&
2066 anargs->lowest_freq < anargs->highest_freq))) {
2067 csound->Warning(csound,
2068 Str("lowest freq. %f out of bounds, "
2069 "forced to default: %f"), anargs->lowest_freq,
2070 ATSA_LFREQ);
2071 anargs->lowest_freq = ATSA_LFREQ;
2072 }
2073 /* check highest frequency */
2074 if (UNLIKELY(!
2075 (anargs->highest_freq > anargs->lowest_freq &&
2076 anargs->highest_freq <= anargs->srate * 0.5))) {
2077 csound->Warning(csound,
2078 Str("highest freq. %f out of bounds, "
2079 "forced to default: %f"), anargs->highest_freq,
2080 ATSA_HFREQ);
2081 anargs->highest_freq = ATSA_HFREQ;
2082 }
2083 /* frequency deviation */
2084 if (UNLIKELY(!(anargs->freq_dev > 0.0f && anargs->freq_dev < 1.0f))) {
2085 csound->Warning(csound, Str("freq. dev. %f out of bounds, "
2086 "should be > 0.0 and <= 1.0, "
2087 "forced to default: %f"),
2088 anargs->freq_dev, ATSA_FREQDEV);
2089 anargs->freq_dev = ATSA_FREQDEV;
2090 }
2091 /* window cycles */
2092 if (UNLIKELY(!(anargs->win_cycles >= 1 && anargs->win_cycles <= 8))) {
2093 csound->Warning(csound, Str("windows cycles %d out of bounds, "
2094 "should be between 1 and 8, "
2095 "forced to default: %d"),
2096 anargs->win_cycles, ATSA_WCYCLES);
2097 anargs->win_cycles = ATSA_WCYCLES;
2098 }
2099 /* window type */
2100 if (UNLIKELY(!(anargs->win_type >= 0 && anargs->win_type <= 3))) {
2101 csound->Warning(csound, Str("window type %d out of bounds, "
2102 "should be between 0 and 3, "
2103 "forced to default: %d"),
2104 anargs->win_type, ATSA_WTYPE);
2105 anargs->win_type = ATSA_WTYPE;
2106 }
2107 /* hop size */
2108 if (UNLIKELY(!(anargs->hop_size > 0.0 && anargs->hop_size <= 1.0))) {
2109 csound->Warning(csound, Str("hop size %f out of bounds, "
2110 "should be > 0.0 and <= 1.0, "
2111 "forced to default: %f"),
2112 anargs->hop_size, ATSA_HSIZE);
2113 anargs->hop_size = ATSA_HSIZE;
2114 }
2115 /* lowest mag */
2116 if (UNLIKELY(!(anargs->lowest_mag <= 0.0))) {
2117 csound->Warning(csound, Str("lowest magnitude %f out of bounds, "
2118 "should be >= 0.0 and <= 1.0, "
2119 "forced to default: %f"),
2120 anargs->lowest_mag, ATSA_LMAG);
2121 anargs->lowest_mag = ATSA_LMAG;
2122 }
2123 /* set some values before checking next set of parameters */
2124 anargs->first_smp = (int) floor(anargs->start * (float) anargs->srate);
2125 anargs->total_samps = (int) floor(anargs->duration * (float) anargs->srate);
2126 /* fundamental cycles */
2127 anargs->cycle_smp =
2128 (int) floor((double) anargs->win_cycles * (double) anargs->srate /
2129 (double) anargs->lowest_freq);
2130 /* window size */
2131 anargs->win_size =
2132 (anargs->cycle_smp % 2 ==
2133 0) ? anargs->cycle_smp + 1 : anargs->cycle_smp;
2134 /* calculate hop samples */
2135 anargs->hop_smp = (int)floor((float) anargs->win_size * anargs->hop_size);
2136 /* compute total number of frames */
2137 anargs->frames = compute_frames(anargs);
2138 /* check that we have enough frames for the analysis */
2139 if (UNLIKELY(!(anargs->frames >= ATSA_MFRAMES))) {
2140 csound->ErrorMsg(csound,
2141 Str("atsa: %d frames are not enough for analysis, "
2142 "need at least %d"), anargs->frames, ATSA_MFRAMES);
2143 return NULL;
2144 }
2145 /* check other user parameters */
2146 /* track length */
2147 if (UNLIKELY(!(anargs->track_len >= 1 && anargs->track_len < anargs->frames))) {
2148 i_tmp = (ATSA_TRKLEN < anargs->frames) ? ATSA_TRKLEN : anargs->frames - 1;
2149 csound->Warning(csound,
2150 Str("track length %d out of bounds, forced to: %d"),
2151 anargs->track_len, i_tmp);
2152 anargs->track_len = i_tmp;
2153 }
2154 /* min. segment length */
2155 if (UNLIKELY(!(anargs->min_seg_len >= 1 &&
2156 anargs->min_seg_len < anargs->frames))) {
2157 i_tmp =
2158 (ATSA_MSEGLEN < anargs->frames) ? ATSA_MSEGLEN : anargs->frames - 1;
2159 csound->Warning(csound,
2160 Str("min. segment length %d out of bounds, "
2161 "forced to: %d"), anargs->min_seg_len, i_tmp);
2162 anargs->min_seg_len = i_tmp;
2163 }
2164 /* min. gap length */
2165 if (UNLIKELY(!(anargs->min_gap_len >= 0 &&
2166 anargs->min_gap_len < anargs->frames))) {
2167 i_tmp =
2168 (ATSA_MGAPLEN < anargs->frames) ? ATSA_MGAPLEN : anargs->frames - 1;
2169 csound->Warning(csound,
2170 Str("min. gap length %d out of bounds, forced to: %d"),
2171 anargs->min_gap_len, i_tmp);
2172 anargs->min_gap_len = i_tmp;
2173 }
2174 /* SMR threshold */
2175 if (UNLIKELY(!(anargs->SMR_thres >= 0.0 &&
2176 anargs->SMR_thres < ATSA_MAX_DB_SPL))) {
2177 csound->Warning(csound, Str("SMR threshold %f out of bounds, "
2178 "should be >= 0.0 and < %f dB SPL, "
2179 "forced to default: %f"),
2180 anargs->SMR_thres, ATSA_MAX_DB_SPL, ATSA_SMRTHRES);
2181 anargs->SMR_thres = ATSA_SMRTHRES;
2182 }
2183 /* min. seg. SMR */
2184 if (UNLIKELY(!
2185 (anargs->min_seg_SMR >= anargs->SMR_thres &&
2186 anargs->min_seg_SMR < ATSA_MAX_DB_SPL))) {
2187 csound->Warning(csound,
2188 Str("min. seg. SMR %f out of bounds, "
2189 "should be >= %f and < %f dB SPL, "
2190 "forced to default: %f"), anargs->min_seg_SMR,
2191 anargs->SMR_thres, ATSA_MAX_DB_SPL, ATSA_MSEGSMR);
2192 anargs->min_seg_SMR = ATSA_MSEGSMR;
2193 }
2194 /* last peak contribution */
2195 if (UNLIKELY(!(anargs->last_peak_cont >= 0.0 &&
2196 anargs->last_peak_cont <= 1.0))) {
2197 csound->Warning(csound, Str("last peak contribution %f out of bounds, "
2198 "should be >= 0.0 and <= 1.0, "
2199 "forced to default: %f"),
2200 anargs->last_peak_cont, ATSA_LPKCONT);
2201 anargs->last_peak_cont = ATSA_LPKCONT;
2202 }
2203 /* SMR cont. */
2204 if (UNLIKELY(!(anargs->SMR_cont >= 0.0 && anargs->SMR_cont <= 1.0))) {
2205 csound->Warning(csound, Str("SMR contribution %f out of bounds, "
2206 "should be >= 0.0 and <= 1.0, "
2207 "forced to default: %f"),
2208 anargs->SMR_cont, ATSA_SMRCONT);
2209 anargs->SMR_cont = ATSA_SMRCONT;
2210 }
2211 /* continue computing parameters */
2212 /* fft size */
2213 anargs->fft_size = ppp2(2 * anargs->win_size);
2214
2215 /* allocate memory for sound, we read the whole sound in memory */
2216 bufs = (mus_sample_t **) csound->Malloc(csound, sizeof(mus_sample_t *));
2217 bufs[0] =
2218 (mus_sample_t *) csound->Malloc(csound, sflen * sizeof(mus_sample_t));
2219 /* make our window */
2220 window = make_window(csound, anargs->win_type, anargs->win_size);
2221 /* get window norm */
2222 norm = window_norm(window, anargs->win_size);
2223 /* fft mag for computing frequencies */
2224 anargs->fft_mag = (double) anargs->srate / (double) anargs->fft_size;
2225 /* lowest fft bin for analysis */
2226 anargs->lowest_bin = (int)floor(anargs->lowest_freq / anargs->fft_mag);
2227 /* highest fft bin for analisis */
2228 anargs->highest_bin = (int)floor(anargs->highest_freq / anargs->fft_mag);
2229 /* allocate an array analysis frames in memory */
2230 ana_frames =
2231 (ATS_FRAME *) csound->Malloc(csound,
2232 anargs->frames * sizeof(ATS_FRAME));
2233 /* alocate memory to store mid-point window sample numbers */
2234 win_samps = (int *) csound->Malloc(csound, anargs->frames * sizeof(int));
2235 /* center point of window */
2236 M_2 = (anargs->win_size-1)/2; /* Was (int)floor((anargs->win_size - 1) / 2) */
2237 /* first point in fft buffer to write */
2238 first_point = anargs->fft_size - M_2;
2239 /* half a window from first sample */
2240 filptr = anargs->first_smp - M_2;
2241 /* read sound into memory */
2242 atsa_sound_read_noninterleaved(sf, bufs, 1, sflen);
2243
2244 /* make our fft-struct */
2245 fft.size = anargs->fft_size;
2246 fft.rate = anargs->srate;
2247 fft.data =
2248 (MYFLT *) csound->Malloc(csound,
2249 (anargs->fft_size + 2) * sizeof(MYFLT));
2250
2251 /* main loop */
2252 for (frame_n = 0; frame_n < anargs->frames; frame_n++) {
2253 /* clear fft arrays */
2254 for (k = 0; k < (fft.size + 2); k++)
2255 fft.data[k] = (MYFLT) 0;
2256 /* multiply by window */
2257 for (k = 0; k < anargs->win_size; k++) {
2258 if ((filptr >= 0) && (filptr < sflen))
2259 fft.data[(k + first_point) % anargs->fft_size] =
2260 (MYFLT) window[k] * (MYFLT) bufs[0][filptr];
2261 filptr++;
2262 }
2263 /* we keep sample numbers of window midpoints in win_samps array */
2264 win_samps[frame_n] = filptr - M_2 - 1;
2265 /* move file pointer back */
2266 filptr = filptr - anargs->win_size + anargs->hop_smp;
2267 /* take the fft */
2268 csound->RealFFTnp2(csound, fft.data, fft.size);
2269 /* peak detection */
2270 peaks_size = 0;
2271 peaks =
2272 peak_detection(csound, &fft, anargs->lowest_bin, anargs->highest_bin,
2273 anargs->lowest_mag, norm, &peaks_size);
2274 /* peak tracking */
2275 if (peaks != NULL) {
2276 /* evaluate peaks SMR (masking curves) */
2277 evaluate_smr(peaks, peaks_size);
2278 if (frame_n) {
2279 /* initialise or update tracks */
2280 if ((tracks =
2281 update_tracks(csound, tracks, &tracks_size, anargs->track_len,
2282 frame_n, ana_frames,
2283 anargs->last_peak_cont)) != NULL) {
2284 /* do peak matching */
2285 unmatched_peaks =
2286 peak_tracking(csound, tracks, &tracks_size, peaks, &peaks_size,
2287 anargs->freq_dev, 2.0 * anargs->SMR_cont,
2288 &n_partials);
2289 /* kill unmatched peaks from previous frame */
2290 if (unmatched_peaks[0].peaks != NULL) {
2291 for (k = 0; k < unmatched_peaks[0].n_peaks; k++) {
2292 cpy_peak = unmatched_peaks[0].peaks[k];
2293 cpy_peak.amp = cpy_peak.smr = 0.0;
2294 peaks = push_peak(csound, &cpy_peak, peaks, &peaks_size);
2295 }
2296 csound->Free(csound, unmatched_peaks[0].peaks);
2297 }
2298 /* give birth to peaks from new frame */
2299 if (unmatched_peaks[1].peaks != NULL) {
2300 for (k = 0; k < unmatched_peaks[1].n_peaks; k++) {
2301 tracks =
2302 push_peak(csound, &unmatched_peaks[1].peaks[k], tracks,
2303 &tracks_size);
2304 unmatched_peaks[1].peaks[k].amp =
2305 unmatched_peaks[1].peaks[k].smr = 0.0;
2306 ana_frames[frame_n - 1].peaks =
2307 push_peak(csound, &unmatched_peaks[1].peaks[k],
2308 ana_frames[frame_n - 1].peaks,
2309 &ana_frames[frame_n - 1].n_peaks);
2310 }
2311 csound->Free(csound, unmatched_peaks[1].peaks);
2312 }
2313 }
2314 else {
2315 /* give number to all peaks */
2316 qsort(peaks, peaks_size, sizeof(ATS_PEAK), peak_frq_inc);
2317 for (k = 0; k < peaks_size; k++)
2318 peaks[k].track = n_partials++;
2319 }
2320 }
2321 else {
2322 /* give number to all peaks */
2323 qsort(peaks, peaks_size, sizeof(ATS_PEAK), peak_frq_inc);
2324 for (k = 0; k < peaks_size; k++)
2325 peaks[k].track = n_partials++;
2326 }
2327 /* attach peaks to ana_frames */
2328 ana_frames[frame_n].peaks = peaks;
2329 ana_frames[frame_n].n_peaks = n_partials;
2330 ana_frames[frame_n].time =
2331 (double) (win_samps[frame_n] -
2332 anargs->first_smp) / (double) anargs->srate;
2333 /* free memory */
2334 csound->Free(csound, unmatched_peaks);
2335 }
2336 else {
2337 /* if no peaks found, initialise empty frame */
2338 ana_frames[frame_n].peaks = NULL;
2339 ana_frames[frame_n].n_peaks = 0;
2340 ana_frames[frame_n].time =
2341 (double) (win_samps[frame_n] -
2342 anargs->first_smp) / (double) anargs->srate;
2343 }
2344 }
2345 /* free up some memory */
2346 csound->Free(csound, window);
2347 csound->Free(csound, tracks);
2348 csound->Free(csound, fft.data);
2349 /* init sound */
2350 csound->Message(csound, "%s", Str("Initializing ATS data..."));
2351 sound = (ATS_SOUND *) csound->Malloc(csound, sizeof(ATS_SOUND));
2352 init_sound(csound, sound, anargs->srate,
2353 (int) (anargs->hop_size * anargs->win_size), anargs->win_size,
2354 anargs->frames, anargs->duration, n_partials,
2355 ((anargs->type == 3 || anargs->type == 4) ? 1 : 0));
2356 /* store values from frames into the arrays */
2357 for (k = 0; k < n_partials; k++) {
2358 for (frame = 0; frame < sound->frames; frame++) {
2359 sound->time[k][frame] = ana_frames[frame].time;
2360 for (i = 0; i < ana_frames[frame].n_peaks; i++)
2361 if (ana_frames[frame].peaks[i].track == k) {
2362 sound->amp[k][frame] = ana_frames[frame].peaks[i].amp;
2363 sound->frq[k][frame] = ana_frames[frame].peaks[i].frq;
2364 sound->pha[k][frame] = ana_frames[frame].peaks[i].pha;
2365 sound->smr[k][frame] = ana_frames[frame].peaks[i].smr;
2366 }
2367 }
2368 }
2369 csound->Message(csound, "%s", Str("done!\n"));
2370 /* free up ana_frames memory */
2371 /* first, free all peaks in each slot of ana_frames... */
2372 for (k = 0; k < anargs->frames; k++)
2373 csound->Free(csound, ana_frames[k].peaks);
2374 /* ...then free ana_frames */
2375 csound->Free(csound, ana_frames);
2376 /* optimise sound */
2377 optimize_sound(csound, anargs, sound);
2378 /* compute residual */
2379 if (UNLIKELY(anargs->type == 3 || anargs->type == 4)) {
2380 csound->Message(csound, "%s", Str("Computing residual..."));
2381 compute_residual(csound, bufs, sflen, resfile, sound, win_samps,
2382 anargs->srate);
2383 csound->Message(csound, "%s", Str("done!\n"));
2384 }
2385 /* free the rest of the memory */
2386 csound->Free(csound, win_samps);
2387 csound->Free(csound, bufs[0]);
2388 csound->Free(csound, bufs);
2389 /* analyse residual */
2390 if (UNLIKELY(anargs->type == 3 || anargs->type == 4)) {
2391 #ifdef WIN32
2392 char buffer[160];
2393 char * tmp = getenv("TEMP");
2394 strNcpy(buffer, tmp, 160);
2395 // MKG 2014 Jan 29: No linkage for strlcat with MinGW here.
2396 // snd corrected
2397 //strlcat(buffer, ATSA_RES_FILE, 160);
2398 strncat(buffer, ATSA_RES_FILE, 159-strlen(buffer)); buffer[159]='\0';
2399 csound->Message(csound, "%s", Str("Analysing residual..."));
2400 residual_analysis(csound, buffer, sound);
2401 #else
2402 csound->Message(csound, "%s", Str("Analysing residual..."));
2403 residual_analysis(csound, ATSA_RES_FILE, sound);
2404 #endif
2405 csound->Message(csound, "%s", Str("done!\n"));
2406 }
2407 csound->Message(csound, "%s", Str("tracking completed.\n"));
2408 return (sound);
2409 }
2410
2411 /* int compute_frames(ANARGS *anargs)
2412 * computes number of analysis frames from the user's parameters
2413 * returns the number of frames
2414 * anargs: pointer to analysis parameters
2415 */
compute_frames(ANARGS * anargs)2416 static int compute_frames(ANARGS *anargs)
2417 {
2418 int n_frames =
2419 (int) floor((float) anargs->total_samps / (float) anargs->hop_smp);
2420
2421 while ((n_frames++ * anargs->hop_smp - anargs->hop_smp +
2422 anargs->first_smp) < (anargs->first_smp + anargs->total_samps));
2423 return (n_frames);
2424 }
2425
2426 /* ------------------------------------------------------------------------ */
2427
2428 /* private function prototypes */
2429 static int find_next_val_arr(double *arr, int beg, int size);
2430 static int find_next_zero_arr(double *arr, int beg, int size);
2431 static int find_prev_val_arr(double *arr, int beg);
2432 static void fill_sound_gaps(CSOUND *csound, ATS_SOUND *sound, int min_gap_len);
2433 static void trim_partials(CSOUND *csound, ATS_SOUND *sound, int min_seg_len,
2434 float min_seg_smr);
2435 static void set_av(CSOUND *csound, ATS_SOUND *sound);
2436
2437 /* various conversion functions
2438 * to deal with dB and dB SPL
2439 * they take and return double floats
2440 */
amp2db(double amp)2441 static inline double amp2db(double amp)
2442 {
2443 return (20.0 * log10(amp));
2444 }
2445
db2amp(double db)2446 static inline double db2amp(double db)
2447 {
2448 return (pow(10.0, db / 20.0));
2449 }
2450
amp2db_spl(double amp)2451 static inline double amp2db_spl(double amp)
2452 {
2453 return (amp2db(amp) + ATSA_MAX_DB_SPL);
2454 }
2455
2456
2457 /*
2458 static inline double db2amp_spl(double db_spl)
2459 {
2460 return (db2amp(db_spl - ATSA_MAX_DB_SPL));
2461 }
2462 */
2463
2464 /* ppp2
2465 * ====
2466 * returns the closest power of two
2467 * greater than num
2468 */
ppp2(int num)2469 static inline unsigned int ppp2(int num)
2470 {
2471 unsigned int tmp = 2;
2472
2473 while (tmp < (unsigned int) num)
2474 tmp = tmp << 1;
2475 return (tmp);
2476 }
2477
2478 /* optimize_sound
2479 * ==============
2480 * optimises an ATS_SOUND in memory before saving
2481 * anargs: pointer to analysis parameters
2482 * sound: pointer to ATS_SOUND structure
2483 */
optimize_sound(CSOUND * csound,ANARGS * anargs,ATS_SOUND * sound)2484 static void optimize_sound(CSOUND *csound, ANARGS *anargs, ATS_SOUND *sound)
2485 {
2486 double ampmax = 0.0, frqmax = 0.0;
2487 int frame, partial;
2488
2489 for (frame = 0; frame < sound->frames; frame++)
2490 for (partial = 0; partial < sound->partials; partial++) {
2491 if (ampmax < sound->amp[partial][frame])
2492 ampmax = sound->amp[partial][frame];
2493 if (frqmax < sound->frq[partial][frame])
2494 frqmax = sound->frq[partial][frame];
2495 }
2496 sound->ampmax = ampmax;
2497 sound->frqmax = frqmax;
2498
2499 fill_sound_gaps(csound, sound, anargs->min_gap_len);
2500 trim_partials(csound, sound, anargs->min_seg_len, anargs->min_seg_SMR);
2501 set_av(csound, sound);
2502 /* finally set slot to 1 */
2503 sound->optimized = 1;
2504 }
2505
2506 /* fill_sound_gaps
2507 * ===============
2508 * fills gaps in ATS_SOUND partials by interpolation
2509 * sound: pointer to ATS_SOUND
2510 * min_gap_len: minimum gap length, gaps shorter or equal to this
2511 * value will be filled in by interpolation
2512 */
fill_sound_gaps(CSOUND * csound,ATS_SOUND * sound,int min_gap_len)2513 static void fill_sound_gaps(CSOUND *csound, ATS_SOUND *sound, int min_gap_len)
2514 {
2515 int i, j, k, next_val, next_zero, prev_val, gap_size;
2516 double f_inc, a_inc, s_inc, mag = TWOPI / (double) sound->srate;
2517
2518 csound->Message(csound, "%s", Str("Filling sound gaps..."));
2519 for (i = 0; i < sound->partials; i++) {
2520 /* first we fix the freq gap before attack */
2521 next_val = find_next_val_arr(sound->frq[i], 0, sound->frames);
2522 if (next_val > 0) {
2523 for (j = 0; j < next_val; j++) {
2524 sound->frq[i][j] = sound->frq[i][next_val];
2525 }
2526 }
2527 /* fix the freq gap at end */
2528 prev_val = find_prev_val_arr(sound->frq[i], sound->frames - 1);
2529 if (prev_val != NIL && prev_val < sound->frames - 1) {
2530 for (j = prev_val; j < sound->frames; j++) {
2531 sound->frq[i][j] = sound->frq[i][prev_val];
2532 }
2533 }
2534 /* now we fix inner gaps of frq, pha, and amp */
2535 k = find_next_val_arr(sound->amp[i], 0, sound->frames);
2536 while (k < sound->frames && k != NIL) {
2537 /* find next gap: we consider gaps in amplitude, */
2538 /* we fix frequency and phase in parallel */
2539 next_zero = find_next_zero_arr(sound->amp[i], k, sound->frames);
2540 if (next_zero != NIL) {
2541 prev_val = next_zero - 1;
2542 next_val = find_next_val_arr(sound->amp[i], next_zero, sound->frames);
2543 /* check if we didn't get to end of array */
2544 if (next_val == NIL)
2545 break;
2546 gap_size = next_val - prev_val;
2547 /* csound->Message(csound,
2548 "par: %d prev_val: %d next_val: %d gap_size %d\n",
2549 i, prev_val, next_val, gap_size); */
2550 /* check validity of found gap */
2551 if (gap_size <= min_gap_len) {
2552 /* csound->Message(csound, "Filling gap of par: %d\n", i); */
2553 f_inc =
2554 (sound->frq[i][next_val] - sound->frq[i][prev_val]) / gap_size;
2555 a_inc =
2556 (sound->amp[i][next_val] - sound->amp[i][prev_val]) / gap_size;
2557 s_inc =
2558 (sound->smr[i][next_val] - sound->smr[i][prev_val]) / gap_size;
2559 for (j = next_zero; j < next_val; j++) {
2560 sound->frq[i][j] = sound->frq[i][j - 1] + f_inc;
2561 sound->amp[i][j] = sound->amp[i][j - 1] + a_inc;
2562 sound->smr[i][j] = sound->smr[i][j - 1] + s_inc;
2563 sound->pha[i][j] =
2564 sound->pha[i][j - 1] -
2565 (sound->frq[i][j] * sound->frame_size * mag);
2566 /* wrap phase */
2567 while (sound->pha[i][j] > PI) {
2568 sound->pha[i][j] -= TWOPI;
2569 }
2570 while (sound->pha[i][j] < (PI * (-1.0))) {
2571 sound->pha[i][j] += TWOPI;
2572 }
2573 }
2574 /* gap fixed, find next gap */
2575 k = next_val;
2576 }
2577 else {
2578 /* gap not valid, move to next one */
2579 /* csound->Message(csound, "jumping to next_val: %d\n", next_val); */
2580 k = next_val;
2581 }
2582 }
2583 else {
2584 break;
2585 }
2586 }
2587 }
2588 csound->Message(csound, "%s", Str("done!\n"));
2589 }
2590
2591 /* trim_partials
2592 * =============
2593 * trim short segments from ATS_SOUND partials
2594 * sound: pointer to ATS_SOUND
2595 * min_seg_len: minimum segment length, segments shorter or equal
2596 * to this value will be candidates for trimming
2597 * min_seg_smr: minimum segment average SMR, segment candidates
2598 * should have an average SMR below this value to be trimmed
2599 */
trim_partials(CSOUND * csound,ATS_SOUND * sound,int min_seg_len,float min_seg_smr)2600 static void trim_partials(CSOUND *csound, ATS_SOUND *sound, int min_seg_len,
2601 float min_seg_smr)
2602 {
2603 int i, j, k, seg_beg, seg_end, seg_size, count = 0;
2604 double val = 0.0, smr_av = 0.0;
2605
2606 csound->Message(csound, "%s", Str("Trimming short partials..."));
2607 for (i = 0; i < sound->partials; i++) {
2608 k = 0;
2609 while (k < sound->frames) {
2610 /* find next segment */
2611 seg_beg = find_next_val_arr(sound->amp[i], k, sound->frames);
2612 if (seg_beg != NIL) {
2613 seg_end = find_next_zero_arr(sound->amp[i], seg_beg, sound->frames);
2614 /* check if we didn't get to end of array */
2615 if (seg_end == NIL)
2616 seg_end = sound->frames;
2617 seg_size = seg_end - seg_beg;
2618 /* csound->Message(csound,
2619 "par: %d seg_beg: %d seg_end: %d seg_size %d\n",
2620 i, seg_beg, seg_end, seg_size); */
2621 if (seg_size <= min_seg_len) {
2622 for (j = seg_beg; j < seg_end; j++) {
2623 if (sound->smr[i][j] > 0.0) {
2624 val += sound->smr[i][j];
2625 count++;
2626 }
2627 }
2628 if (count > 0)
2629 smr_av = val / (double) count;
2630 if (smr_av < min_seg_smr) {
2631 /* trim segment, only amplitude and SMR data */
2632 /* csound->Message(csound, "Trimming par: %d\n", i); */
2633 for (j = seg_beg; j < seg_end; j++) {
2634 sound->amp[i][j] = 0.0;
2635 sound->smr[i][j] = 0.0;
2636 }
2637 }
2638 k = seg_end;
2639 }
2640 else {
2641 /* segment not valid, move to the next one */
2642 /* csound->Message(csound, "jumping to seg_end: %d\n", seg_end); */
2643 k = seg_end;
2644 }
2645 }
2646 else {
2647 break;
2648 }
2649 }
2650 }
2651 csound->Message(csound, "%s", Str("done!\n"));
2652 }
2653
2654 /* auxiliary functions to fill_sound_gaps and trim_partials */
find_next_val_arr(double * arr,int beg,int size)2655 static int find_next_val_arr(double *arr, int beg, int size)
2656 {
2657 int j, next_val = NIL;
2658
2659 for (j = beg; j < size; j++)
2660 if (arr[j] > 0.0) {
2661 next_val = j;
2662 break;
2663 }
2664 return (next_val);
2665 }
2666
find_next_zero_arr(double * arr,int beg,int size)2667 static int find_next_zero_arr(double *arr, int beg, int size)
2668 {
2669 int j, next_zero = NIL;
2670
2671 for (j = beg; j < size; j++)
2672 if (arr[j] == 0.0) {
2673 next_zero = j;
2674 break;
2675 }
2676 return (next_zero);
2677 }
2678
find_prev_val_arr(double * arr,int beg)2679 static int find_prev_val_arr(double *arr, int beg)
2680 {
2681 int j, prev_val = NIL;
2682
2683 for (j = beg; j >= 0; j--)
2684 if (arr[j] != 0.0) {
2685 prev_val = j;
2686 break;
2687 }
2688 return (prev_val);
2689 }
2690
2691 /* set_av
2692 * ======
2693 * sets the av structure slot of an ATS_SOUND,
2694 * it computes the average freq. and SMR for each partial
2695 * sound: pointer to ATS_SOUND structure
2696 */
set_av(CSOUND * csound,ATS_SOUND * sound)2697 static void set_av(CSOUND *csound, ATS_SOUND *sound)
2698 {
2699 int i, j, count;
2700 double val;
2701
2702 csound->Message(csound, "%s", Str("Computing averages..."));
2703 for (i = 0; i < sound->partials; i++) {
2704 /* smr */
2705 val = 0.0;
2706 count = 0;
2707 for (j = 0; j < sound->frames; j++) {
2708 if (sound->smr[i][j] > 0.0) {
2709 val += sound->smr[i][j];
2710 count++;
2711 }
2712 }
2713 if (count > 0) {
2714 sound->av[i].smr = val / (double) count;
2715 }
2716 else {
2717 sound->av[i].smr = 0.0;
2718 }
2719 /* csound->Message(csound, "par: %d smr_av: %f\n",
2720 i, (float)sound->av[i].smr); */
2721 /* frq */
2722 val = 0.0;
2723 count = 0;
2724 for (j = 0; j < sound->frames; j++) {
2725 if (sound->frq[i][j] > 0.0) {
2726 val += sound->frq[i][j];
2727 count++;
2728 }
2729 }
2730 if (count > 0) {
2731 sound->av[i].frq = val / (double) count;
2732 }
2733 else {
2734 sound->av[i].frq = 0.0;
2735 }
2736 /* set track# */
2737 sound->av[i].track = i;
2738 }
2739 csound->Message(csound, "%s", Str("done!\n"));
2740 }
2741
2742 /* init_sound
2743 * ==========
2744 * initialises a new sound allocating memory
2745 */
init_sound(CSOUND * csound,ATS_SOUND * sound,int sampling_rate,int frame_size,int window_size,int frames,double duration,int partials,int use_noise)2746 static void init_sound(CSOUND *csound, ATS_SOUND *sound, int sampling_rate,
2747 int frame_size, int window_size, int frames,
2748 double duration, int partials, int use_noise)
2749 {
2750 int i /* , j*/;
2751
2752 if (UNLIKELY(partials==0)) {
2753 csound->Die(csound, "%s", Str("No partials to track -- stopping\n"));
2754 }
2755 sound->srate = sampling_rate;
2756 sound->frame_size = frame_size;
2757 sound->window_size = window_size;
2758 sound->frames = frames;
2759 sound->dur = duration;
2760 sound->partials = partials;
2761 sound->av =
2762 (ATS_PEAK *) csound->Malloc(csound, partials * sizeof(ATS_PEAK));
2763 sound->optimized = NIL;
2764 sound->time = (void *) csound->Malloc(csound, partials * sizeof(void *));
2765 sound->frq = (void *) csound->Malloc(csound, partials * sizeof(void *));
2766 sound->amp = (void *) csound->Malloc(csound, partials * sizeof(void *));
2767 sound->pha = (void *) csound->Malloc(csound, partials * sizeof(void *));
2768 sound->smr = (void *) csound->Malloc(csound, partials * sizeof(void *));
2769 sound->res = (void *) csound->Malloc(csound, partials * sizeof(void *));
2770 /* allocate memory for time, amp, frq, smr, and res data */
2771 for (i = 0; i < partials; i++) {
2772 sound->time[i] =
2773 (double *) csound->Malloc(csound, frames * sizeof(double));
2774 sound->amp[i] =
2775 (double *) csound->Calloc(csound, frames * sizeof(double));
2776 sound->frq[i] =
2777 (double *) csound->Calloc(csound, frames * sizeof(double));
2778 sound->pha[i] =
2779 (double *) csound->Calloc(csound, frames * sizeof(double));
2780 sound->smr[i] =
2781 (double *) csound->Calloc(csound, frames * sizeof(double));
2782 sound->res[i] =
2783 (double *) csound->Calloc(csound, frames * sizeof(double));
2784 }
2785 /* init all array values with 0.0 */
2786 /* for (i = 0; i < partials; i++) */
2787 /* for (j = 0; j < frames; j++) { */
2788 /* sound->amp[i][j] = 0.0; */
2789 /* sound->frq[i][j] = 0.0; */
2790 /* sound->pha[i][j] = 0.0; */
2791 /* sound->smr[i][j] = 0.0; */
2792 /* sound->res[i][j] = 0.0; */
2793 /* } */
2794 if (use_noise) {
2795 sound->band_energy =
2796 (double **) csound->Malloc(csound,
2797 ATSA_CRITICAL_BANDS * sizeof(double *));
2798 for (i = 0; i < ATSA_CRITICAL_BANDS; i++)
2799 sound->band_energy[i] =
2800 (double *) csound->Malloc(csound, frames * sizeof(double));
2801 }
2802 else
2803 sound->band_energy = NULL;
2804 }
2805
2806 /* free_sound
2807 * ==========
2808 * frees sound's memory
2809 */
free_sound(CSOUND * csound,ATS_SOUND * sound)2810 static void free_sound(CSOUND *csound, ATS_SOUND *sound)
2811 {
2812 int k;
2813
2814 if (sound != NULL) {
2815 csound->Free(csound, sound->av);
2816 /* data */
2817 for (k = 0; k < sound->partials; k++) {
2818 csound->Free(csound, sound->time[k]);
2819 csound->Free(csound, sound->amp[k]);
2820 csound->Free(csound, sound->frq[k]);
2821 csound->Free(csound, sound->pha[k]);
2822 csound->Free(csound, sound->smr[k]);
2823 csound->Free(csound, sound->res[k]);
2824 }
2825 /* pointers to data */
2826 csound->Free(csound, sound->time);
2827 csound->Free(csound, sound->frq);
2828 csound->Free(csound, sound->amp);
2829 csound->Free(csound, sound->pha);
2830 csound->Free(csound, sound->smr);
2831 csound->Free(csound, sound->res);
2832 /* check if we have residual data
2833 * and free its memory up
2834 */
2835 if (sound->band_energy != NULL) {
2836 for (k = 0; k < ATSA_CRITICAL_BANDS; k++)
2837 csound->Free(csound, sound->band_energy[k]);
2838 csound->Free(csound, sound->band_energy);
2839 }
2840 csound->Free(csound, sound);
2841 }
2842 }
2843
2844 /* module interface */
2845
atsa_init_(CSOUND * csound)2846 int atsa_init_(CSOUND *csound)
2847 {
2848 int retval = csound->AddUtility(csound, "atsa", atsa_main);
2849
2850 if (!retval) {
2851 retval =
2852 csound->SetUtilityDescription(csound, "atsa",
2853 Str("Soundfile analysis for ATS opcodes"));
2854 }
2855 return retval;
2856 }
2857
2858