Mercurial > hg > vamp-fanchirp
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 |
