comparison FChTransformF0gram.cpp @ 10:af59167b3d35 perf

size_t -> int throughout: permits better optimisation for tight for-loops, making the whole thing about 10% faster
author Chris Cannam
date Tue, 02 Oct 2018 13:21:15 +0100
parents 54dd64b6cfc0
children fc8f351d2cd6
comparison
equal deleted inserted replaced
9:54dd64b6cfc0 10:af59167b3d35
476 // See OutputDescriptor documentation for the possibilities here. 476 // See OutputDescriptor documentation for the possibilities here.
477 // Every plugin must have at least one output. 477 // Every plugin must have at least one output.
478 478
479 /* f0 values of F0gram grid as string values */ 479 /* f0 values of F0gram grid as string values */
480 vector<string> f0values; 480 vector<string> f0values;
481 size_t ind = 0; 481 int ind = 0;
482 char f0String[10]; 482 char f0String[10];
483 while (ind < m_num_f0s) { 483 while (ind < m_num_f0s) {
484 sprintf(f0String, "%4.2f", m_f0s[ind]); 484 sprintf(f0String, "%4.2f", m_f0s[ind]);
485 f0values.push_back(f0String); 485 f0values.push_back(f0String);
486 ind++; 486 ind++;
550 550
551 void 551 void
552 FChTransformF0gram::design_GLogS() { 552 FChTransformF0gram::design_GLogS() {
553 553
554 // total number & initial quantity of f0s 554 // total number & initial quantity of f0s
555 m_glogs_init_f0s = (size_t)(((double)m_f0_params.num_f0s_per_oct)*log2(5.0))+1; 555 m_glogs_init_f0s = (int)(((double)m_f0_params.num_f0s_per_oct)*log2(5.0))+1;
556 m_glogs_num_f0s = (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct + m_glogs_init_f0s; 556 m_glogs_num_f0s = (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct + m_glogs_init_f0s;
557 557
558 // Initialize arrays 558 // Initialize arrays
559 m_glogs_f0 = new double[m_glogs_num_f0s]; 559 m_glogs_f0 = new double[m_glogs_num_f0s];
560 m_glogs = new double[m_glogs_num_f0s*m_warp_params.num_warps]; 560 m_glogs = new double[m_glogs_num_f0s*m_warp_params.num_warps];
561 m_glogs_n = new size_t[m_glogs_num_f0s]; 561 m_glogs_n = new int[m_glogs_num_f0s];
562 m_glogs_index = new size_t[m_glogs_num_f0s]; 562 m_glogs_index = new int[m_glogs_num_f0s];
563 563
564 // Compute f0 values 564 // Compute f0 values
565 m_glogs_harmonic_count = 0; 565 m_glogs_harmonic_count = 0;
566 double factor = (double)(m_warp_params.nsamps_twarp/2)/(double)(m_warp_params.nsamps_twarp/2+1); 566 double factor = (double)(m_warp_params.nsamps_twarp/2)/(double)(m_warp_params.nsamps_twarp/2+1);
567 for (size_t i = 0; i < m_glogs_num_f0s; i++) { 567 for (int i = 0; i < m_glogs_num_f0s; i++) {
568 m_glogs_f0[i] = (m_f0_params.f0min/5.0)*pow(2.0,(double)i/(double)m_f0_params.num_f0s_per_oct); 568 m_glogs_f0[i] = (m_f0_params.f0min/5.0)*pow(2.0,(double)i/(double)m_f0_params.num_f0s_per_oct);
569 // for every f0 compute number of partials less or equal than m_fmax. 569 // for every f0 compute number of partials less or equal than m_fmax.
570 m_glogs_n[i] = m_fmax*factor/m_glogs_f0[i]; 570 m_glogs_n[i] = m_fmax*factor/m_glogs_f0[i];
571 m_glogs_index[i] = m_glogs_harmonic_count; 571 m_glogs_index[i] = m_glogs_harmonic_count;
572 m_glogs_harmonic_count += m_glogs_n[i]; 572 m_glogs_harmonic_count += m_glogs_n[i];
573 } 573 }
574 574
575 // Initialize arrays for interpolation 575 // Initialize arrays for interpolation
576 m_glogs_posint = new size_t[m_glogs_harmonic_count]; 576 m_glogs_posint = new int[m_glogs_harmonic_count];
577 m_glogs_posfrac = new double[m_glogs_harmonic_count]; 577 m_glogs_posfrac = new double[m_glogs_harmonic_count];
578 m_glogs_interp = new double[m_glogs_harmonic_count]; 578 m_glogs_interp = new double[m_glogs_harmonic_count];
579 579
580 // Compute int & frac of interpolation positions 580 // Compute int & frac of interpolation positions
581 size_t aux_index = 0; 581 int aux_index = 0;
582 double aux_pos; 582 double aux_pos;
583 for (size_t i = 0; i < m_glogs_num_f0s; i++) { 583 for (int i = 0; i < m_glogs_num_f0s; i++) {
584 for (size_t j = 1; j <= m_glogs_n[i]; j++) { 584 for (int j = 1; j <= m_glogs_n[i]; j++) {
585 // indice en el vector de largo t_warp/2+1 donde el ultimo valor corresponde a f=m_fmax 585 // indice en el vector de largo t_warp/2+1 donde el ultimo valor corresponde a f=m_fmax
586 aux_pos = ((double)j*m_glogs_f0[i])*((double)(m_warp_params.nsamps_twarp/2+1))/m_fmax; 586 aux_pos = ((double)j*m_glogs_f0[i])*((double)(m_warp_params.nsamps_twarp/2+1))/m_fmax;
587 m_glogs_posint[aux_index] = (size_t)aux_pos; 587 m_glogs_posint[aux_index] = (int)aux_pos;
588 m_glogs_posfrac[aux_index] = aux_pos - (double)m_glogs_posint[aux_index]; 588 m_glogs_posfrac[aux_index] = aux_pos - (double)m_glogs_posint[aux_index];
589 aux_index++; 589 aux_index++;
590 } 590 }
591 } 591 }
592 592
593 // Third harmonic attenuation 593 // Third harmonic attenuation
594 double aux_third_harmonic; 594 double aux_third_harmonic;
595 m_glogs_third_harmonic_posint = new size_t[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct]; 595 m_glogs_third_harmonic_posint = new int[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
596 m_glogs_third_harmonic_posfrac = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct]; 596 m_glogs_third_harmonic_posfrac = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
597 for (size_t i = 0; i < (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct; i++) { 597 for (int i = 0; i < (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct; i++) {
598 aux_third_harmonic = (double)i + (double)m_glogs_init_f0s - ((double)m_f0_params.num_f0s_per_oct)*log2(3.0); 598 aux_third_harmonic = (double)i + (double)m_glogs_init_f0s - ((double)m_f0_params.num_f0s_per_oct)*log2(3.0);
599 m_glogs_third_harmonic_posint[i] = (size_t)aux_third_harmonic; 599 m_glogs_third_harmonic_posint[i] = (int)aux_third_harmonic;
600 m_glogs_third_harmonic_posfrac[i] = aux_third_harmonic - (double)(m_glogs_third_harmonic_posint[i]); 600 m_glogs_third_harmonic_posfrac[i] = aux_third_harmonic - (double)(m_glogs_third_harmonic_posint[i]);
601 } 601 }
602 m_glogs_third_harmonic = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct]; 602 m_glogs_third_harmonic = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
603 603
604 // Fifth harmonic attenuation 604 // Fifth harmonic attenuation
605 double aux_fifth_harmonic; 605 double aux_fifth_harmonic;
606 m_glogs_fifth_harmonic_posint = new size_t[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct]; 606 m_glogs_fifth_harmonic_posint = new int[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
607 m_glogs_fifth_harmonic_posfrac = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct]; 607 m_glogs_fifth_harmonic_posfrac = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
608 for (size_t i = 0; i < (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct; i++) { 608 for (int i = 0; i < (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct; i++) {
609 aux_fifth_harmonic = (double)i + (double)m_glogs_init_f0s - ((double)m_f0_params.num_f0s_per_oct)*log2(5.0); 609 aux_fifth_harmonic = (double)i + (double)m_glogs_init_f0s - ((double)m_f0_params.num_f0s_per_oct)*log2(5.0);
610 m_glogs_fifth_harmonic_posint[i] = (size_t)aux_fifth_harmonic; 610 m_glogs_fifth_harmonic_posint[i] = (int)aux_fifth_harmonic;
611 m_glogs_fifth_harmonic_posfrac[i] = aux_fifth_harmonic - (double)(m_glogs_fifth_harmonic_posint[i]); 611 m_glogs_fifth_harmonic_posfrac[i] = aux_fifth_harmonic - (double)(m_glogs_fifth_harmonic_posint[i]);
612 } 612 }
613 m_glogs_fifth_harmonic = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct]; 613 m_glogs_fifth_harmonic = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
614 614
615 // Normalization & attenuation windows 615 // Normalization & attenuation windows
616 m_glogs_f0_preference_weights = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct]; 616 m_glogs_f0_preference_weights = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct];
617 m_glogs_median_correction = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct]; 617 m_glogs_median_correction = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct];
618 m_glogs_sigma_correction = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct]; 618 m_glogs_sigma_correction = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct];
619 m_glogs_hf_smoothing_window = new double[m_warp_params.nsamps_twarp/2+1]; 619 m_glogs_hf_smoothing_window = new double[m_warp_params.nsamps_twarp/2+1];
620 double MIDI_value; 620 double MIDI_value;
621 for (size_t i = 0; i < m_f0_params.num_octs*m_f0_params.num_f0s_per_oct; i++) { 621 for (int i = 0; i < m_f0_params.num_octs*m_f0_params.num_f0s_per_oct; i++) {
622 MIDI_value = 69.0 + 12.0 * log2(m_glogs_f0[i + m_glogs_init_f0s]/440.0); 622 MIDI_value = 69.0 + 12.0 * log2(m_glogs_f0[i + m_glogs_init_f0s]/440.0);
623 m_glogs_f0_preference_weights[i] = 1.0/sqrt(2.0*M_PI*m_f0_params.prefer_stdev*m_f0_params.prefer_stdev)*exp(-(MIDI_value-m_f0_params.prefer_mean)*(MIDI_value-m_f0_params.prefer_mean)/(2.0*m_f0_params.prefer_stdev*m_f0_params.prefer_stdev)); 623 m_glogs_f0_preference_weights[i] = 1.0/sqrt(2.0*M_PI*m_f0_params.prefer_stdev*m_f0_params.prefer_stdev)*exp(-(MIDI_value-m_f0_params.prefer_mean)*(MIDI_value-m_f0_params.prefer_mean)/(2.0*m_f0_params.prefer_stdev*m_f0_params.prefer_stdev));
624 m_glogs_f0_preference_weights[i] = (0.01 + m_glogs_f0_preference_weights[i]) / (1.01); 624 m_glogs_f0_preference_weights[i] = (0.01 + m_glogs_f0_preference_weights[i]) / (1.01);
625 625
626 m_glogs_median_correction[i] = m_glogs_params.median_poly_coefs[0]*(i+1.0)*(i+1.0) + m_glogs_params.median_poly_coefs[1]*(i+1.0) + m_glogs_params.median_poly_coefs[2]; 626 m_glogs_median_correction[i] = m_glogs_params.median_poly_coefs[0]*(i+1.0)*(i+1.0) + m_glogs_params.median_poly_coefs[1]*(i+1.0) + m_glogs_params.median_poly_coefs[2];
627 m_glogs_sigma_correction[i] = 1.0 / (m_glogs_params.sigma_poly_coefs[0]*(i+1.0)*(i+1.0) + m_glogs_params.sigma_poly_coefs[1]*(i+1.0) + m_glogs_params.sigma_poly_coefs[2]); 627 m_glogs_sigma_correction[i] = 1.0 / (m_glogs_params.sigma_poly_coefs[0]*(i+1.0)*(i+1.0) + m_glogs_params.sigma_poly_coefs[1]*(i+1.0) + m_glogs_params.sigma_poly_coefs[2]);
628 } 628 }
629 629
630 double smooth_width = 1000.0; // hertz. 630 double smooth_width = 1000.0; // hertz.
631 double smooth_aux = (double)(m_warp_params.nsamps_twarp/2+1)*(m_fmax-smooth_width)/m_fmax; 631 double smooth_aux = (double)(m_warp_params.nsamps_twarp/2+1)*(m_fmax-smooth_width)/m_fmax;
632 for (size_t i = 0; i < m_warp_params.nsamps_twarp/2+1; i++) { 632 for (int i = 0; i < m_warp_params.nsamps_twarp/2+1; i++) {
633 if (i < smooth_aux) { 633 if (i < smooth_aux) {
634 m_glogs_hf_smoothing_window[i] = 1.0; 634 m_glogs_hf_smoothing_window[i] = 1.0;
635 } else { 635 } else {
636 m_glogs_hf_smoothing_window[i] = ((double)i - (double)m_warp_params.nsamps_twarp/2.0)*(-1.0/((double)(m_warp_params.nsamps_twarp/2+1)-smooth_aux)); 636 m_glogs_hf_smoothing_window[i] = ((double)i - (double)m_warp_params.nsamps_twarp/2.0)*(-1.0/((double)(m_warp_params.nsamps_twarp/2+1)-smooth_aux));
637 } 637 }
670 // equivalent to: m_warpings.nsamps_torig = m_warp_params.fact_over_samp * m_blockSize; 670 // equivalent to: m_warpings.nsamps_torig = m_warp_params.fact_over_samp * m_blockSize;
671 671
672 // time instants of the original signal frame 672 // time instants of the original signal frame
673 double t_orig[m_warpings.nsamps_torig]; 673 double t_orig[m_warpings.nsamps_torig];
674 //float * t_orig = new float [m_warpings.nsamps_torig]; 674 //float * t_orig = new float [m_warpings.nsamps_torig];
675 for (size_t ind = 0; ind < m_warpings.nsamps_torig; ind++) { 675 for (int ind = 0; ind < m_warpings.nsamps_torig; ind++) {
676 t_orig[ind] = ((double)(ind + 1) - (double)m_warpings.nsamps_torig / 2.0) / m_warpings.fs_orig; 676 t_orig[ind] = ((double)(ind + 1) - (double)m_warpings.nsamps_torig / 2.0) / m_warpings.fs_orig;
677 } 677 }
678 678
679 // linear chirps warping definition as relative frequency deviation 679 // linear chirps warping definition as relative frequency deviation
680 //double * freq_relative = new double [m_warpings.nsamps_torig * m_warp_params.num_warps]; 680 //double * freq_relative = new double [m_warpings.nsamps_torig * m_warp_params.num_warps];
682 double *freq_relative = new double [m_warpings.nsamps_torig * m_warp_params.num_warps]; 682 double *freq_relative = new double [m_warpings.nsamps_torig * m_warp_params.num_warps];
683 define_warps_linear_chirps(freq_relative, t_orig); 683 define_warps_linear_chirps(freq_relative, t_orig);
684 684
685 // maximum relative frequency deviation 685 // maximum relative frequency deviation
686 double freq_relative_max = 0; 686 double freq_relative_max = 0;
687 for (size_t i = 0; i < m_warpings.nsamps_torig; i++) 687 for (int i = 0; i < m_warpings.nsamps_torig; i++)
688 for (size_t j = 0; j < m_warp_params.num_warps; j++) 688 for (int j = 0; j < m_warp_params.num_warps; j++)
689 if (freq_relative_max < freq_relative[j * m_warpings.nsamps_torig + i]) 689 if (freq_relative_max < freq_relative[j * m_warpings.nsamps_torig + i])
690 freq_relative_max = freq_relative[j * m_warpings.nsamps_torig + i]; 690 freq_relative_max = freq_relative[j * m_warpings.nsamps_torig + i];
691 691
692 // sampling frequency of warped signal to be free of aliasing up to fmax 692 // sampling frequency of warped signal to be free of aliasing up to fmax
693 m_warpings.fs_warp = 2 * m_fmax * freq_relative_max; 693 m_warpings.fs_warp = 2 * m_fmax * freq_relative_max;
694 694
695 // time instants of the warped signal frame 695 // time instants of the warped signal frame
696 double t_warp[m_warp_params.nsamps_twarp]; 696 double t_warp[m_warp_params.nsamps_twarp];
697 for (size_t ind = 0; ind < m_warp_params.nsamps_twarp; ind++) { 697 for (int ind = 0; ind < m_warp_params.nsamps_twarp; ind++) {
698 t_warp[ind] = ((double)((int)(ind + 1)- (int)m_warp_params.nsamps_twarp / 2)) / (double)m_warpings.fs_warp; 698 t_warp[ind] = ((double)((int)(ind + 1)- (int)m_warp_params.nsamps_twarp / 2)) / (double)m_warpings.fs_warp;
699 } 699 }
700 700
701 // design of warpings for efficient interpolation 701 // design of warpings for efficient interpolation
702 design_warps(freq_relative, t_orig, t_warp); 702 design_warps(freq_relative, t_orig, t_warp);
706 * FILES FOR DEBUGGING 706 * FILES FOR DEBUGGING
707 */ 707 */
708 708
709 /* 709 /*
710 output << "chirp_rates" << endl; 710 output << "chirp_rates" << endl;
711 for (size_t j = 0; j < m_warp_params.num_warps; j++){ 711 for (int j = 0; j < m_warp_params.num_warps; j++){
712 output << m_warpings.chirp_rates[j]; 712 output << m_warpings.chirp_rates[j];
713 output << " "; 713 output << " ";
714 } 714 }
715 output << endl << "freq_relative" << endl; 715 output << endl << "freq_relative" << endl;
716 716
717 for (size_t i = 0; i < m_warpings.nsamps_torig; i++){ 717 for (int i = 0; i < m_warpings.nsamps_torig; i++){
718 for (size_t j = 0; j < m_warp_params.num_warps; j++){ 718 for (int j = 0; j < m_warp_params.num_warps; j++){
719 output << freq_relative[j * m_warpings.nsamps_torig + i]; 719 output << freq_relative[j * m_warpings.nsamps_torig + i];
720 output << " "; 720 output << " ";
721 } 721 }
722 output << endl; 722 output << endl;
723 } 723 }
724 724
725 output << endl << "t_orig" << endl; 725 output << endl << "t_orig" << endl;
726 726
727 for (size_t i = 0; i < m_warpings.nsamps_torig; i++){ 727 for (int i = 0; i < m_warpings.nsamps_torig; i++){
728 output << t_orig[i] << endl ; 728 output << t_orig[i] << endl ;
729 } 729 }
730 */ 730 */
731 731
732 delete [] freq_relative; 732 delete [] freq_relative;
747 given by the desired frequency deviation, to do this, the interpolation 747 given by the desired frequency deviation, to do this, the interpolation
748 instants are stored in a structure as an integer index and a fractional value 748 instants are stored in a structure as an integer index and a fractional value
749 hypothesis: sampling frequency at the central point equals the original 749 hypothesis: sampling frequency at the central point equals the original
750 */ 750 */
751 751
752 m_warpings.pos_int = new size_t[m_warp_params.num_warps * m_warp_params.nsamps_twarp]; 752 m_warpings.pos_int = new int[m_warp_params.num_warps * m_warp_params.nsamps_twarp];
753 m_warpings.pos_frac = new double[m_warp_params.num_warps * m_warp_params.nsamps_twarp]; 753 m_warpings.pos_frac = new double[m_warp_params.num_warps * m_warp_params.nsamps_twarp];
754 754
755 // vector of phase values 755 // vector of phase values
756 double *phi = new double[m_warpings.nsamps_torig]; 756 double *phi = new double[m_warpings.nsamps_torig];
757 double aux; 757 double aux;
758 758
759 // warped positions 759 // warped positions
760 double *pos1 = new double[m_warp_params.nsamps_twarp*m_warp_params.num_warps]; 760 double *pos1 = new double[m_warp_params.nsamps_twarp*m_warp_params.num_warps];
761 761
762 for (size_t i = 0; i < m_warp_params.num_warps; i++) { 762 for (int i = 0; i < m_warp_params.num_warps; i++) {
763 763
764 // integration of relative frequency to obtain phase values 764 // integration of relative frequency to obtain phase values
765 cumtrapz(t_orig, freq_relative + i*(m_warpings.nsamps_torig), m_warpings.nsamps_torig, phi); 765 cumtrapz(t_orig, freq_relative + i*(m_warpings.nsamps_torig), m_warpings.nsamps_torig, phi);
766 766
767 // centering of phase values to force original frequency in the middle 767 // centering of phase values to force original frequency in the middle
768 aux = phi[m_warpings.nsamps_torig/2]; 768 aux = phi[m_warpings.nsamps_torig/2];
769 for (size_t j = 0; j < m_warpings.nsamps_torig; j++) { 769 for (int j = 0; j < m_warpings.nsamps_torig; j++) {
770 phi[j] -= aux; 770 phi[j] -= aux;
771 } //for 771 } //for
772 772
773 // interpolation of phase values to obtain warped positions 773 // interpolation of phase values to obtain warped positions
774 interp1(phi, t_orig, m_warpings.nsamps_torig, t_warp, pos1 + i*m_warp_params.nsamps_twarp, m_warp_params.nsamps_twarp); 774 interp1(phi, t_orig, m_warpings.nsamps_torig, t_warp, pos1 + i*m_warp_params.nsamps_twarp, m_warp_params.nsamps_twarp);
779 // % integer corresponding to previous sample index in "c" 779 // % integer corresponding to previous sample index in "c"
780 // warps.pos1_int = (pos1_int - uint32(1)); 780 // warps.pos1_int = (pos1_int - uint32(1));
781 // % fractional value that defines the warped position 781 // % fractional value that defines the warped position
782 // warps.pos1_frac = (double(pos1)' - double(pos1_int)); 782 // warps.pos1_frac = (double(pos1)' - double(pos1_int));
783 783
784 for (size_t j = 0; j < m_warp_params.nsamps_twarp*m_warp_params.num_warps; j++) { 784 for (int j = 0; j < m_warp_params.nsamps_twarp*m_warp_params.num_warps; j++) {
785 // previous sample index 785 // previous sample index
786 pos1[j] = pos1[j]*m_warpings.fs_orig + m_warpings.nsamps_torig/2 + 1; 786 pos1[j] = pos1[j]*m_warpings.fs_orig + m_warpings.nsamps_torig/2 + 1;
787 m_warpings.pos_int[j] = (size_t) pos1[j]; 787 m_warpings.pos_int[j] = (int) pos1[j];
788 m_warpings.pos_frac[j] = pos1[j] - (double)(m_warpings.pos_int[j]); 788 m_warpings.pos_frac[j] = pos1[j] - (double)(m_warpings.pos_int[j]);
789 } //for 789 } //for
790 790
791 delete [] phi; 791 delete [] phi;
792 delete [] pos1; 792 delete [] pos1;
804 m_warpings.chirp_rates = new double [m_warp_params.num_warps]; 804 m_warpings.chirp_rates = new double [m_warp_params.num_warps];
805 // WARNING m_warp_params.num_warps must be odd 805 // WARNING m_warp_params.num_warps must be odd
806 m_warpings.chirp_rates[0] = -m_warp_params.alpha_max; 806 m_warpings.chirp_rates[0] = -m_warp_params.alpha_max;
807 double increment = (double) m_warp_params.alpha_max / ((m_warp_params.num_warps - 1) / 2); 807 double increment = (double) m_warp_params.alpha_max / ((m_warp_params.num_warps - 1) / 2);
808 808
809 for (size_t ind = 1; ind < m_warp_params.num_warps; ind++) { 809 for (int ind = 1; ind < m_warp_params.num_warps; ind++) {
810 m_warpings.chirp_rates[ind] = m_warpings.chirp_rates[ind - 1] + increment; 810 m_warpings.chirp_rates[ind] = m_warpings.chirp_rates[ind - 1] + increment;
811 } 811 }
812 // force zero value 812 // force zero value
813 m_warpings.chirp_rates[(int) ((m_warp_params.num_warps - 1) / 2)] = 0; 813 m_warpings.chirp_rates[(int) ((m_warp_params.num_warps - 1) / 2)] = 0;
814 814
824 double increment = logMax / ((m_warp_params.num_warps - 1) / 2.0f); 824 double increment = logMax / ((m_warp_params.num_warps - 1) / 2.0f);
825 double exponent = 0; 825 double exponent = 0;
826 826
827 // fill positive values 827 // fill positive values
828 int ind_log = middle_point; 828 int ind_log = middle_point;
829 for (size_t ind = 0; ind < (m_warp_params.num_warps + 1) / 2; ind++) { 829 for (int ind = 0; ind < (m_warp_params.num_warps + 1) / 2; ind++) {
830 m_warpings.chirp_rates[ind_log] = pow(10, exponent) - 1; 830 m_warpings.chirp_rates[ind_log] = pow(10, exponent) - 1;
831 exponent += increment; 831 exponent += increment;
832 ind_log++; 832 ind_log++;
833 } 833 }
834 // fill negative values 834 // fill negative values
835 for (size_t ind = 0; ind < (m_warp_params.num_warps - 1) / 2; ind++) { 835 for (int ind = 0; ind < (m_warp_params.num_warps - 1) / 2; ind++) {
836 m_warpings.chirp_rates[ind] = -m_warpings.chirp_rates[m_warp_params.num_warps - 1 - ind]; 836 m_warpings.chirp_rates[ind] = -m_warpings.chirp_rates[m_warp_params.num_warps - 1 - ind];
837 } 837 }
838 } 838 }
839 839
840 // compute relative frequency deviation 840 // compute relative frequency deviation
841 for (size_t i = 0; i < m_warpings.nsamps_torig; i++) 841 for (int i = 0; i < m_warpings.nsamps_torig; i++)
842 for (size_t j = 0; j < m_warp_params.num_warps; j++) 842 for (int j = 0; j < m_warp_params.num_warps; j++)
843 freq_relative[j * m_warpings.nsamps_torig + i] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j]; 843 freq_relative[j * m_warpings.nsamps_torig + i] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j];
844 //freq_relative[i * m_warpings.nsamps_torig + j] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j]; 844 //freq_relative[i * m_warpings.nsamps_torig + j] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j];
845 //freq_relative[i][j] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j]; 845 //freq_relative[i][j] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j];
846 } 846 }
847 847
853 // in_window = (float*) fftw_malloc(sizeof (float) * tamanoVentana); 853 // in_window = (float*) fftw_malloc(sizeof (float) * tamanoVentana);
854 // p = fftw_plan_dft_1d(tamanoVentana, in, out, FFTW_FORWARD, FFTW_ESTIMATE); 854 // p = fftw_plan_dft_1d(tamanoVentana, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
855 double *lp_LPFWindow_aux = new double[m_blockSize/2+1]; 855 double *lp_LPFWindow_aux = new double[m_blockSize/2+1];
856 mp_LPFWindow = new double[m_blockSize/2+1]; 856 mp_LPFWindow = new double[m_blockSize/2+1];
857 857
858 size_t i_max = (size_t) ((2.0*m_fmax/m_fs) * ( (double)m_blockSize / 2.0 + 1.0 )); 858 int i_max = (int) ((2.0*m_fmax/m_fs) * ( (double)m_blockSize / 2.0 + 1.0 ));
859 for (size_t i = 0; i < m_blockSize/2+1; i++) { 859 for (int i = 0; i < m_blockSize/2+1; i++) {
860 if (i >= i_max) { 860 if (i >= i_max) {
861 lp_LPFWindow_aux[i] = 0.0; 861 lp_LPFWindow_aux[i] = 0.0;
862 } else { 862 } else {
863 lp_LPFWindow_aux[i] = 1.0; 863 lp_LPFWindow_aux[i] = 1.0;
864 } 864 }
865 } 865 }
866 LPF_time = (double*)fftw_malloc(sizeof ( double) * m_warpings.nsamps_torig); 866 LPF_time = (double*)fftw_malloc(sizeof ( double) * m_warpings.nsamps_torig);
867 //memset((char*)LPF_time, 0, m_warpings.nsamps_torig * sizeof(double)); 867 //memset((char*)LPF_time, 0, m_warpings.nsamps_torig * sizeof(double));
868 // sustituyo el memset por un for: 868 // sustituyo el memset por un for:
869 for (size_t i = 0; i < m_warpings.nsamps_torig; i++) { 869 for (int i = 0; i < m_warpings.nsamps_torig; i++) {
870 LPF_time[i] = 0.0; 870 LPF_time[i] = 0.0;
871 } 871 }
872 #ifdef DEBUG 872 #ifdef DEBUG
873 printf(" Corrio primer memset...\n"); 873 printf(" Corrio primer memset...\n");
874 #endif 874 #endif
875 LPF_frequency = (fftw_complex*)fftw_malloc(sizeof ( fftw_complex) * (m_warpings.nsamps_torig/2 + 1)); //tamaño de la fft cuando la entrada es real 875 LPF_frequency = (fftw_complex*)fftw_malloc(sizeof ( fftw_complex) * (m_warpings.nsamps_torig/2 + 1)); //tamaño de la fft cuando la entrada es real
876 //memset((char*)LPF_frequency, 0, sizeof(fftw_complex) * (m_warpings.nsamps_torig/2 + 1)); 876 //memset((char*)LPF_frequency, 0, sizeof(fftw_complex) * (m_warpings.nsamps_torig/2 + 1));
877 // sustituyo el memset por un for: 877 // sustituyo el memset por un for:
878 for (size_t i = 0; i < (m_warpings.nsamps_torig/2 + 1); i++) { 878 for (int i = 0; i < (m_warpings.nsamps_torig/2 + 1); i++) {
879 LPF_frequency[i][0] = 0.0; 879 LPF_frequency[i][0] = 0.0;
880 LPF_frequency[i][1] = 0.0; 880 LPF_frequency[i][1] = 0.0;
881 } 881 }
882 // for (int i=0; i<(m_blockSize/2+1); i++) { 882 // for (int i=0; i<(m_blockSize/2+1); i++) {
883 // LPF_frequency[i] = new fftw_complex; 883 // LPF_frequency[i] = new fftw_complex;
884 // } 884 // }
885 plan_forward_LPF = fftw_plan_dft_r2c_1d(m_blockSize, LPF_time, LPF_frequency, FFTW_ESTIMATE); 885 plan_forward_LPF = fftw_plan_dft_r2c_1d(m_blockSize, LPF_time, LPF_frequency, FFTW_ESTIMATE);
886 plan_backward_LPF = fftw_plan_dft_c2r_1d(m_warpings.nsamps_torig, LPF_frequency, LPF_time, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT); 886 plan_backward_LPF = fftw_plan_dft_c2r_1d(m_warpings.nsamps_torig, LPF_frequency, LPF_time, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
887 887
888 size_t winWidth = 11; 888 int winWidth = 11;
889 double *lp_hanningWindow = new double[winWidth]; 889 double *lp_hanningWindow = new double[winWidth];
890 double accum=0; 890 double accum=0;
891 for (size_t i = 0; i < winWidth; i++) { 891 for (int i = 0; i < winWidth; i++) {
892 lp_hanningWindow[i]=0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)winWidth+1.0))); 892 lp_hanningWindow[i]=0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)winWidth+1.0)));
893 accum+=lp_hanningWindow[i]; 893 accum+=lp_hanningWindow[i];
894 894
895 } 895 }
896 for (size_t i = 0; i < winWidth; i++) { //window normalization 896 for (int i = 0; i < winWidth; i++) { //window normalization
897 lp_hanningWindow[i]=lp_hanningWindow[i]/accum; 897 lp_hanningWindow[i]=lp_hanningWindow[i]/accum;
898 } 898 }
899 for (size_t i = 0; i < m_blockSize/2+1; i++) { 899 for (int i = 0; i < m_blockSize/2+1; i++) {
900 //if (((i-(winWidth-1)/2)<0)||(i+(winWidth-1))/2>m_blockSize/2-1) {//consideramos winWidth impar, si la ventana sale del arreglo se rellena con el valor origianl 900 //if (((i-(winWidth-1)/2)<0)||(i+(winWidth-1))/2>m_blockSize/2-1) {//consideramos winWidth impar, si la ventana sale del arreglo se rellena con el valor origianl
901 if ( (i > (i_max + (winWidth-1)/2)) || (i <= (i_max - (winWidth-1)/2)) ) { 901 if ( (i > (i_max + (winWidth-1)/2)) || (i <= (i_max - (winWidth-1)/2)) ) {
902 mp_LPFWindow[i]=lp_LPFWindow_aux[i]; 902 mp_LPFWindow[i]=lp_LPFWindow_aux[i];
903 } else { 903 } else {
904 accum=0; 904 accum=0;
905 for (size_t j = -((winWidth-1)/2); j <= (winWidth-1)/2; j++) { 905 for (int j = -((winWidth-1)/2); j <= (winWidth-1)/2; j++) {
906 accum+=lp_LPFWindow_aux[i-j]*lp_hanningWindow[j+(winWidth-1)/2]; 906 accum+=lp_LPFWindow_aux[i-j]*lp_hanningWindow[j+(winWidth-1)/2];
907 } 907 }
908 mp_LPFWindow[i]=accum; 908 mp_LPFWindow[i]=accum;
909 } 909 }
910 } 910 }
913 delete[] lp_hanningWindow; 913 delete[] lp_hanningWindow;
914 } 914 }
915 915
916 void FChTransformF0gram::apply_LPF() { 916 void FChTransformF0gram::apply_LPF() {
917 fftw_execute(plan_forward_LPF); 917 fftw_execute(plan_forward_LPF);
918 for (size_t i = 0; i < m_blockSize/2+1; i++) { 918 for (int i = 0; i < m_blockSize/2+1; i++) {
919 LPF_frequency[i][0]*=mp_LPFWindow[i]; 919 LPF_frequency[i][0]*=mp_LPFWindow[i];
920 LPF_frequency[i][1]*=mp_LPFWindow[i]; 920 LPF_frequency[i][1]*=mp_LPFWindow[i];
921 } 921 }
922 fftw_execute(plan_backward_LPF); 922 fftw_execute(plan_backward_LPF);
923 923
1009 printf(" m_warpings.nsamps_torig = %d.\n",m_warpings.nsamps_torig); 1009 printf(" m_warpings.nsamps_torig = %d.\n",m_warpings.nsamps_torig);
1010 printf(" m_warp_params.num_warps = %d.\n",m_warp_params.num_warps); 1010 printf(" m_warp_params.num_warps = %d.\n",m_warp_params.num_warps);
1011 printf(" m_glogs_harmonic_count = %d.\n",m_glogs_harmonic_count); 1011 printf(" m_glogs_harmonic_count = %d.\n",m_glogs_harmonic_count);
1012 #endif 1012 #endif
1013 1013
1014 // size_t n = m_nfft/2 + 1; 1014 // int n = m_nfft/2 + 1;
1015 // double *tbuf = in_window; 1015 // double *tbuf = in_window;
1016 1016
1017 for (size_t i = 0; i < m_blockSize; i++) { 1017 for (int i = 0; i < m_blockSize; i++) {
1018 LPF_time[i] = (double)(inputBuffers[0][i]) * m_timeWindow[i]; 1018 LPF_time[i] = (double)(inputBuffers[0][i]) * m_timeWindow[i];
1019 } 1019 }
1020 1020
1021 // #ifdef DEBUG 1021 // #ifdef DEBUG
1022 // printf(" HASTA ACÁ ANDA!!!\n"); 1022 // printf(" HASTA ACÁ ANDA!!!\n");
1030 feature.hasTimestamp = false; 1030 feature.hasTimestamp = false;
1031 1031
1032 1032
1033 /* Solo a modo de prueba, voy a poner la salida del filtrado en «in» y 1033 /* Solo a modo de prueba, voy a poner la salida del filtrado en «in» y
1034 voy a mostrar la FFT de eso, para ver el efecto del filtrado. */ 1034 voy a mostrar la FFT de eso, para ver el efecto del filtrado. */
1035 // for (size_t i = 0; i < m_nfft; i++) { 1035 // for (int i = 0; i < m_nfft; i++) {
1036 // in[i][0] = tbuf[i]; 1036 // in[i][0] = tbuf[i];
1037 // in[i][1] = 0; 1037 // in[i][1] = 0;
1038 // } 1038 // }
1039 // fftw_execute(planFFT); 1039 // fftw_execute(planFFT);
1040 // double real, imag; 1040 // double real, imag;
1041 // for (size_t i=0; i<n; ++i) { // preincremento?? ver version de nacho 1041 // for (int i=0; i<n; ++i) { // preincremento?? ver version de nacho
1042 // real = out[i][0]; 1042 // real = out[i][0];
1043 // imag = out[i][1]; 1043 // imag = out[i][1];
1044 // feature.values.push_back(real*real + imag*imag); 1044 // feature.values.push_back(real*real + imag*imag);
1045 // } 1045 // }
1046 // fs[0].push_back(feature); 1046 // fs[0].push_back(feature);
1047 1047
1048 // float real; 1048 // float real;
1049 // float imag; 1049 // float imag;
1050 // for (size_t i=0; i<m_blockSize/2+1; i++) { 1050 // for (int i=0; i<m_blockSize/2+1; i++) {
1051 // real = (float)(LPF_frequency[i][0]); 1051 // real = (float)(LPF_frequency[i][0]);
1052 // imag = (float)(LPF_frequency[i][1]); 1052 // imag = (float)(LPF_frequency[i][1]);
1053 // feature.values.push_back(real*real+imag*imag); 1053 // feature.values.push_back(real*real+imag*imag);
1054 // //feature.values.push_back((float)(mp_LPFWindow[i])); 1054 // //feature.values.push_back((float)(mp_LPFWindow[i]));
1055 // } 1055 // }
1056 1056
1057 // ---------------------------------------------------------------------------------------------- 1057 // ----------------------------------------------------------------------------------------------
1058 // Hanning window & FFT for all warp directions 1058 // Hanning window & FFT for all warp directions
1059 1059
1060 double max_glogs = -DBL_MAX; 1060 double max_glogs = -DBL_MAX;
1061 size_t ind_max_glogs = 0; 1061 int ind_max_glogs = 0;
1062 1062
1063 for (size_t i_warp = 0; i_warp < m_warp_params.num_warps; i_warp++) { 1063 for (int i_warp = 0; i_warp < m_warp_params.num_warps; i_warp++) {
1064 // Interpolate 1064 // Interpolate
1065 interp1q(LPF_time, (m_warpings.pos_int) + i_warp*m_warp_params.nsamps_twarp, m_warpings.pos_frac + i_warp*m_warp_params.nsamps_twarp, x_warping, m_warp_params.nsamps_twarp); 1065 interp1q(LPF_time, (m_warpings.pos_int) + i_warp*m_warp_params.nsamps_twarp, m_warpings.pos_frac + i_warp*m_warp_params.nsamps_twarp, x_warping, m_warp_params.nsamps_twarp);
1066 1066
1067 // Apply window 1067 // Apply window
1068 for (size_t i = 0; i < m_warp_params.nsamps_twarp; i++) { 1068 for (int i = 0; i < m_warp_params.nsamps_twarp; i++) {
1069 x_warping[i] *= mp_HanningWindow[i]; 1069 x_warping[i] *= mp_HanningWindow[i];
1070 } 1070 }
1071 1071
1072 // Transform 1072 // Transform
1073 fftw_execute(plan_forward_xwarping); 1073 fftw_execute(plan_forward_xwarping);
1074 1074
1075 // Copy result 1075 // Copy result
1076 //memcpy(m_absFanChirpTransform + i_warp*(m_warp_params.nsamps_twarp/2+1), m_auxFanChirpTransform, (m_warp_params.nsamps_twarp/2+1)*sizeof(fftw_complex)); asi como esta no funciona 1076 //memcpy(m_absFanChirpTransform + i_warp*(m_warp_params.nsamps_twarp/2+1), m_auxFanChirpTransform, (m_warp_params.nsamps_twarp/2+1)*sizeof(fftw_complex)); asi como esta no funciona
1077 double *aux_abs_fcht = m_absFanChirpTransform + i_warp*(m_warp_params.nsamps_twarp/2+1); 1077 double *aux_abs_fcht = m_absFanChirpTransform + i_warp*(m_warp_params.nsamps_twarp/2+1);
1078 for (size_t i = 0; i < (m_warp_params.nsamps_twarp/2+1); i++) { 1078 for (int i = 0; i < (m_warp_params.nsamps_twarp/2+1); i++) {
1079 aux_abs_fcht[i] = log10(1.0 + 10.0*sqrt(m_auxFanChirpTransform[i][0]*m_auxFanChirpTransform[i][0]+m_auxFanChirpTransform[i][1]*m_auxFanChirpTransform[i][1])); 1079 aux_abs_fcht[i] = log10(1.0 + 10.0*sqrt(m_auxFanChirpTransform[i][0]*m_auxFanChirpTransform[i][0]+m_auxFanChirpTransform[i][1]*m_auxFanChirpTransform[i][1]));
1080 // smoothing high frequency values 1080 // smoothing high frequency values
1081 //aux_abs_fcht[i] *= m_glogs_hf_smoothing_window[i]; 1081 //aux_abs_fcht[i] *= m_glogs_hf_smoothing_window[i];
1082 } 1082 }
1083 1083
1084 // ----------------------------------------------------------------------------------------- 1084 // -----------------------------------------------------------------------------------------
1085 // GLogS 1085 // GLogS
1086 interp1q(aux_abs_fcht, m_glogs_posint, m_glogs_posfrac, m_glogs_interp, m_glogs_harmonic_count); 1086 interp1q(aux_abs_fcht, m_glogs_posint, m_glogs_posfrac, m_glogs_interp, m_glogs_harmonic_count);
1087 size_t glogs_ind = 0; 1087 int glogs_ind = 0;
1088 for (size_t i = 0; i < m_glogs_num_f0s; i++) { 1088 for (int i = 0; i < m_glogs_num_f0s; i++) {
1089 double glogs_accum = 0; 1089 double glogs_accum = 0;
1090 for (size_t j = 1; j <= m_glogs_n[i]; j++) { 1090 for (int j = 1; j <= m_glogs_n[i]; j++) {
1091 glogs_accum += m_glogs_interp[glogs_ind++]; 1091 glogs_accum += m_glogs_interp[glogs_ind++];
1092 } 1092 }
1093 m_glogs[i + i_warp*m_glogs_num_f0s] = glogs_accum/(double)m_glogs_n[i]; 1093 m_glogs[i + i_warp*m_glogs_num_f0s] = glogs_accum/(double)m_glogs_n[i];
1094 } 1094 }
1095 1095
1096 // Sub/super harmonic correction 1096 // Sub/super harmonic correction
1097 interp1q(m_glogs + i_warp*m_glogs_num_f0s, m_glogs_third_harmonic_posint, m_glogs_third_harmonic_posfrac, m_glogs_third_harmonic, (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct); 1097 interp1q(m_glogs + i_warp*m_glogs_num_f0s, m_glogs_third_harmonic_posint, m_glogs_third_harmonic_posfrac, m_glogs_third_harmonic, (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct);
1098 interp1q(m_glogs + i_warp*m_glogs_num_f0s, m_glogs_fifth_harmonic_posint, m_glogs_fifth_harmonic_posfrac, m_glogs_fifth_harmonic, (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct); 1098 interp1q(m_glogs + i_warp*m_glogs_num_f0s, m_glogs_fifth_harmonic_posint, m_glogs_fifth_harmonic_posfrac, m_glogs_fifth_harmonic, (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct);
1099 for (size_t i = m_glogs_num_f0s-1; i >= m_glogs_init_f0s; i--) { 1099 for (int i = m_glogs_num_f0s-1; i >= m_glogs_init_f0s; i--) {
1100 m_glogs[i + i_warp*m_glogs_num_f0s] -= MAX(MAX(m_glogs[i-m_f0_params.num_f0s_per_oct + i_warp*m_glogs_num_f0s],m_glogs_third_harmonic[i-m_glogs_init_f0s]),m_glogs_fifth_harmonic[i-m_glogs_init_f0s]); 1100 m_glogs[i + i_warp*m_glogs_num_f0s] -= MAX(MAX(m_glogs[i-m_f0_params.num_f0s_per_oct + i_warp*m_glogs_num_f0s],m_glogs_third_harmonic[i-m_glogs_init_f0s]),m_glogs_fifth_harmonic[i-m_glogs_init_f0s]);
1101 //m_glogs[i] -= MAX(m_glogs[i-m_f0_params.num_f0s_per_oct],m_glogs_third_harmonic[i-m_glogs_init_f0s]); 1101 //m_glogs[i] -= MAX(m_glogs[i-m_f0_params.num_f0s_per_oct],m_glogs_third_harmonic[i-m_glogs_init_f0s]);
1102 } 1102 }
1103 for (size_t i = m_glogs_init_f0s; i < m_glogs_num_f0s-m_f0_params.num_f0s_per_oct; i++) { 1103 for (int i = m_glogs_init_f0s; i < m_glogs_num_f0s-m_f0_params.num_f0s_per_oct; i++) {
1104 m_glogs[i + i_warp*m_glogs_num_f0s] -= 0.3*m_glogs[i+m_f0_params.num_f0s_per_oct + i_warp*m_glogs_num_f0s]; 1104 m_glogs[i + i_warp*m_glogs_num_f0s] -= 0.3*m_glogs[i+m_f0_params.num_f0s_per_oct + i_warp*m_glogs_num_f0s];
1105 // Median, sigma $ weights correction 1105 // Median, sigma $ weights correction
1106 m_glogs[i + i_warp*m_glogs_num_f0s] = (m_glogs[i + i_warp*m_glogs_num_f0s]-m_glogs_median_correction[i-m_glogs_init_f0s])*m_glogs_sigma_correction[i-m_glogs_init_f0s]*m_glogs_f0_preference_weights[i-m_glogs_init_f0s]; 1106 m_glogs[i + i_warp*m_glogs_num_f0s] = (m_glogs[i + i_warp*m_glogs_num_f0s]-m_glogs_median_correction[i-m_glogs_init_f0s])*m_glogs_sigma_correction[i-m_glogs_init_f0s]*m_glogs_f0_preference_weights[i-m_glogs_init_f0s];
1107 } 1107 }
1108 1108
1109 // Look for maximum value to determine best direction 1109 // Look for maximum value to determine best direction
1110 for (size_t i = m_glogs_init_f0s; i < m_glogs_num_f0s-m_f0_params.num_f0s_per_oct; i++) { 1110 for (int i = m_glogs_init_f0s; i < m_glogs_num_f0s-m_f0_params.num_f0s_per_oct; i++) {
1111 if (m_glogs[i + i_warp*m_glogs_num_f0s] > max_glogs) { 1111 if (m_glogs[i + i_warp*m_glogs_num_f0s] > max_glogs) {
1112 max_glogs = m_glogs[i + i_warp*m_glogs_num_f0s]; 1112 max_glogs = m_glogs[i + i_warp*m_glogs_num_f0s];
1113 ind_max_glogs = i_warp; 1113 ind_max_glogs = i_warp;
1114 } 1114 }
1115 } 1115 }
1116 } 1116 }
1117 1117
1118 // ---------------------------------------------------------------------------------------------- 1118 // ----------------------------------------------------------------------------------------------
1119 1119
1120 for (size_t i=m_glogs_init_f0s; i< m_glogs_num_f0s - m_f0_params.num_f0s_per_oct; i++) { 1120 for (int i=m_glogs_init_f0s; i< m_glogs_num_f0s - m_f0_params.num_f0s_per_oct; i++) {
1121 //for (size_t i=0; i<(m_warp_params.nsamps_twarp/2+1); i++) { 1121 //for (int i=0; i<(m_warp_params.nsamps_twarp/2+1); i++) {
1122 //feature.values.push_back((float)(m_warpings.pos_int[i])+ (float)(m_warpings.pos_frac[i])); 1122 //feature.values.push_back((float)(m_warpings.pos_int[i])+ (float)(m_warpings.pos_frac[i]));
1123 //feature.values.push_back((float)(phi[i]*100000.0)); 1123 //feature.values.push_back((float)(phi[i]*100000.0));
1124 //feature.values.push_back((float)(t_orig[i])); 1124 //feature.values.push_back((float)(t_orig[i]));
1125 //feature.values.push_back((float)(pos1[i])); 1125 //feature.values.push_back((float)(pos1[i]));
1126 //feature.values.push_back((float)x_warping[i]); 1126 //feature.values.push_back((float)x_warping[i]);
1127 //feature.values.push_back(m_absFanChirpTransform[i + ind_max_glogs*(m_warp_params.nsamps_twarp/2+1)]); 1127 //feature.values.push_back(m_absFanChirpTransform[i + ind_max_glogs*(m_warp_params.nsamps_twarp/2+1)]);
1128 //feature.values.push_back((float)m_glogs[i+(long)ind_max_glogs*(long)m_glogs_num_f0s]); 1128 //feature.values.push_back((float)m_glogs[i+(long)ind_max_glogs*(long)m_glogs_num_f0s]);
1129 switch (m_f0gram_mode) { 1129 switch (m_f0gram_mode) {
1130 case 1: 1130 case 1:
1131 max_glogs = -DBL_MAX; 1131 max_glogs = -DBL_MAX;
1132 for (size_t i_warp = 0; i_warp < m_warp_params.num_warps; i_warp++) { 1132 for (int i_warp = 0; i_warp < m_warp_params.num_warps; i_warp++) {
1133 if (m_glogs[i + i_warp*m_glogs_num_f0s] > max_glogs) { 1133 if (m_glogs[i + i_warp*m_glogs_num_f0s] > max_glogs) {
1134 max_glogs = m_glogs[i + i_warp*m_glogs_num_f0s]; 1134 max_glogs = m_glogs[i + i_warp*m_glogs_num_f0s];
1135 ind_max_glogs = i_warp; 1135 ind_max_glogs = i_warp;
1136 } 1136 }
1137 } 1137 }
1138 feature.values.push_back((float)max_glogs); 1138 feature.values.push_back((float)max_glogs);
1139 break; 1139 break;
1140 case 0: 1140 case 0:
1141 feature.values.push_back((float)m_glogs[i+(size_t)ind_max_glogs*(size_t)m_glogs_num_f0s]); 1141 feature.values.push_back((float)m_glogs[i+(int)ind_max_glogs*(int)m_glogs_num_f0s]);
1142 break; 1142 break;
1143 } 1143 }
1144 //feature.values.push_back((float)m_glogs_hf_smoothing_window[i]); 1144 //feature.values.push_back((float)m_glogs_hf_smoothing_window[i]);
1145 } 1145 }
1146 1146
1164 } 1164 }
1165 1165
1166 void 1166 void
1167 FChTransformF0gram::design_time_window() { 1167 FChTransformF0gram::design_time_window() {
1168 1168
1169 size_t transitionWidth = (size_t)m_blockSize/128 + 1;; 1169 int transitionWidth = (int)m_blockSize/128 + 1;;
1170 m_timeWindow = new double[m_blockSize]; 1170 m_timeWindow = new double[m_blockSize];
1171 double *lp_transitionWindow = new double[transitionWidth]; 1171 double *lp_transitionWindow = new double[transitionWidth];
1172 1172
1173 //memset(m_timeWindow, 1.0, m_blockSize); 1173 //memset(m_timeWindow, 1.0, m_blockSize);
1174 for (size_t i = 0; i < m_blockSize; i++) { 1174 for (int i = 0; i < m_blockSize; i++) {
1175 m_timeWindow[i] = 1.0; 1175 m_timeWindow[i] = 1.0;
1176 } 1176 }
1177 1177
1178 for (size_t i = 0; i < transitionWidth; i++) { 1178 for (int i = 0; i < transitionWidth; i++) {
1179 lp_transitionWindow[i]=0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)transitionWidth+1.0))); 1179 lp_transitionWindow[i]=0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)transitionWidth+1.0)));
1180 } 1180 }
1181 1181
1182 for (size_t i = 0; i < transitionWidth/2; i++) { 1182 for (int i = 0; i < transitionWidth/2; i++) {
1183 m_timeWindow[i] = lp_transitionWindow[i]; 1183 m_timeWindow[i] = lp_transitionWindow[i];
1184 m_timeWindow[m_blockSize-1-i] = lp_transitionWindow[transitionWidth-1-i]; 1184 m_timeWindow[m_blockSize-1-i] = lp_transitionWindow[transitionWidth-1-i];
1185 } 1185 }
1186 1186
1187 #ifdef DEBUG 1187 #ifdef DEBUG