#include #include #include #include #include #if TEST_USE_LINEAR_HARMONIC_ENVELOPE #include #endif #include #include #include #include #include #include #include "spp.h" static double harmonic_get_frequency(const void *harmonic) { return ((const struct narrowband_harmonic*)harmonic) -> frequency; } void spp_analyze(struct spp_frame * restrict spp_frames, positive long frame_count, const float * restrict samples) { // XXX: Handle this in other places #if TEST_SPP_WINDOW_TYPE == SPP_WINDOW_2048 static thread_local float window_function[DFT_HALF(SPP_WINDOW_SIZE)]; generate_window(window_function, SPP_WINDOW_SIZE, WINDOW_BLACKMANHARRIS); #elif TEST_SPP_WINDOW_TYPE == SPP_WINDOW_PREPAD static thread_local float window_function[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; generate_window(window_function, SPP_PADDED_WINDOW_SIZE, WINDOW_BLACKMANHARRIS); #else #error #endif for (positive long frame_index = 0; frame_index < frame_count; frame_index++) { // XXX: Make this true instead #if TEST_SPP_WINDOW_TYPE == SPP_WINDOW_2048 stdft(samples + frame_index * HOP_SIZE, window_function, spp_frames[frame_index].spectrum, SPP_WINDOW_SIZE, SPP_WINDOW_ZEROPADDING_FACTOR, false); #elif TEST_SPP_WINDOW_TYPE == SPP_WINDOW_PREPAD stdft(samples + frame_index * HOP_SIZE, window_function, spp_frames[frame_index].spectrum, SPP_PADDED_WINDOW_SIZE, 1, false); #else #error #endif }; } #if TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SCALE void spp_synthesize(float * restrict buffer, positive long sample_count, const struct spp_frame * restrict spp_frames, positive long frame_count, float pitch_scale) { #elif TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SET // NOTE: In hertz, not cents void spp_synthesize(float * restrict buffer, positive long sample_count, const struct spp_frame * restrict spp_frames, positive long frame_count, float pitch) { #else #error #endif static thread_local float window_function[DFT_HALF(SPP_WINDOW_SIZE)]; generate_window(window_function, SPP_WINDOW_SIZE, WINDOW_BLACKMANHARRIS); for (positive long itr = 0; itr < sample_count; itr++) { buffer[itr] = 0.0f; }; const float analysis_frequency_step = (float)SAMPLE_RATE / (float)SPP_PADDED_WINDOW_SIZE; const float synthesis_frequency_step = (float)SAMPLE_RATE / (float)SPP_WINDOW_SIZE; float fundamental_phase = 0.0f; // XXX: Finish implementing this #if TEST_SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS == SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS_YES static thread_local positive int prev_peak_sources[MAX_HARMONICS]; static thread_local float harmonic_phase_correction[MAX_HARMONICS]; // XXX: Also consider a potential bug in the original where prev_peak_sources is not reset #if TEST_SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE == SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE_INITIALIZATION_ONLY for (positive int itr = 0; itr < MAX_HARMONICS; itr++) { harmonic_phase_correction[itr] = 0.0f; }; #elif TEST_SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE == SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE_NEW_HARMONICS positive int prev_new_harmonic_count = 0; #else #error #endif #elif TEST_SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS == SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS_NO // Intentionally blank #else #error #endif #if TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_2C_SPECTRAL_GAP_FILLING positive int spectral_gap_start = 0; #endif for (positive long frame_index = 0; frame_index < frame_count; frame_index++) { static thread_local complex float spectrum[SPP_WINDOW_SIZE]; static thread_local float window[SPP_WINDOW_SIZE]; for (positive int itr = 0; itr < DFT_HALF(SPP_WINDOW_SIZE); itr++) { spectrum[itr] = 0.0f; }; #if TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_SPLINE static thread_local spline_scratchspace_item_t spectrum_spline_scratchspace[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; static thread_local float spline_positions[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; static thread_local float spline_values[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; for (positive int itr = 0; itr < DFT_HALF(SPP_PADDED_WINDOW_SIZE); itr++) { spline_positions[itr] = analysis_frequency_step * itr; }; for (positive int itr = 0; itr < DFT_HALF(SPP_PADDED_WINDOW_SIZE); itr++) { #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR // XXX: DC/Nyquist #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC spline_values[itr] = logf(AMPLITUDE(spp_frames[frame_index].spectrum[itr])); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR spline_values[itr] = AMPLITUDE(spp_frames[frame_index].spectrum[itr]); #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR spline_values[itr] = creal(spp_frames[frame_index].spectrum[itr]); #else #error #endif }; #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR static thread_local struct spline_segment spectrum_spline_amplitude[SPLINE_SIZE(DFT_HALF(SPP_PADDED_WINDOW_SIZE))]; create_spline(spline_positions, spline_values, spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), SPLINE_NATURAL, spectrum_spline_scratchspace); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR static thread_local struct spline_segment spectrum_spline_real_part[SPLINE_SIZE(DFT_HALF(SPP_PADDED_WINDOW_SIZE))]; create_spline(spline_positions, spline_values, spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), SPLINE_NATURAL, spectrum_spline_scratchspace); #else #error #endif for (positive long itr = 0; itr < DFT_HALF(SPP_PADDED_WINDOW_SIZE); itr++) { #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE3 == SPP_INTERPOLATION3_PHASE_UNWRAPPED if (itr > 0) { spline_values[itr] = spline_values[itr - 1] + princarg(PHASE(spp_frames[frame_index].spectrum[itr]) - spline_values[itr - 1]); } else { spline_values[itr] = PHASE(spp_frames[frame_index].spectrum[itr]); }; #elif TEST_SPP_INTERPOLATION_TYPE3 == SPP_INTERPOLATION3_PHASE_NOT_UNWRAPPED spline_values[itr] = PHASE(spp_frames[frame_index].spectrum[itr]); #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR spline_values[itr] = cimag(spp_frames[frame_index].spectrum[itr]); #else #error #endif }; #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR static thread_local struct spline_segment spectrum_spline_phase[SPLINE_SIZE(DFT_HALF(SPP_PADDED_WINDOW_SIZE))]; create_spline(spline_positions, spline_values, spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), SPLINE_NATURAL, spectrum_spline_scratchspace); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR static thread_local struct spline_segment spectrum_spline_imaginary_part[SPLINE_SIZE(DFT_HALF(SPP_PADDED_WINDOW_SIZE))]; create_spline(spline_positions, spline_values, spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), SPLINE_NATURAL, spectrum_spline_scratchspace); #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_LINEAR #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR static thread_local struct line_segment spectrum_spline_amplitude[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; static thread_local struct line_segment spectrum_spline_phase[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR static thread_local struct line_segment spectrum_spline_real_part[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; static thread_local struct line_segment spectrum_spline_imaginary_part[DFT_HALF(SPP_PADDED_WINDOW_SIZE)]; #else #error #endif for (positive int itr = 0; itr < DFT_HALF(SPP_PADDED_WINDOW_SIZE); itr++) { #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR spectrum_spline_amplitude[itr].position = analysis_frequency_step * itr; spectrum_spline_phase[itr].position = analysis_frequency_step * itr; #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR spectrum_spline_real_part[itr].position = analysis_frequency_step * itr; spectrum_spline_imaginary_part[itr].position = analysis_frequency_step * itr; #else #error #endif #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR // XXX: DC/Nyquist #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC spectrum_spline_amplitude[itr].value = logf(AMPLITUDE(spp_frames[frame_index].spectrum[itr])); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR spectrum_spline_amplitude[itr].value = AMPLITUDE(spp_frames[frame_index].spectrum[itr]); #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR spectrum_spline_real_part[itr].value = creal(spp_frames[frame_index].spectrum[itr]); #else #error #endif #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE3 == SPP_INTERPOLATION3_PHASE_UNWRAPPED if (itr > 0) { spectrum_spline_phase[itr - 1].value = spectrum_spline_phase[itr].value + princarg(PHASE(spp_frames[frame_index].spectrum[itr]) - spectrum_spline_phase[itr - 1].value); } else { spectrum_spline_phase[itr].value = PHASE(spp_frames[frame_index].spectrum[itr]); }; #elif TEST_SPP_INTERPOLATION_TYPE3 == SPP_INTERPOLATION3_PHASE_NOT_UNWRAPPED spline_values[itr] = PHASE(spp_frames[frame_index].spectrum[itr]); #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR spectrum_spline_imaginary_part[itr].value = cimag(spp_frames[frame_index].spectrum[itr]); #else #error #endif }; #else #error #endif // NOTE: Assumes voiced #if TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SCALE float new_fundamental_frequency = spp_frames[frame_index].fundamental_frequency * pitch_scale; #elif TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SET float new_fundamental_frequency = pitch; #else #error #endif positive int new_harmonic_count = convert_float_int(floorf(NYQUIST_FREQUENCY / new_fundamental_frequency), 1, MAX_HARMONICS); #if TEST_SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS == SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS_YES && TEST_SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE == SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE_NEW_HARMONICS for (positive int itr = prev_new_harmonic_count; itr < new_harmonic_count; itr++) { harmonic_phase_correction[itr] = 0.0f; }; #endif // Strictly speaking, this is incorrect, but the original could have had a bug or implementation quirk like this. Besides, it requires almost no additional work to implement this variant. #if TEST_SPP_PHASE_PROPGATION_TIME == SPP_PHASE_PROPGATION_TIME_BEFORE #if TEST_SPP_PHASE_PROPGATION_METHOD == SPP_PHASE_PROPGATION_METHOD_CONSTANT fundamental_phase = princarg(fundamental_phase + new_fundamental_frequency * HOP_TIME * M_PIF * 2.0f); #elif TEST_SPP_PHASE_PROPGATION_METHOD == SPP_PHASE_PROPGATION_METHOD_LINEAR #if TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SCALE float next_new_fundamental_frequency; if (frame_index < frame_count - 1) { next_new_fundamental_frequency = spp_frames[frame_index + 1].fundamental_frequency * pitch_scale; } else { next_new_fundamental_frequency = spp_frames[frame_index].fundamental_frequency * pitch_scale; }; #elif TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SET float next_new_fundamental_frequency = pitch; #else #error #endif fundamental_phase = princarg(fundamental_phase + (new_fundamental_frequency + next_new_fundamental_frequency) * HOP_TIME * M_PIF); // XXX: Also consider the quadratic phase propagation technique from Bonada's 1997 Master's Thesis #else #error #endif #elif TEST_SPP_PHASE_PROPGATION_TIME == SPP_PHASE_PROPGATION_TIME_AFTER // Intentionally blank #else #error #endif // XXX: Consider ideal-frequency versions for both of these. Same at interpolation. #if TEST_SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_SPLINE static thread_local struct spline_segment harmonic_amplitude_envelope[SPLINE_SIZE(MAX_HARMONICS)]; static thread_local float harmonic_amplitude_envelope_positions[MAX_HARMONICS]; static thread_local float harmonic_amplitude_envelope_values[MAX_HARMONICS]; for (positive int itr = 0; itr < spp_frames[frame_index].harmonics.harmonic_count; itr++) { harmonic_amplitude_envelope_positions[itr] = spp_frames[frame_index].harmonics.harmonics[itr].frequency; harmonic_amplitude_envelope_values[itr] = logf(spp_frames[frame_index].harmonics.harmonics[itr].amplitude); }; spline_scratchspace_item_t harmonic_amplitude_envelope_scratchspace[MAX_HARMONICS]; create_spline(harmonic_amplitude_envelope_positions, harmonic_amplitude_envelope_values, harmonic_amplitude_envelope, spp_frames[frame_index].harmonics.harmonic_count, SPLINE_NATURAL, harmonic_amplitude_envelope_scratchspace); #elif TEST_SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_LINEAR static thread_local struct line_segment harmonic_amplitude_envelope[MAX_HARMONICS]; for (positive int itr = 0; itr < spp_frames[frame_index].harmonics.harmonic_count; itr++) { harmonic_amplitude_envelope[itr].position = spp_frames[frame_index].harmonics.harmonics[itr].frequency; harmonic_amplitude_envelope[itr].value = logf(spp_frames[frame_index].harmonics.harmonics[itr].amplitude); }; #else #error #endif #if TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_SPLINE static thread_local struct spline_segment harmonic_phase_envelope[SPLINE_SIZE(MAX_HARMONICS)]; static thread_local float harmonic_phase_envelope_positions[MAX_HARMONICS]; static thread_local float harmonic_phase_envelope_values[MAX_HARMONICS]; for (positive int itr = 0; itr < spp_frames[frame_index].harmonics.harmonic_count; itr++) { harmonic_phase_envelope_positions[itr] = spp_frames[frame_index].harmonics.harmonics[itr].frequency; harmonic_phase_envelope_values[itr] = spp_frames[frame_index].mfpa_envelope[itr]; }; spline_scratchspace_item_t harmonic_phase_envelope_scratchspace[MAX_HARMONICS]; create_spline(harmonic_phase_envelope_positions, harmonic_phase_envelope_values, harmonic_phase_envelope, spp_frames[frame_index].harmonics.harmonic_count, SPLINE_NATURAL, harmonic_phase_envelope_scratchspace); #elif TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_LINEAR static thread_local struct line_segment harmonic_phase_envelope[MAX_HARMONICS]; for (positive int itr = 0; itr < spp_frames[frame_index].harmonics.harmonic_count; itr++) { harmonic_phase_envelope[itr].position = spp_frames[frame_index].harmonics.harmonics[itr].frequency; harmonic_phase_envelope[itr].value = spp_frames[frame_index].mfpa_envelope[itr]; }; #elif TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_MAPPING_CLOSEST // Intentionally empty #elif TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_MAPPING_BIJECTIVE // Intentionally empty #else #error #endif #if TEST_SPP_PEAK_PICKING_METHOD == SPP_PEAK_PICKING_BIJECTIVE #if TEST_SPP_PEAK_PICKING_BIJECTIVE_REPEAT_LAST == SPP_PEAK_PICKING_BIJECT_REPEAT_LAST_YES // XXX: Actually implement this for (positive int synthesis_harmonic_index = 0; synthesis_harmonic_index < new_harmonic_count; synthesis_harmonic_index++) { #elif TEST_SPP_PEAK_PICKING_BIJECTIVE_REPEAT_LAST == SPP_PEAK_PICKING_BIJECT_REPEAT_LAST_NO for (positive int synthesis_harmonic_index = 0; synthesis_harmonic_index < min(new_harmonic_count, spp_frames[frame_index].harmonics.harmonic_count); synthesis_harmonic_index++) { #else #error #endif #else for (positive int synthesis_harmonic_index = 0; synthesis_harmonic_index < new_harmonic_count; synthesis_harmonic_index++) { #endif #if TEST_SPP_SYNTHESIS_HARMONIC_DISTRIBUTION == SPP_SYNTHESIS_HARMONIC_DISTRIBUTION_IDEAL float synthesis_harmonic_frequency = new_fundamental_frequency * (synthesis_harmonic_index + 1); #elif TEST_SPP_SYNTHESIS_HARMONIC_DISTRIBUTION == SPP_SYNTHESIS_HARMONIC_DISTRIBUTION_SCALE_ANALYSIS static_assert(TEST_SPP_PEAK_PICKING_METHOD == SPP_PEAK_PICKING_BIJECTIVE); float synthesis_harmonic_frequency = spp_frames[frame_index].harmonics.harmonics[synthesis_harmonic_index] * (new_fundamental_frequency / spp_frames[frame_index].fundamental_frequency) * (float)(synthesis_harmonic_index + 1); #elif TEST_SPP_SYNTHESIS_HARMONIC_DISTRIBUTION == SPP_SYNTHESIS_HARMONIC_DISTRIBUTION_ADD_ANALYSIS static_assert(TEST_SPP_PEAK_PICKING_METHOD == SPP_PEAK_PICKING_BIJECTIVE); float synthesis_harmonic_frequency = spp_frames[frame_index].harmonics.harmonics[synthesis_harmonic_index] + (new_fundamental_frequency - spp_frames[frame_index].fundamental_frequency) * (float)(synthesis_harmonic_index + 1); #else #error #endif #if TEST_SPP_PEAK_PICKING_METHOD == SPP_PEAK_PICKING_BIJECTIVE positive int analysis_harmonic_index = synthesis_harmonic_index; #elif TEST_SPP_PEAK_PICKING_METHOD == SPP_PEAK_PICKING_CLOSEST // XXX: Further distinguish between without and with timbre scaling factor when that is reintroduced #if TEST_SPP_PEAK_PICKING_CLOSEST_TYPE == SPP_PEAK_PICKING_CLOSEST_ROUND positive int analysis_harmonic_index = (const struct narrowband_harmonic*)search_closest(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), synthesis_harmonic_frequency, harmonic_get_frequency) - spp_frames[frame_index].harmonics.harmonics; #elif TEST_SPP_PEAK_PICKING_CLOSEST_TYPE == SPP_PEAK_PICKING_CLOSEST_FLOOR positive int analysis_harmonic_index = (const struct narrowband_harmonic*)search_floor(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), synthesis_harmonic_frequency, harmonic_get_frequency) - spp_frames[frame_index].harmonics.harmonics; #elif TEST_SPP_PEAK_PICKING_CLOSEST_TYPE == SPP_PEAK_PICKING_CLOSEST_CEIL positive int analysis_harmonic_index = (const struct narrowband_harmonic*)search_ceil(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), synthesis_harmonic_frequency, harmonic_get_frequency) - spp_frames[frame_index].harmonics.harmonics; #else #error #endif // XXX: Also consider weighted frequency+(log)amplitude method? #elif TEST_SPP_PEAK_PICKING_METHOD == SPP_PEAK_PICKING_BONADA2008 // XXX: Hysteresis #if TEST_SPP_PEAK_PICKING_CLOSEST_TYPE == SPP_PEAK_PICKING_CLOSEST_ROUND positive int analysis_harmonic_index = (const struct narrowband_harmonic*)search_closest(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), synthesis_harmonic_frequency, harmonic_get_frequency) - spp_frames[frame_index].harmonics.harmonics; #elif TEST_SPP_PEAK_PICKING_CLOSEST_TYPE == SPP_PEAK_PICKING_CLOSEST_FLOOR positive int analysis_harmonic_index = (const struct narrowband_harmonic*)search_floor(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), synthesis_harmonic_frequency, harmonic_get_frequency) - spp_frames[frame_index].harmonics.harmonics; #elif TEST_SPP_PEAK_PICKING_CLOSEST_TYPE == SPP_PEAK_PICKING_CLOSEST_CEIL positive int analysis_harmonic_index = (const struct narrowband_harmonic*)search_ceil(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), synthesis_harmonic_frequency, harmonic_get_frequency) - spp_frames[frame_index].harmonics.harmonics; #else #error #endif if (analysis_harmonic_index > 0 && logdeltaf(spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index - 1].amplitude, spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].amplitude) > SPP_PEAK_PICKING_BONADA2008_VALLEY_THRESHOLD) { analysis_harmonic_index--; } else if (analysis_harmonic_index < spp_frames[frame_index].harmonics.harmonic_count - 1 && logdeltaf(spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index + 1].amplitude, spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].amplitude) > SPP_PEAK_PICKING_BONADA2008_VALLEY_THRESHOLD) { analysis_harmonic_index++; }; #else #error #endif #if TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_IDEAL float harmonic_region_center = spp_frames[frame_index].fundamental_frequency * (analysis_harmonic_index + 1); #elif TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_ANALYSIS float harmonic_region_center = spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index]; #else #error #endif // XXX: Consider different rounding modes here too #if TEST_SPP_REGION_DELIMITATION_METHOD == SPP_REGION_DELIMITATION_MIDDLE #if TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_IDEAL unused float harmonic_region_start = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 0.5f); unused float harmonic_region_end = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.5f); #elif TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_ANALYSIS unused float harmonic_region_start; if (analysis_harmonic_index > 0) { harmonic_region_start = (spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index] + spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index - 1]) / 2.0f; } else { harmonic_region_start = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 0.5f); }; unused float harmonic_region_end; if (analysis_harmonic_index < spp_frames[frame_index].harmonics.harmonic_index - 1) { harmonic_region_end = (spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index] + spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index + 1]) / 2.0f; } else { harmonic_region_end = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.5f); }; #else #error #endif #elif TEST_SPP_REGION_DELIMITATION_METHOD == SPP_REGION_DELIMITATION_LOWEST unused float harmonic_region_start; unused float harmonic_region_end; if (analysis_harmonic_index > 0) { #if TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_IDEAL positive int start_bin = convert_float_int(roundf((spp_frames[frame_index].fundamental_frequency * analysis_harmonic_index) / analysis_frequency_step), 1, SPP_PADDED_WINDOW_SIZE - 2); positive int end_bin = convert_float_int(roundf((spp_frames[frame_index].fundamental_frequency * (analysis_harmonic_index + 1)) / analysis_frequency_step), start_bin + 1, SPP_PADDED_WINDOW_SIZE - 2); #elif TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_ANALYSIS positive int start_bin = convert_float_int(roundf(spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index - 1].frequency / analysis_frequency_step), 1, SPP_PADDED_WINDOW_SIZE - 2); positive int end_bin = convert_float_int(roundf(spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].frequency / analysis_frequency_step), start_bin + 1, SPP_PADDED_WINDOW_SIZE - 2); #else #error #endif float lowest_bin_frequency; float lowest_bin_value = INFINITY; for (positive int itr = start_bin; itr < end_bin; itr++) { if (AMPLITUDE(spp_frames[frame_index].spectrum[itr]) < lowest_bin_value) { #if TEST_SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE == SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE_NONE lowest_bin_frequency = analysis_frequency_step * itr; lowest_bin_value = AMPLITUDE(spp_frames[frame_index].spectrum[itr]); #elif TEST_SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE == SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE_QUADRATIC float _logamp[3]; float *logamp = _logamp + 1; for (int offset = -1; offset <= +1; offset++) { logamp[offset] = logf(AMPLITUDE(spp_frames[frame_index].spectrum[itr + offset])); }; float new_value = expf(logamp[0] - powf(logamp[-1] - logamp[+1], 2.0f) / ((logamp[-1] - logamp[0] * 2.0f + logamp[+1]) * 8.0f)); if (new_value < lowest_bin_value) { lowest_bin_frequency = ((float)itr + (logamp[-1] - logamp[+1]) / ((logamp[-1] - logamp[0] * 2.0f + logamp[+1]) * 2.0f)) * analysis_frequency_step; lowest_bin_value = new_value; }; #else #error #endif }; }; harmonic_region_start = lowest_bin_frequency; } else { harmonic_region_start = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 0.5f); }; if (analysis_harmonic_index < spp_frames[frame_index].harmonics.harmonic_index - 1) { #if TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_IDEAL positive int start_bin = convert_float_int(roundf((spp_frames[frame_index].fundamental_frequency * (analysis_harmonic_index + 1)) / analysis_frequency_step), 1, SPP_PADDED_WINDOW_SIZE - 2); positive int end_bin = convert_float_int(roundf((spp_frames[frame_index].fundamental_frequency * (analysis_harmonic_index + 2)) / analysis_frequency_step), start_bin + 1, SPP_PADDED_WINDOW_SIZE - 2); #elif TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_ANALYSIS positive int start_bin = convert_float_int(roundf(spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].frequency / analysis_frequency_step), 1, SPP_PADDED_WINDOW_SIZE - 2); positive int end_bin = convert_float_int(roundf(spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index + 1].frequency / analysis_frequency_step), start_bin + 1, SPP_PADDED_WINDOW_SIZE - 2); #else #error #endif float lowest_bin_frequency; float lowest_bin_value = INFINITY; for (positive int itr = start_bin; itr < end_bin; itr++) { if (AMPLITUDE(spp_frames[frame_index].spectrum[itr]) < lowest_bin_value) { #if TEST_SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE == SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE_NONE lowest_bin_frequency = analysis_frequency_step * itr; lowest_bin_value = AMPLITUDE(spp_frames[frame_index].spectrum[itr]); #elif TEST_SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE == SPP_REGION_DELIMITATION_LOWEST_INTERPOLATION_MODE_QUADRATIC float _logamp[3]; float *logamp = _logamp + 1; for (int offset = -1; offset <= +1; offset++) { logamp[offset] = logf(AMPLITUDE(spp_frames[frame_index].spectrum[itr + offset])); }; float new_value = expf(logamp[0] - powf(logamp[-1] - logamp[+1], 2.0f) / ((logamp[-1] - logamp[0] * 2.0f + logamp[+1]) * 8.0f)); if (new_value < lowest_bin_value) { lowest_bin_frequency = ((float)itr + (logamp[-1] - logamp[+1]) / ((logamp[-1] - logamp[0] * 2.0f + logamp[+1]) * 2.0f)) * analysis_frequency_step; lowest_bin_value = new_value; }; #else #error #endif }; }; harmonic_region_end = lowest_bin_frequency; } else { harmonic_region_end = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.5f); }; #elif TEST_SPP_REGION_DELIMITATION_METHOD == SPP_REGION_DELIMITATION_MIDDLE_SCALED #if TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_IDEAL unused float harmonic_region_start = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.0f - SPP_REGION_DELIMITATION_MIDDLE_SCALING_FACTOR / 2.0f); unused float harmonic_region_end = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.0f + SPP_REGION_DELIMITATION_MIDDLE_SCALING_FACTOR / 2.0f); #elif TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_ANALYSIS unused float harmonic_region_start; if (analysis_harmonic_index > 0) { harmonic_region_start = spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index] - (spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index] - spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index - 1]) * (SPP_REGION_DELIMITATION_MIDDLE_SCALING_FACTOR / 2.0f); } else { harmonic_region_start = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.0f - SPP_REGION_DELIMITATION_MIDDLE_SCALING_FACTOR / 2.0f); }; unused float harmonic_region_end; if (analysis_harmonic_index < spp_frames[frame_index].harmonics.harmonic_index - 1) { harmonic_region_end = spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index] + (spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index + 1] - spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index]) * (SPP_REGION_DELIMITATION_MIDDLE_SCALING_FACTOR / 2.0f); } else { harmonic_region_end = spp_frames[frame_index].fundamental_frequency * ((float)analysis_harmonic_index + 1.0f + SPP_REGION_DELIMITATION_MIDDLE_SCALING_FACTOR / 2.0f); }; #else #error #endif #else #error #endif #if TEST_SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_SPLINE float amplitude_scale = expf(interpolate_spline(harmonic_amplitude_envelope, spp_frames[frame_index].harmonics.harmonic_count, synthesis_harmonic_frequency)) / spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].amplitude; #elif TEST_SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_AMPLITUDE_ENVELOPE_INTERPOLATION_LINEAR float amplitude_scale = expf(interpolate_line_segments(harmonic_amplitude_envelope, spp_frames[frame_index].harmonics.harmonic_count, synthesis_harmonic_frequency)) / spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].amplitude; #else #error #endif #if TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_SPLINE float phase_rotation = interpolate_spline(harmonic_phase_envelope, spp_frames[frame_index].harmonics.harmonic_count, synthesis_harmonic_frequency); #elif TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_LINEAR float phase_rotation = interpolate_line_segments(harmonic_phase_envelope, spp_frames[frame_index].harmonics.harmonic_count, synthesis_harmonic_frequency); #elif TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_MAPPING_CLOSEST float phase_rotation = spp_frames[frame_index].mfpa_envelope[analysis_harmonic_index]; #elif TEST_SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_TYPE == SPP_HARMONIC_PHASE_ENVELOPE_INTERPOLATION_MAPPING_BIJECTIVE float phase_rotation = spp_frames[frame_index].mfpa_envelope[min(itr, spp_frames[frame_index].harmonics.harmonic_count - 1)]; #else #error #endif phase_rotation -= spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index].phase; phase_rotation += (fundamental_phase - spp_frames[frame_index].mfpa_envelope[0]) * (float)(synthesis_harmonic_index + 1); #if TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_1B_ADAPTIVE_REGION_SEGMENTATION float harmonic_region_start; if (analysis_harmonic_index > 0) { harmonic_region_start = max(harmonic_region_center - new_fundamental_frequency / 2.0f, spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index - 1] + SPP_REGION_RENDERING_STRATEGY_1B_ADAPTIVE_REGION_SEGMENTATION_MINIMUM_HARMONIC_DISTANCE); } else { harmonic_region_start = max(harmonic_region_center - new_fundamental_frequency / 2.0f, 0.0f); }; float harmonic_region_start; if (analysis_harmonic_index < spp_frames[frame_index].harmonics.harmonic_count - 1) { harmonic_region_start = min(harmonic_region_center + new_fundamental_frequency / 2.0f, spp_frames[frame_index].harmonics.harmonics[analysis_harmonic_index + 1] - SPP_REGION_RENDERING_STRATEGY_1B_ADAPTIVE_REGION_SEGMENTATION_MINIMUM_HARMONIC_DISTANCE); } else { harmonic_region_start = min(harmonic_region_center + new_fundamental_frequency / 2.0f, NYQUIST_FREQUENCY); }; #endif #if TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_0A_BASIC_STRATEGY || TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_1B_ADAPTIVE_REGION_SEGMENTATION || TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_2C_SPECTRAL_GAP_FILLING float synthesis_region_start = synthesis_harmonic_frequency - (harmonic_region_center - harmonic_region_start); float synthesis_region_end = synthesis_harmonic_frequency + (harmonic_region_end - harmonic_region_center); #if TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_FLOOR || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_CEIL positive int synthesis_region_start_bin = convert_float_int(floorf(synthesis_region_start / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #elif TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_FLOOR || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_CEIL positive int synthesis_region_start_bin = convert_float_int(ceilf(synthesis_region_start / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #else #error #endif #if TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_FLOOR || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_FLOOR positive int synthesis_region_end_bin = convert_float_int(floorf(synthesis_region_end / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #elif TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_CEIL || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_CEIL positive int synthesis_region_end_bin = convert_float_int(ceilf(synthesis_region_end / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #else #error #endif #if TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_2C_SPECTRAL_GAP_FILLING for (positive int itr = spectral_gap_start; itr < synthesis_region_start_bin; itr++) { // XXX: Remember to adjust according to the timbre scaling (and equalization?) when that is re-added float analysis_frequency = itr * synthesis_frequency_step + (harmonic_region_start - synthesis_region_start); const struct narrowband_harmonic *closest_harmonic = search_closest(spp_frames[frame_index].harmonics.harmonics, spp_frames[frame_index].harmonics.harmonic_count, sizeof(struct narrowband_harmonic), analysis_frequency, harmonic_get_frequency); if (fabs(analysis_frequency, closest_harmonic -> frequency) < SPP_REGION_RENDERING_STRATEGY_2C_SPECTRAL_GAP_FILLING_MINIMUM_HARMONIC_DISTANCE) { continue; }; complex float value; #if TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_SPLINE #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC float amplitude = expf(interpolate_spline(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency)); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR float amplitude = interpolate_spline(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); #else #error #endif float phase = interpolate_spline(spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = Z(amplitude, phase); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR float real_part = interpolate_spline(spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); float imaginary_part = interpolate_spline(spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = real_part + imaginary_part * I; #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_LINEAR #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC float amplitude = expf(interpolate_line_segments(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency)); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR float amplitude = interpolate_line_segments(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); #else #error #endif float phase = interpolate_line_segments(spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = Z(amplitude, phase); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR float real_part = interpolate_spline(spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); float imaginary_part = interpolate_line_segments(spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = real_part + imaginary_part * I; #else #error #endif #else #error #endif spectrum[itr] += value; }; #endif for (positive int itr = synthesis_region_start_bin; itr <= synthesis_region_end_bin; itr++) { float analysis_frequency = itr * synthesis_frequency_step + (harmonic_region_start - synthesis_region_start); complex float value; #if TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_SPLINE #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC float amplitude = expf(interpolate_spline(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency)); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR float amplitude = interpolate_spline(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); #else #error #endif float phase = interpolate_spline(spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = Z(amplitude, phase); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR float real_part = interpolate_spline(spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); float imaginary_part = interpolate_spline(spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = real_part + imaginary_part * I; #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_LINEAR #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC float amplitude = expf(interpolate_line_segments(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency)); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR float amplitude = interpolate_line_segments(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); #else #error #endif float phase = interpolate_line_segments(spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = Z(amplitude, phase); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR float real_part = interpolate_spline(spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); float imaginary_part = interpolate_line_segments(spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = real_part + imaginary_part * I; #else #error #endif #else #error #endif spectrum[itr] += value * Z(amplitude_scale, phase_rotation); }; #if TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_2C_SPECTRAL_GAP_FILLING spectral_gap_start = synthesis_region_end_bin + 1; #endif // XXX: Implement strategy 3(D) #elif TEST_SPP_REGION_RENDERING_STRATEGY == SPP_REGION_RENDERING_STRATEGY_4E_FREQUENCY_WARPING // XXX: Can these restrictions be loosened? static_assert(TEST_SPP_SYNTHESIS_HARMONIC_DISTRIBUTION == SPP_SYNTHESIS_HARMONIC_DISTRIBUTION_IDEAL); static_assert(TEST_SPP_ANALYSIS_HARMONIC_DISTRIBUTION == SPP_ANALYSIS_HARMONIC_DISTRIBUTION_IDEAL); static_assert(TEST_SPP_REGION_DELIMITATION_METHOD == SPP_REGION_DELIMITATION_MIDDLE); float synthesis_region_start = synthesis_harmonic_frequency - new_fundamental_frequency / 2.0f; float synthesis_region_end = synthesis_harmonic_frequency + new_fundamental_frequency / 2.0f; #if TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_FLOOR || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_CEIL positive int synthesis_region_start_bin = convert_float_int(floorf(synthesis_region_start / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #elif TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_FLOOR || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_CEIL positive int synthesis_region_start_bin = convert_float_int(ceilf(synthesis_region_start / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #else #error #endif #if TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_FLOOR || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_FLOOR positive int synthesis_region_end_bin = convert_float_int(floorf(synthesis_region_end / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #elif TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_FLOOR_CEIL || TEST_SPP_REGION_RENDERING_ORDER == TEST_SPP_REGION_RENDERING_ORDER_CEIL_CEIL positive int synthesis_region_end_bin = convert_float_int(ceilf(synthesis_region_end / synthesis_frequency_step), 0, SPP_WINDOW_SIZE - 1); #else #error #endif for (positive int itr = synthesis_region_start_bin; itr <= synthesis_region_end_bin; itr++) { #if TEST_SPP_REGION_RENDERING_STRATEGY_4E_FREQUENCY_WARPING_SCALING_TYPE == SPP_REGION_RENDERING_STRATEGY_4E_FREQUENCY_WARPING_TYPE_LINEAR_SCALING float analysis_frequency = ((itr * synthesis_frequency_step) - synthesis_region_start) / (new_fundamental_frequency / spp_frames[frame_index].fundamental_frequency) + harmonic_region_start; #elif TEST_SPP_REGION_RENDERING_STRATEGY_4E_FREQUENCY_WARPING_SCALING_TYPE == SPP_REGION_RENDERING_STRATEGY_4E_FREQUENCY_WARPING_TYPE_NONLINEAR_SCALING float offset = (itr * synthesis_frequency_step) - synthesis_harmonic_frequency; float analysis_frequency = harmonic_region_center + offset + (offset * (spp_frames[frame_index].fundamental_frequency / new_fundamental_frequency) - offset) * (fabs(offset) / new_fundamental_frequency) * 2.0f; #else #error #endif complex float value; #if TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_SPLINE #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC float amplitude = expf(interpolate_spline(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency)); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR float amplitude = interpolate_spline(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); #else #error #endif float phase = interpolate_spline(spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = Z(amplitude, phase); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR float real_part = interpolate_spline(spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); float imaginary_part = interpolate_spline(spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = real_part + imaginary_part * I; #else #error #endif #elif TEST_SPP_INTERPOLATION_TYPE == SPP_INTERPOLATION_LINEAR #if TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_POLAR #if TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LOGARITHMIC float amplitude = expf(interpolate_line_segments(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency)); #elif TEST_SPP_INTERPOLATION_TYPE4 == SPP_INTERPOLATION4_LINEAR float amplitude = interpolate_line_segments(spectrum_spline_amplitude, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); #else #error #endif float phase = interpolate_line_segments(spectrum_spline_phase, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = Z(amplitude, phase); #elif TEST_SPP_INTERPOLATION_TYPE2 == SPP_INTERPOLATION2_RECTANGULAR float real_part = interpolate_spline(spectrum_spline_real_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); float imaginary_part = interpolate_line_segments(spectrum_spline_imaginary_part, DFT_HALF(SPP_PADDED_WINDOW_SIZE), analysis_frequency); value = real_part + imaginary_part * I; #else #error #endif #else #error #endif spectrum[itr] += value * Z(amplitude_scale, phase_rotation); }; // XXX: Implement magnitude/phase gapless version of 4(E) // XXX: Implement patent version of 4(E) // XXX: Implement PHF/PHP #endif }; reconstruct_full_spectrum(spectrum, spectrum, SPP_WINDOW_SIZE); idft(spectrum, window, SPP_WINDOW_SIZE); idftshift(window, SPP_WINDOW_SIZE); unwindow(window, window_function, SPP_WINDOW_SIZE); overlap_and_add(buffer, sample_count, frame_index, window + (SPP_WINDOW_SIZE / 2) - HOP_SIZE); #if TEST_SPP_PHASE_PROPGATION_TIME == SPP_PHASE_PROPGATION_TIME_BEFORE // Intentionally blank #elif TEST_SPP_PHASE_PROPGATION_TIME == SPP_PHASE_PROPGATION_TIME_AFTER #if TEST_SPP_PHASE_PROPGATION_METHOD == SPP_PHASE_PROPGATION_METHOD_CONSTANT fundamental_phase = princarg(fundamental_phase + new_fundamental_frequency * HOP_TIME * M_PIF * 2.0f); #elif TEST_SPP_PHASE_PROPGATION_METHOD == SPP_PHASE_PROPGATION_METHOD_LINEAR #if TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SCALE float next_new_fundamental_frequency; if (frame_index < frame_count - 1) { next_new_fundamental_frequency = spp_frames[frame_index + 1].fundamental_frequency * pitch_scale; } else { next_new_fundamental_frequency = spp_frames[frame_index].fundamental_frequency * pitch_scale; }; #elif TEST_PITCH_TRANSPOSITION_TYPE == PITCH_TRANSPOSITION_SET float next_new_fundamental_frequency = pitch; #else #error #endif fundamental_phase = princarg(fundamental_phase + (new_fundamental_frequency + next_new_fundamental_frequency) * HOP_TIME * M_PIF); // XXX: Also consider the quadratic phase propagation technique from Bonada's 1997 Master's Thesis #else #error #endif #else #error #endif #if TEST_SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS == SPP_USE_PHASE_DISCONTINUITY_HYSTERESIS_YES && TEST_SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE == SPP_PHASE_DISCONTINUITY_HYSTERESIS_PHASE_RESET_TYPE_NEW_HARMONICS prev_new_harmonic_count = new_harmonic_count; #endif }; }