1/* 2 * Copyright (C) 2000-2004 K�re Sj�lander <kare@speech.kth.se> 3 * Copyright (C) 1997 Philippe Langlais <felipe@speech.kth.se> 4 * 5 * This file is part of the Snack Sound Toolkit. 6 * The latest version can be found at http://www.speech.kth.se/snack/ 7 * 8 * This program is free software; you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation; either version 2 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program; if not, write to the Free Software 20 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 21 */ 22 23#include <stdlib.h> 24#include <stdio.h> 25#include <string.h> 26#include <math.h> 27#include "tcl.h" 28#include "snack.h" 29 30extern int Get_f0(Sound *s, Tcl_Interp *interp, int objc, 31 Tcl_Obj *CONST objv[]); 32 33/* ********************************************************************* */ 34/* LES PARAMETRES GLOBAUX DU DETECTEUR */ 35/* ********************************************************************* */ 36 37#define INFO 1 38 39/* #define CORRECTION pour faire une adaptation de la fo */ 40 41 42#define SEUIL_NRJ 40 /* on coupe si on a pas XX% de la dyn. de l'energie */ 43#define SEUIL_DPZ 50 /* on coupe si on a plus que YY% de la dyn de dpz */ 44#define EPSILON 10 /* seuil pour le calcul de la dpz */ 45 46#define SEUIL_VOISE 7 47#define AVANCE 2 48 49#define cst_pics_amdf 5 /* max. number of amdf values per trame */ 50#define PITCH_VARIABILITY 25 /* percentage of pitch-variability between 2 values */ 51#define FILTRE_PASSE_BAS 5 /* number of time it should apply the high frequency filter */ 52 53 54#define POURCENT(n,t) ( (t)? (((n)*100) /(t)) : 0 ) 55#define CARRE(a) ((a)*(a)) 56#define TO_FREQ(p) ((p)? (cst_freq_ech / (p)) : 0) 57#define minimum(a,b) (((a)>(b))? (b) : (a)) 58#define maximum(a,b) (((a)>(b))? (a) : (b)) 59 60#ifndef SEEK_SET 61#define SEEK_SET 0 62#endif 63#ifndef SEEK_END 64#define SEEK_END 2 65#endif 66 67#define PI_2 6.28318530717958 68#define MAX_ENTIER 2147483 69 70 71#define NON_VOISEE(t) (Vois[t] < SEUIL_VOISE) 72#define VOISEE(t) (Vois[t] >= SEUIL_VOISE) 73 74 75/* ********************************************************************* */ 76/* Quelques declarations anodines */ 77/* ********************************************************************* */ 78 79typedef struct 80 { 81 int total; 82 int rang; 83 } RESULT; 84 85typedef struct bidon 86 { 87 int debut,fin,ancrage; 88 struct bidon *suiv,*pred; 89 } *ZONE; 90 91 92/* ********************************************************************* */ 93/* LES VARIABLES GLOBALES DU DETECTEUR */ 94/* ********************************************************************* */ 95 96 97static int min_fo,max_fo,min_nrj,max_nrj,min_dpz,max_dpz; 98static int nb_coupe=0,debug,quick; 99static int seuil_dpz,seuil_nrj; 100 101/* 102 103 hamming window length = 2.5 * frequency / min_fo (16kHz,60Hz = 667) 104 window skip = frequency / 100 (16kHz = 160) 105 106 */ 107 108static int cst_freq_ech, 109 cst_freq_coupure, 110 cst_length_hamming, 111 cst_step_hamming, 112 cst_point_par_trame, 113 cst_step_min, 114 cst_step_max; 115 116 117static RESULT *(Coeff_Amdf[cst_pics_amdf]); 118static ZONE zone; 119static double *Hamming; 120static int max_amdf,min_amdf,amplitude_amdf; 121static short *Nrj,*Dpz,*Vois,*Fo; 122static float *Signal; 123static int **Resultat; 124 125 126/* *********************************************************************** */ 127static void init(int frequence,int cst_min_fo,int cst_max_fo) 128{ 129 cst_freq_coupure = 800; 130 cst_freq_ech = frequence; 131 cst_length_hamming = ( (int) (2.5 * cst_freq_ech) / cst_min_fo ); 132 cst_step_hamming = cst_point_par_trame = ((int) (cst_freq_ech / 100)); 133 cst_step_min = ( cst_freq_ech / cst_max_fo); 134 cst_step_max = ( cst_freq_ech / cst_min_fo); 135 136 if (debug > 1) 137 printf("sampling:%d, hamming size=%d, hamming overlap=%d\n", 138 cst_freq_ech,cst_length_hamming,cst_step_hamming); 139 } 140 141 /* *********************************************************************** */ 142/* 143static void filtre_passe_bas(int frequence,int longueur) 144 { 145 int i; 146 double coeff, 147 delai=0.0; 148 coeff = ( PI_2 * frequence ) / cst_freq_ech; 149 for (i=0;i<longueur;i++) 150 Signal[i] = (short) (delai = (double) (Signal[i] * coeff) + (delai * (1-coeff))); 151 } 152 */ 153/* *********************************************************************** */ 154 155static void libere_zone(ZONE zone) 156 { 157 ZONE l; 158 159 while (zone) 160 { 161 l=zone->suiv; 162 ckfree((char *) zone); 163 zone = l; 164 } 165 166 } 167 168/* *********************************************************************** */ 169static void libere_coeff_amdf(void) 170 { 171 int i; 172 for (i=0;i<cst_pics_amdf;i++) ckfree((char *) Coeff_Amdf[i]); 173 } 174 175/* *********************************************************************** */ 176 177static int voisement_par_profondeur_des_pics(int imin,int *result,int length) 178 { 179 int gauche,droite,i; 180 181 for (i=imin; i>0 && result[i] <= result[i-1]; i--); /* emergence gauche */ 182 gauche = result[i] - result[imin]; 183 for (i=imin; i<length-1 && result[i] <= result[i+1]; i++); /* emergence droite */ 184 droite = result[i] - result[imin]; 185 186 return minimum(droite,gauche); 187 } 188 189/* *********************************************************************** */ 190static void ranger_minimum(int trame, int valeur,int rang) 191 { 192 int i,j; 193 194 for (i=0; (i<cst_pics_amdf) && (valeur>=Coeff_Amdf[i][trame].total);i++); 195 if (i<cst_pics_amdf) 196 { 197 for (j=cst_pics_amdf-1;j>i;j--) 198 Coeff_Amdf[j][trame]=Coeff_Amdf[j-1][trame]; 199 Coeff_Amdf[i][trame].total = valeur; 200 Coeff_Amdf[i][trame].rang = rang; 201 } 202 } 203 204 205/* *********************************************************************** */ 206 207static void retiens_n_pics(int *result,int trame,int length,int *profondeur) 208 { 209#define homothetie(y,minval,amp) (int) ( (amp)? (((y)-(minval)) * 200 / (amp)) : 0 ) 210 211 int i,prof,profondeur_locale=0; 212 int *result_local; 213 int minVal,maxVal; 214 215 for (i=0;i<cst_pics_amdf;i++) /* ----- init ------ */ 216 { 217 Coeff_Amdf[i][trame].total = MAX_ENTIER ; 218 Coeff_Amdf[i][trame].rang = -1; 219 } 220 221 for (i=0; i<length-1; ) /* ----- minimum ----- */ 222 { 223 while ( (i<length-1) && (result[i]>result[i+1]) ) i++; /* on descend */ 224 if (i && i<length-1) ranger_minimum(trame,result[i],i+cst_step_min); 225 while ( (i<length-1) && (result[i]<=result[i+1]) ) i++; /* on monte */ 226 } 227 228 /* -------- Recherche du voisement par profondeur du pic minimum --------- */ 229 230 prof = *profondeur = 0; 231 232 { 233 result_local = (int *) ckalloc(sizeof(int)*length); 234 maxVal=0; 235 minVal = MAX_ENTIER; 236 for (i=0;i<length;i++) 237 { 238 if (result[i] > maxVal) maxVal = result[i]; 239 if (result[i] < minVal) minVal = result[i]; 240 } 241 if (debug > 1) printf("DYN AMDF[%d] : %d - %d (%d) ",trame,minVal,maxVal,maxVal-minVal); 242 } 243 244 for (i=0;i<length;i++) 245 { 246 result_local[i] = homothetie(result[i],minVal,maxVal-minVal); 247 result[i] = homothetie(result[i],min_amdf,amplitude_amdf); 248 } 249 250 for (i=0; i<cst_pics_amdf; i++) 251 if (Coeff_Amdf[i][trame].rang != -1) 252 { 253 prof = voisement_par_profondeur_des_pics( Coeff_Amdf[i][trame].rang-cst_step_min, 254 result, 255 length); 256 if ( prof > *profondeur) *profondeur = prof; 257 prof = voisement_par_profondeur_des_pics( Coeff_Amdf[i][trame].rang-cst_step_min, 258 result_local, 259 length); 260 if ( prof > profondeur_locale) profondeur_locale = prof; 261 } 262 263 Vois[trame] = profondeur_locale; 264 ckfree((char *) result_local); 265 if (debug>1) printf("----> %d\n",*profondeur); 266#undef homothetie 267 } 268 269/* *********************************************************************** */ 270 271static void calcul_voisement(int nb_trames) 272{ 273 int trame, profondeur, length; 274 275 amplitude_amdf = max_amdf - min_amdf; 276 277 length = cst_step_max-cst_step_min+1; 278 279 for (trame = 0; trame < nb_trames; trame++) 280 { 281 if (quick && (Nrj[trame] < seuil_nrj) && (Dpz[trame] > seuil_dpz)) Vois[trame] = 0; 282 else 283 { 284 retiens_n_pics(Resultat[trame],trame,length,&profondeur); 285 Vois[trame] = profondeur; 286 } 287 } 288 289 } 290 291/* *********************************************************************** */ 292 293static void precalcul_hamming(void) 294 { 295 double pas; 296 int i; 297 298 pas = PI_2 / cst_length_hamming; 299 for (i=0;i<cst_length_hamming;i++) 300 { 301 Hamming[i] = 0.54 - 0.46 * cos(i*pas); 302 } 303 } 304 305/*�������������������������������������������������������������������������*/ 306 307static void 308amdf(Sound *s, int i, int *Hammer, int *result, int nrj, int start) 309{ 310#define FACTEUR 50 311 312 int decal,j,l,somme,k; 313 double coeff, delai; 314 static double odelai[FILTRE_PASSE_BAS]; 315 316 Snack_GetSoundData(s, start+i, Signal, cst_length_hamming); 317 318 if (i == 0) { 319 for (k = 0; k < FILTRE_PASSE_BAS; k++) { 320 odelai[k] = 0.0; 321 } 322 } 323 324 for (k = 0; k < FILTRE_PASSE_BAS; k++) { 325 delai = odelai[k]; 326 coeff = ( PI_2 * cst_freq_coupure ) / cst_freq_ech; 327 for (j = 0; j < cst_length_hamming; j++) { 328 Signal[j] = (float) (delai = (double) (Signal[j] * coeff) 329 + (delai * (1-coeff))); 330 } 331 odelai[k] = (double) Signal[cst_step_hamming-1]; 332 } 333 334 for (j=0; j<cst_length_hamming; j++) { 335 Hammer[j] = (int) (Hamming[j]*Signal[j]); 336 } 337 338 /* ---- la double boucle couteuse --- */ 339 for (decal=cst_step_min; decal <= cst_step_max; decal++) { 340 for (somme=l=0,j=decal; j<cst_length_hamming; j++,l++) { 341 somme += abs( Hammer[j] - Hammer[l]); 342 } 343 result[decal-cst_step_min] = (int) ( (somme * FACTEUR) / ((cst_length_hamming-decal))); 344 } 345 346#undef FACTEUR 347 } 348/* *********************************************************************** */ 349static int 350parametre_amdf(Sound *s, Tcl_Interp *interp, int start, int longueur, 351 int *nb_trames,int *Hammer) 352{ 353 int j,i,length,trame; 354 355 max_amdf = 0; 356 min_amdf = MAX_ENTIER; 357 358 length = cst_step_max-cst_step_min+1; 359 360 for (trame=i=0; i < longueur; i += cst_step_hamming,trame++) { 361 if (i > s->length - cst_length_hamming) break; 362 if (i > longueur - cst_length_hamming/2) break; 363 if (quick && (Nrj[trame] < seuil_nrj) && (Dpz[trame] > seuil_dpz)) { 364 } else { 365 amdf(s, (int) i,Hammer,Resultat[trame],(Nrj[trame])? Nrj[trame]:1,start); 366 for (j=0; j<length;j++) { 367 if (Resultat[trame][j]>max_amdf) max_amdf = Resultat[trame][j]; 368 369 if (Resultat[trame][j]<min_amdf) min_amdf = Resultat[trame][j]; 370 } 371 } 372 if ((trame % 20) == 19) { 373 int res = Snack_ProgressCallback(s->cmdPtr, interp, "Computing pitch", 374 0.05 + 0.95 * (double) i / longueur); 375 if (res != TCL_OK) { 376 return TCL_ERROR; 377 } 378 } 379 } 380 Snack_ProgressCallback(s->cmdPtr, interp, "Computing pitch", 1.0); 381 *nb_trames = (int) trame; 382 383 if (debug) printf("min_amdf=%d, max_amdf=%d\n",min_amdf,max_amdf); 384 return TCL_OK; 385} 386 387/* ************************************************************************ */ 388/* LA MAGOUILLE DE LA RECHERCHE DU MEILLEUR CHEMIN */ 389/* ************************************************************************ */ 390 391static int interpolation(int x1,int y1,int x2,int y2,int x) 392 { 393 int a,b,delta; 394 395 /* y = ax+b */ 396 397 if (x==x1) return y1; 398 if (x==x2) return y2; 399 400 if ((delta = x1-x2)) 401 { 402 a = y1 -y2 ; 403 b = x1 * y2 - x2 * y1; 404 405 return (a*x +b) / delta; 406 } 407 return 0; 408 } 409/* *********************************************************************** */ 410 411static int ecart_correct(int avant,int apres) 412 { 413 int a,b; 414 415 a = cst_freq_ech / maximum(avant,apres); 416 b = cst_freq_ech / minimum(avant,apres); 417 418 419 return ( b <= ((a * (100+PITCH_VARIABILITY)) / 100) ) ; 420 421 } 422 423/* *********************************************************************** 424 selection du pic le plus proche d'une reference 425 ����������������������������������������������������������������������� */ 426 427void trier(int trame,int reference,RESULT *table) 428 { 429#define TABLEAU(t) (abs(table[t].rang-reference)) 430 int t,bulle=1; 431 432 433 for (t=0;t<cst_pics_amdf;t++) 434 table[t] = Coeff_Amdf[t][trame]; 435 436 while (bulle) 437 for (bulle=t=0;t<cst_pics_amdf-1;t++) 438 if ( (table[t].rang == -1 && table[t+1].rang != -1) || 439 (TABLEAU(t) > TABLEAU(t+1) && table[t+1].rang != -1) ) 440 { 441 RESULT temp; 442 temp = table[t+1]; 443 table[t+1] = table[t]; 444 table[t] = temp; 445 bulle = 1; 446 } 447 448#undef TABLEAU 449 } 450 451/* *********************************************************************** */ 452 453static void extension_zone(int accrochage,int debut,int fin,int to,short *tableau) 454 { 455#define POURCENTAGE 25 456#define SEUIL 20 457#define MAUVAIS_CHEMIN() 458#define ECART_CORRECT(a,b) ( (maximum(a,b) - minimum(a,b)) < minimum(SEUIL,(POURCENTAGE * minimum(a,b) / 100)) ) 459 460 int trame,avant,i,m,j; 461 RESULT table[AVANCE+1][cst_pics_amdf]; 462 RESULT normal[cst_pics_amdf]; 463 464 tableau[accrochage] = avant = to; 465 466/* ---------------on etend a gauche ---------------------------*/ 467 468 for (trame=accrochage;trame>debut;) 469 { 470 table[0][0].rang = avant; 471 for (i=1;i<=AVANCE && trame-i >=debut; i++) 472 trier(trame-i,table[i-1][0].rang,table[i]); 473 trier( m=maximum(trame-AVANCE , debut) ,avant,normal); 474 475 i--; 476 477 if ( table[i][0].rang != -1 && normal[0].rang != -1 && table[i][0].rang && normal[0].rang && avant && 478 abs(TO_FREQ(table[i][0].rang) - TO_FREQ(avant)) > abs(TO_FREQ(normal[0].rang) - TO_FREQ(avant)) 479 ) 480 { 481 for (j=1;j<=i; j++) 482 tableau[trame-j] = interpolation(trame,avant,m,normal[0].rang,trame-j); 483 trame -= i; 484 avant = normal[0].rang; 485 } 486 else 487 { 488 if ( table[1][0].rang && ECART_CORRECT(TO_FREQ(avant),TO_FREQ(table[1][0].rang))) 489 avant = tableau[--trame] = table[1][0].rang; 490 else tableau[--trame] = 0; 491 } 492 } 493 494 495/* ------------------ on etend a droite -----------------------*/ 496 497 avant = tableau[accrochage]; 498 for (trame=accrochage;trame<fin;) 499 { 500 table[0][0].rang = avant; 501 for (i=1;i<=AVANCE && trame+i <= fin; i++) 502 trier(trame+i,table[i-1][0].rang,table[i]); 503 trier( m = minimum(trame+AVANCE , fin) ,avant,normal); 504 505 i--; 506 507 if ( table[i][0].rang != -1 && normal[0].rang != -1 && table[i][0].rang && normal[0].rang && avant && 508 abs(TO_FREQ(table[i][0].rang) - TO_FREQ(avant)) > abs(TO_FREQ(normal[0].rang) - TO_FREQ(avant)) 509 ) 510 { 511 for (j=1;j<=i; j++) 512 tableau[trame+j] = interpolation(trame,avant,m,normal[0].rang,trame+j); 513 trame += i; 514 515 avant = normal[0].rang; 516 } 517 else 518 { 519 if ( table[1][0].rang && ECART_CORRECT(TO_FREQ(avant),TO_FREQ(table[1][0].rang))) 520 avant = tableau[++trame] = table[1][0].rang; 521 else tableau[++trame] = 0; 522 } 523 } 524 525#undef MAUVAIS_CHEMIN 526#undef SEUIL 527#undef ECART_CORRECT 528#undef POURCENTAGE 529 } 530/* *********************************************************************** */ 531 532static void recupere_trou(ZONE l,short *tableau) 533 { 534 int trame,avant; 535 536 537 for (trame=l->debut;trame<l->fin;) 538 { 539 for (; trame<=l->fin && tableau[trame];trame++); 540 for (avant=trame; trame<=l->fin && !tableau[trame];trame++); 541 while ( avant > l->debut && trame<l->fin && 542 (!tableau[trame] || !ecart_correct(tableau[avant-1],tableau[trame]) ) 543 ) trame++; 544 if (avant > l->debut && trame<l->fin && ecart_correct(tableau[avant-1],tableau[trame]) ) 545 { 546 int debut,fin,i; 547 548 debut = avant-1; 549 fin = trame; 550 for (i=debut+1;i<fin;i++) tableau[i] = interpolation(debut,tableau[debut],fin,tableau[fin],i); 551 } 552 } 553 554 } 555/* *********************************************************************** */ 556static int point_accrochage(int debut,int fin,int *accrochage,int To_Moyenne) 557{ 558#define POURCENTAGE 30 559#define PROCHE_MOYENNE(e,bande) \ 560 ( \ 561 ( (e) >= (To_Moyenne-bande) ) && \ 562 ( (e) <= (To_Moyenne+bande) ) \ 563 ) 564int bande,valeur,trame,pic; 565 566 bande = To_Moyenne * POURCENTAGE / 100; 567 valeur = MAX_ENTIER; 568 569 for (pic=0;pic<=1;pic++) 570 for (trame=debut;trame<=fin;trame++) 571 if ( abs(Coeff_Amdf[pic][trame].rang-To_Moyenne) < abs(valeur-To_Moyenne) ) 572 { 573 valeur = Coeff_Amdf[pic][trame].rang; 574 *accrochage = trame; 575 } 576 577 578 if (PROCHE_MOYENNE(valeur,bande)) return valeur; 579 580 return 0; 581 582#undef POURCENTAGE 583#undef PROCHE_MOYENNE 584} 585 586 587/* *********************************************************************** */ 588static int cherche_chemin(ZONE l,int To_moyen) 589 { 590 int accrochage; 591 592 if ( !( l->ancrage=point_accrochage(l->debut,l->fin,&accrochage,To_moyen) )) 593 return 0; 594 595 extension_zone(accrochage,l->debut,l->fin,l->ancrage,Fo); 596 recupere_trou(l,Fo); 597 return 1; 598 } 599 600/* ************************************************************************ */ 601 602static void calcul_courbe_fo(int nb_trames,int *To_Moyen) 603 { 604 int t,memo,debut; 605 ZONE l; 606 607 608 /* on met tout a zero au debut */ 609 610 memo = *To_Moyen; 611 612 for (debut=0; debut<nb_trames; debut++) Fo[debut] = 0; 613 614 for (l=zone;l;l=l->suiv) 615 { 616 if ( !cherche_chemin(l,*To_Moyen) ) 617 for (t=l->debut;t<=l->fin;t++) Fo[t]=0; 618 else 619 { 620#ifdef CORRECTION 621 int cum=0,nb=0; 622 623 for (t=l->fin; (t>= l->debut) ; t--) 624 { 625 cum += Fo[t]; 626 nb++; 627 } 628 *To_Moyen = ((2*cum) + (nb * memo)) / (3 * nb); 629 if (debug) printf("correction moyenne : %d\n",*To_Moyen); 630#endif 631 } 632 } 633 634 /* -------- on recalcule la moyenne sur le chemin choisi --------------*/ 635 636 min_fo = MAX_ENTIER; 637 max_fo = 0; 638 639 for ((*To_Moyen)=debut=t=0;t<nb_trames;t++) 640 if (Fo[t]) 641 { 642 (*To_Moyen) += Fo[t]; 643 Fo[t] = cst_freq_ech / Fo[t]; 644 if (max_fo < Fo[t]) max_fo = Fo[t]; 645 if (min_fo > Fo[t]) min_fo = Fo[t]; 646 debut++; 647 } 648 if (debut) (*To_Moyen) /= debut; 649 650 if (debug) printf("MOYENNE FINALE : %d (fo=%d)\n",*To_Moyen, cst_freq_ech / *To_Moyen); 651 } 652 653/* ********************************************************************** */ 654static void calcul_fo_moyen(int nb_trames, int *To_Moyenne) 655 { 656 int trame,nb,bulle; 657 RESULT *table; 658 int To_Moyenne_Corrigee; 659 660#define POURCENTAGE 30 661#define ACCEPTABLE(valeur,moyenne,pourcentage) \ 662 ( \ 663 ((valeur) > (moyenne-pourcentage)) && \ 664 ((valeur) < (moyenne+pourcentage)) \ 665 ) 666#define TABLEAU(t) (abs(table[t].rang-(*To_Moyenne))) 667 668 table = (RESULT *) ckalloc(sizeof(RESULT)*nb_trames); 669 670 for (*To_Moyenne=trame=nb=0;trame<nb_trames;trame++) 671 if ( VOISEE(trame) ) 672 { 673 table[nb++] = Coeff_Amdf[0][trame]; 674 (*To_Moyenne) += Coeff_Amdf[0][trame].rang; 675 } 676 677 *To_Moyenne = (nb)? ((*To_Moyenne) / nb) : 1; 678 679 if (debug) printf("To moyen non corrige : %d (fo=%d) \n",*To_Moyenne,cst_freq_ech / *To_Moyenne); 680 681 /* ------- correction de la valeur de fo ----------- */ 682 683 for (bulle=1;bulle;) 684 { 685 for (bulle=trame=0;trame<nb-1;trame++) 686 if (TABLEAU(trame)>TABLEAU(trame+1)) 687 { 688 RESULT temp; 689 bulle = 1; 690 temp = table[trame]; 691 table[trame] = table[trame+1]; 692 table[trame+1] = temp; 693 } 694 } 695 696 nb -= ((int ) POURCENTAGE * nb) / 100; 697 698 for (To_Moyenne_Corrigee=0,trame=0;trame<nb;trame++) 699 To_Moyenne_Corrigee += table[trame].rang; 700 701 To_Moyenne_Corrigee = (nb)? (To_Moyenne_Corrigee / nb) : 1; 702 703 /* -------- resultats ----------------------------------- */ 704 705 *To_Moyenne = To_Moyenne_Corrigee; 706 707 if (debug) printf("moyenne (a %d%% presque partout): %d (fo=%d)\n",100-POURCENTAGE, *To_Moyenne,cst_freq_ech / *To_Moyenne); 708 709 ckfree((char *) table); 710 711#undef POURCENTAGE 712#undef ACCEPTABLE 713#undef TABLEAU 714 } 715 716/* ************************************************************************ */ 717 718static ZONE calcul_zones_voisees(int nb_trames) 719 { 720 int trame,debut; 721 ZONE z = NULL; 722 723 724 for (trame=0;trame<nb_trames;) 725 { 726 while ( (trame<nb_trames) && (NON_VOISEE(trame)) ) trame++; 727 for (debut = trame; trame<nb_trames && VOISEE(trame);trame++ ); 728 729 if (trame<=nb_trames && debut<trame) 730 { 731 ZONE t,l,pred; 732 733 t = (ZONE) ckalloc(sizeof(struct bidon)); 734 t->debut = debut; 735 t->fin = trame-1; 736 t->ancrage = 0; 737 t->suiv = NULL; 738 739 for (l=z,pred=NULL; l;pred=l,l=l->suiv); 740 741 t->pred = pred; 742 if (pred) pred->suiv = t; 743 else z=t; 744 } 745 } 746 747 return z; 748 } 749/* *********************************************************************** */ 750static int calcul_nrj_dpz(Sound *s, Tcl_Interp *interp, int start,int longueur) 751{ 752 int J,JJ,m,i,j,trame,dpz,sens; 753 double nrj; 754 755 max_nrj = max_dpz = 0; 756 min_nrj = min_dpz = MAX_ENTIER; 757 758 Snack_ProgressCallback(s->cmdPtr, interp, "Computing pitch", 0.0); 759 760 for (trame=i=0; i<longueur; i += cst_step_hamming,trame++) { 761 J = minimum(s->length,(i+cst_length_hamming)); 762 JJ = J-1; 763 if (s->length < i + start + cst_length_hamming) { 764 Snack_GetSoundData(s, i+start, Signal, s->length - i + start); 765 for (j = s->length - i + start; j < cst_length_hamming; j++) { 766 Signal[j] = 0.0f; 767 } 768 } else { 769 Snack_GetSoundData(s, i+start, Signal, cst_length_hamming); 770 } 771 772 /* ---- nrj ---- */ 773 for (nrj=0.0,j=0; j<J-i; j++) 774 nrj += CARRE((double) Signal[j]); 775 776 m = Nrj[trame] = (short) (10 * log10(nrj)); 777 778 if (m > max_nrj) max_nrj = m; 779 if (m < min_nrj) min_nrj = m; 780 781 /* ---- dpz ---- */ 782 for (dpz=0,j=0; j<J-i; j++) { 783 while ( (j<J-i) && (abs((int)Signal[j]) > EPSILON) ) j++; /* looking just close to zero values */ 784 if (j<J-i) dpz ++; 785 sens = ( ((j-1) >= 0) && (Signal[j-1] > Signal[j])); 786 if (sens) while ( (j<JJ-i) && (Signal[j] > Signal[j+1]) ) j++; 787 else while ( (j<JJ-i) && (Signal[j] <= Signal[j+1]) ) j++; 788 } 789 m = Dpz[trame] = dpz; 790 791 if (m > max_dpz) max_dpz = m; 792 if (m < min_dpz) min_dpz = m; 793 794 if ((trame % 300) == 299) { 795 int res = Snack_ProgressCallback(s->cmdPtr, interp, "Computing pitch", 796 0.05 * (double) i / longueur); 797 if (res != TCL_OK) { 798 return TCL_ERROR; 799 } 800 } 801 } 802 803 seuil_nrj = min_nrj + (SEUIL_NRJ*(max_nrj-min_nrj))/100; 804 seuil_dpz = min_dpz + (SEUIL_DPZ*(max_dpz-min_dpz))/100; 805 806 if (debug) printf("dpz <%d,%d>, nrj <%d,%d> => Seuil nrj: %d, Seuil dpz: %d\n",min_dpz,max_dpz,min_nrj,max_nrj,seuil_nrj,seuil_dpz); 807 808 return trame; 809 } 810 811/* ********************************************************************** */ 812/* 813static void adjust_signal(int longueur, int freq, int debug) 814{ 815long i,moy,m; 816 817for (i=0; i<(m=minimum(freq,longueur)); moy += Signal[i++]); 818moy /= (m)? m:1; 819if (debug) fprintf(stderr,"ajustement de %d dans tout le signal\n",moy); 820if (moy) for (i=0; i<longueur; Signal[i++] -= (short) moy); 821} 822*/ 823/* ********************************************************************** */ 824 825int 826pitchCmd(Sound *s, Tcl_Interp *interp, int objc, Tcl_Obj *CONST objv[]) 827{ 828 int longueur, nb_trames, nb_trames_alloc, To_Moyen; 829 int *Hammer; 830 int i; 831 int fmin = 60, fmax = 400, lquick = 1/*, adjust = 0*/, nbframes; 832 short minnrj, maxnrj, minfo, maxfo, mindpz, maxdpz; 833 int arg, startpos = 0, endpos = -1, result, start; 834 Tcl_Obj *list; 835 static CONST84 char *subOptionStrings[] = { 836 "-start", "-end", "-maxpitch", "-minpitch", "-progress", 837 "-framelength", "-method", "-windowlength", NULL 838 }; 839 enum subOptions { 840 START, END, F0MAX, F0MIN, PROGRESS, FRAME, METHOD, WINLEN 841 }; 842 843 if (s->debug > 0) { Snack_WriteLog("Enter pitchCmd\n"); } 844 845 if (s->nchannels != 1) { 846 Tcl_AppendResult(interp, "pitch only works with Mono sounds", 847 (char *) NULL); 848 return TCL_ERROR; 849 } 850 851 for (arg = 2; arg < objc; arg += 2) { 852 char *opt = Tcl_GetStringFromObj(objv[arg], NULL); 853 char *val = Tcl_GetStringFromObj(objv[arg+1], NULL); 854 855 if ((strcmp("-method", opt) == 0) && (strcasecmp("esps", val) == 0)) { 856 Get_f0(s, interp, objc, objv); 857 return TCL_OK; 858 } 859 } 860 861 if (s->cmdPtr != NULL) { 862 Tcl_DecrRefCount(s->cmdPtr); 863 s->cmdPtr = NULL; 864 } 865 866 for (arg = 2; arg < objc; arg += 2) { 867 int index; 868 869 if (Tcl_GetIndexFromObj(interp, objv[arg], subOptionStrings, 870 "option", 0, &index) != TCL_OK) { 871 return TCL_ERROR; 872 } 873 874 if (arg + 1 == objc) { 875 Tcl_AppendResult(interp, "No argument given for ", 876 subOptionStrings[index], " option", (char *) NULL); 877 return TCL_ERROR; 878 } 879 880 switch ((enum subOptions) index) { 881 case START: 882 { 883 if (Tcl_GetIntFromObj(interp, objv[arg+1], &startpos) != TCL_OK) 884 return TCL_ERROR; 885 break; 886 } 887 case END: 888 { 889 if (Tcl_GetIntFromObj(interp, objv[arg+1], &endpos) != TCL_OK) 890 return TCL_ERROR; 891 break; 892 } 893 case F0MAX: 894 { 895 if (Tcl_GetIntFromObj(interp, objv[arg+1], &fmax) != TCL_OK) 896 return TCL_ERROR; 897 if (fmax <= 50) { 898 Tcl_AppendResult(interp, "Maximum pitch must be > 50", NULL); 899 return TCL_ERROR; 900 } 901 break; 902 } 903 case F0MIN: 904 { 905 if (Tcl_GetIntFromObj(interp, objv[arg+1], &fmin) != TCL_OK) 906 return TCL_ERROR; 907 if (fmin <= 50) { 908 Tcl_AppendResult(interp, "Minimum pitch must be > 50", NULL); 909 return TCL_ERROR; 910 } 911 break; 912 } 913 case PROGRESS: 914 { 915 char *str = Tcl_GetStringFromObj(objv[arg+1], NULL); 916 917 if (strlen(str) > 0) { 918 Tcl_IncrRefCount(objv[arg+1]); 919 s->cmdPtr = objv[arg+1]; 920 } 921 break; 922 } 923 case FRAME: 924 { 925 break; 926 } 927 case METHOD: 928 { 929 break; 930 } 931 case WINLEN: 932 { 933 break; 934 } 935 } 936 } 937 if (fmax <= fmin) { 938 Tcl_AppendResult(interp, "maxpitch must be > minpitch", NULL); 939 return TCL_ERROR; 940 } 941 if (startpos < 0) startpos = 0; 942 if (endpos >= (s->length - 1) || endpos == -1) 943 endpos = s->length - 1; 944 if (startpos > endpos) return TCL_OK; 945 946 quick = lquick; 947 init(s->samprate, fmin, fmax); 948 949 start = startpos - (cst_length_hamming / 2); 950 if (start < 0) start = 0; 951 if (endpos - start + 1 < cst_length_hamming) { 952 endpos = cst_length_hamming + start - 1; 953 if (endpos >= s->length) return TCL_OK; 954 } 955 longueur = endpos - start + 1; 956 957 if ((Signal = (float *) ckalloc(cst_length_hamming * sizeof(float))) == NULL) { 958 Tcl_AppendResult(interp, "Couldn't allocate buffer!", NULL); 959 return TCL_ERROR; 960 } 961 962 /* if (adjust) adjust_signal(longueur,s->samprate,s->debug);*/ 963 964 nb_trames = (longueur / cst_step_hamming) + 10; 965 Nrj = (short *) ckalloc(sizeof(short) * nb_trames); 966 Dpz = (short *) ckalloc(sizeof(short) * nb_trames); 967 Vois = (short *) ckalloc(sizeof(short) * nb_trames); 968 Fo = (short *) ckalloc(sizeof(short) * nb_trames); 969 Resultat = (int **) ckalloc(sizeof(int *) * nb_trames); 970 971 for (i = 0; i < nb_trames; i++) { 972 Resultat[i] = (int *) ckalloc(sizeof(int) * (cst_step_max-cst_step_min+1)); 973 } 974 nb_trames_alloc = nb_trames; 975 nb_trames = nbframes = calcul_nrj_dpz(s, interp, start, longueur); 976 977 Hamming = (double *) ckalloc(sizeof(double)*cst_length_hamming); 978 Hammer = (int *) ckalloc(sizeof(int) * cst_length_hamming); 979 for (i=0;i<cst_pics_amdf;i++) { 980 Coeff_Amdf[i] = (RESULT *) ckalloc(sizeof(RESULT)* nb_trames); 981 } 982 983 precalcul_hamming(); 984 985 result = parametre_amdf(s, interp, start, longueur, (int *)&(nbframes), 986 Hammer); 987 988 if (result == TCL_OK) { 989 if (debug) printf("nbframes=%d\n",nbframes); 990 calcul_voisement(nbframes); 991 zone = calcul_zones_voisees(nbframes); 992 calcul_fo_moyen(nbframes,&To_Moyen); 993 calcul_courbe_fo(nbframes,&To_Moyen); 994 995 minfo = (short) min_fo; 996 maxfo = (short) max_fo; 997 minnrj = (short) min_nrj; 998 maxnrj = (short) max_nrj; 999 mindpz = (short) min_dpz; 1000 maxdpz = (short) max_dpz; 1001 1002 if (debug && quick) 1003 printf("%d trames coupees sur %d -> %d %% (seuil nrj = %d, seuil dpz = %d) \n", 1004 nb_coupe,nbframes,POURCENT(nb_coupe,nbframes),seuil_nrj,seuil_dpz); 1005 libere_zone(zone); 1006 for (i=0;i<nb_trames_alloc;i++) 1007 if (Resultat[i]) ckfree((char *) Resultat[i]); 1008 } 1009 ckfree((char *) Hamming); 1010 ckfree((char *) Hammer); 1011 ckfree((char *)Signal); 1012 libere_coeff_amdf(); 1013 ckfree((char *) Resultat); 1014 if (result == TCL_OK) { 1015 int n = cst_length_hamming / (2 * cst_step_hamming) - startpos / cst_step_hamming; 1016 1017 list = Tcl_NewListObj(0, NULL); 1018 for (i = 0; i < n; i++) { 1019 Tcl_ListObjAppendElement(interp, list, Tcl_NewDoubleObj(0.0)); 1020 } 1021 for (i = 0; i < nbframes; i++) { 1022 Tcl_ListObjAppendElement(interp, list, Tcl_NewDoubleObj((double)Fo[i])); 1023 } 1024 Tcl_SetObjResult(interp, list); 1025 } 1026 ckfree((char *) Nrj); 1027 ckfree((char *) Dpz); 1028 ckfree((char *) Vois); 1029 ckfree((char *) Fo); 1030 1031 if (s->debug > 0) { Snack_WriteLog("Exit pitchCmd\n"); } 1032 1033 return TCL_OK; 1034} 1035 1036int 1037cPitch(Sound *s, Tcl_Interp *interp, int **outlist, int *length) 1038{ 1039 int longueur, nb_trames, To_Moyen; 1040 int *Hammer; 1041 int i; 1042 int fmin = 60, fmax = 400, lquick = 1/*, adjust = 0*/, nbframes; 1043 short minnrj, maxnrj, minfo, maxfo, mindpz, maxdpz; 1044 int startpos = 0, endpos = -1, result, start; 1045 1046 if (s->debug > 0) { Snack_WriteLog("Enter pitchCmd\n"); } 1047 1048 if (fmax <= fmin) { 1049 Tcl_AppendResult(interp, "maxpitch must be > minpitch", NULL); 1050 return TCL_ERROR; 1051 } 1052 if (startpos < 0) startpos = 0; 1053 if (endpos >= (s->length - 1) || endpos == -1) 1054 endpos = s->length - 1; 1055 if (startpos > endpos) return TCL_OK; 1056 1057 quick = lquick; 1058 init(s->samprate, fmin, fmax); 1059 1060 start = startpos - (cst_length_hamming / 2); 1061 if (start < 0) start = 0; 1062 longueur = endpos - start + 1; 1063 1064 if ((Signal = (float *) ckalloc(cst_length_hamming * sizeof(float))) == NULL) { 1065 Tcl_AppendResult(interp, "Couldn't allocate buffer!", NULL); 1066 return TCL_ERROR; 1067 } 1068 1069 nb_trames = (longueur / cst_step_hamming) + 10; 1070 Nrj = (short *) ckalloc(sizeof(short) * nb_trames); 1071 Dpz = (short *) ckalloc(sizeof(short) * nb_trames); 1072 Vois = (short *) ckalloc(sizeof(short) * nb_trames); 1073 Fo = (short *) ckalloc(sizeof(short) * nb_trames); 1074 Resultat = (int **) ckalloc(sizeof(int *) * nb_trames); 1075 1076 for (i = 0; i < nb_trames; i++) { 1077 Resultat[i] = (int *) ckalloc(sizeof(int) * (cst_step_max-cst_step_min+1)); 1078 } 1079 1080 nb_trames = nbframes = calcul_nrj_dpz(s, interp, start, longueur); 1081 1082 Hamming = (double *) ckalloc(sizeof(double)*cst_length_hamming); 1083 Hammer = (int *) ckalloc(sizeof(int) * cst_length_hamming); 1084 for (i=0;i<cst_pics_amdf;i++) { 1085 Coeff_Amdf[i] = (RESULT *) ckalloc(sizeof(RESULT)* nb_trames); 1086 } 1087 1088 precalcul_hamming(); 1089 1090 result = parametre_amdf(s, interp, start, longueur, (int *)&(nbframes), 1091 Hammer); 1092 1093 if (result == TCL_OK) { 1094 calcul_voisement(nbframes); 1095 zone = calcul_zones_voisees(nbframes); 1096 calcul_fo_moyen(nbframes,&To_Moyen); 1097 calcul_courbe_fo(nbframes,&To_Moyen); 1098 1099 minfo = (short) min_fo; 1100 maxfo = (short) max_fo; 1101 minnrj = (short) min_nrj; 1102 maxnrj = (short) max_nrj; 1103 mindpz = (short) min_dpz; 1104 maxdpz = (short) max_dpz; 1105 1106 libere_zone(zone); 1107 for (i=0;i<nbframes;i++) 1108 if (Resultat[i]) ckfree((char *) Resultat[i]); 1109 } 1110 ckfree((char *) Hamming); 1111 ckfree((char *) Hammer); 1112 ckfree((char *)Signal); 1113 libere_coeff_amdf(); 1114 ckfree((char *) Resultat); 1115 if (result == TCL_OK) { 1116 int n = cst_length_hamming / (2 * cst_step_hamming) - startpos / cst_step_hamming; 1117 int *tmp = (int *)ckalloc(sizeof(int) * (n + nb_trames)); 1118 for (i = 0; i < n; i++) { 1119 tmp[i] = 0; 1120 } 1121 for (i = n; i < nbframes + n; i++) { 1122 tmp[i] = Fo[i-n]; 1123 } 1124 *outlist = tmp; 1125 *length = nbframes + n; 1126 } 1127 ckfree((char *) Nrj); 1128 ckfree((char *) Dpz); 1129 ckfree((char *) Vois); 1130 ckfree((char *) Fo); 1131 1132 if (s->debug > 0) { Snack_WriteLog("Exit pitchCmd\n"); } 1133 1134 return TCL_OK; 1135} 1136