Something went wrong. Try again.
Monorepo for Aesthetic.Computer aesthetic.computer
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229//-----------------------------------------------------------------------------// Copyright 2012 Masanori Morise// Author: mmorise [at] meiji.ac.jp (Masanori Morise)// Last update: 2021/02/15//// common.cpp includes functions used in at least two files.// (1) Common functions// (2) FFT, IFFT and minimum phase analysis.//// In FFT analysis and minimum phase analysis,// Functions "Initialize*()" allocate the mamory.// Functions "Destroy*()" free the accolated memory.// FFT size is used for initialization, and structs are used to keep the memory.// Functions "GetMinimumPhaseSpectrum()" calculate minimum phase spectrum.// Forward and inverse FFT do not have the function "Get*()",// because forward FFT and inverse FFT can run in one step.////-----------------------------------------------------------------------------#include "world/common.h"
#include <math.h>
#include "world/constantnumbers.h"#include "world/matlabfunctions.h"
namespace {static void SetParametersForLinearSmoothing(int boundary, int fft_size, int fs, double width, const double *power_spectrum, double *mirroring_spectrum, double *mirroring_segment, double *frequency_axis) { for (int i = 0; i < boundary; ++i) mirroring_spectrum[i] = power_spectrum[boundary - i]; for (int i = boundary; i < fft_size / 2 + boundary; ++i) mirroring_spectrum[i] = power_spectrum[i - boundary]; for (int i = fft_size / 2 + boundary; i <= fft_size / 2 + boundary * 2; ++i) mirroring_spectrum[i] = power_spectrum[fft_size / 2 - (i - (fft_size / 2 + boundary))];
mirroring_segment[0] = mirroring_spectrum[0] * fs / fft_size; for (int i = 1; i < fft_size / 2 + boundary * 2 + 1; ++i) mirroring_segment[i] = mirroring_spectrum[i] * fs / fft_size + mirroring_segment[i - 1];
for (int i = 0; i <= fft_size / 2; ++i) frequency_axis[i] = static_cast<double>(i) / fft_size * fs - width / 2.0;}
} // namespace
//-----------------------------------------------------------------------------int GetSuitableFFTSize(int sample) { return static_cast<int>(pow(2.0, static_cast<int>(log(static_cast<double>(sample)) / world::kLog2) + 1.0));}
void DCCorrection(const double *input, double f0, int fs, int fft_size, double *output) { int upper_limit = 2 + static_cast<int>(f0 * fft_size / fs); double *low_frequency_replica = new double[upper_limit]; double *low_frequency_axis = new double[upper_limit];
for (int i = 0; i < upper_limit; ++i) low_frequency_axis[i] = static_cast<double>(i) * fs / fft_size;
int upper_limit_replica = upper_limit - 1; interp1Q(f0 - low_frequency_axis[0], -static_cast<double>(fs) / fft_size, input, upper_limit + 1, low_frequency_axis, upper_limit_replica, low_frequency_replica);
for (int i = 0; i < upper_limit_replica; ++i) output[i] = input[i] + low_frequency_replica[i];
delete[] low_frequency_replica; delete[] low_frequency_axis;}
void LinearSmoothing(const double *input, double width, int fs, int fft_size, double *output) { int boundary = static_cast<int>(width * fft_size / fs) + 1;
// These parameters are set by the other function. double *mirroring_spectrum = new double[fft_size / 2 + boundary * 2 + 1]; double *mirroring_segment = new double[fft_size / 2 + boundary * 2 + 1]; double *frequency_axis = new double[fft_size / 2 + 1]; SetParametersForLinearSmoothing(boundary, fft_size, fs, width, input, mirroring_spectrum, mirroring_segment, frequency_axis);
double *low_levels = new double[fft_size / 2 + 1]; double *high_levels = new double[fft_size / 2 + 1]; double origin_of_mirroring_axis = -(boundary - 0.5) * fs / fft_size; double discrete_frequency_interval = static_cast<double>(fs) / fft_size;
interp1Q(origin_of_mirroring_axis, discrete_frequency_interval, mirroring_segment, fft_size / 2 + boundary * 2 + 1, frequency_axis, fft_size / 2 + 1, low_levels);
for (int i = 0; i <= fft_size / 2; ++i) frequency_axis[i] += width;
interp1Q(origin_of_mirroring_axis, discrete_frequency_interval, mirroring_segment, fft_size / 2 + boundary * 2 + 1, frequency_axis, fft_size / 2 + 1, high_levels);
for (int i = 0; i <= fft_size / 2; ++i) output[i] = (high_levels[i] - low_levels[i]) / width;
delete[] mirroring_spectrum; delete[] mirroring_segment; delete[] frequency_axis; delete[] low_levels; delete[] high_levels;}
void NuttallWindow(int y_length, double *y) { double tmp; for (int i = 0; i < y_length; ++i) { tmp = i / (y_length - 1.0); y[i] = 0.355768 - 0.487396 * cos(2.0 * world::kPi * tmp) + 0.144232 * cos(4.0 * world::kPi * tmp) - 0.012604 * cos(6.0 * world::kPi * tmp); }}
//-----------------------------------------------------------------------------// FFT, IFFT and minimum phase analysisvoid InitializeForwardRealFFT(int fft_size, ForwardRealFFT *forward_real_fft) { forward_real_fft->fft_size = fft_size; forward_real_fft->waveform = new double[fft_size]; forward_real_fft->spectrum = new fft_complex[fft_size]; forward_real_fft->forward_fft = fft_plan_dft_r2c_1d(fft_size, forward_real_fft->waveform, forward_real_fft->spectrum, FFT_ESTIMATE);}
void DestroyForwardRealFFT(ForwardRealFFT *forward_real_fft) { fft_destroy_plan(forward_real_fft->forward_fft); delete[] forward_real_fft->spectrum; delete[] forward_real_fft->waveform;}
void InitializeInverseRealFFT(int fft_size, InverseRealFFT *inverse_real_fft) { inverse_real_fft->fft_size = fft_size; inverse_real_fft->waveform = new double[fft_size]; inverse_real_fft->spectrum = new fft_complex[fft_size]; inverse_real_fft->inverse_fft = fft_plan_dft_c2r_1d(fft_size, inverse_real_fft->spectrum, inverse_real_fft->waveform, FFT_ESTIMATE);}
void DestroyInverseRealFFT(InverseRealFFT *inverse_real_fft) { fft_destroy_plan(inverse_real_fft->inverse_fft); delete[] inverse_real_fft->spectrum; delete[] inverse_real_fft->waveform;}
void InitializeInverseComplexFFT(int fft_size, InverseComplexFFT *inverse_complex_fft) { inverse_complex_fft->fft_size = fft_size; inverse_complex_fft->input = new fft_complex[fft_size]; inverse_complex_fft->output = new fft_complex[fft_size]; inverse_complex_fft->inverse_fft = fft_plan_dft_1d(fft_size, inverse_complex_fft->input, inverse_complex_fft->output, FFT_BACKWARD, FFT_ESTIMATE);}
void DestroyInverseComplexFFT(InverseComplexFFT *inverse_complex_fft) { fft_destroy_plan(inverse_complex_fft->inverse_fft); delete[] inverse_complex_fft->input; delete[] inverse_complex_fft->output;}
void InitializeMinimumPhaseAnalysis(int fft_size, MinimumPhaseAnalysis *minimum_phase) { minimum_phase->fft_size = fft_size; minimum_phase->log_spectrum = new double[fft_size]; minimum_phase->minimum_phase_spectrum = new fft_complex[fft_size]; minimum_phase->cepstrum = new fft_complex[fft_size]; minimum_phase->inverse_fft = fft_plan_dft_r2c_1d(fft_size, minimum_phase->log_spectrum, minimum_phase->cepstrum, FFT_ESTIMATE); minimum_phase->forward_fft = fft_plan_dft_1d(fft_size, minimum_phase->cepstrum, minimum_phase->minimum_phase_spectrum, FFT_FORWARD, FFT_ESTIMATE);}
void GetMinimumPhaseSpectrum(const MinimumPhaseAnalysis *minimum_phase) { // Mirroring for (int i = minimum_phase->fft_size / 2 + 1; i < minimum_phase->fft_size; ++i) minimum_phase->log_spectrum[i] = minimum_phase->log_spectrum[minimum_phase->fft_size - i];
// This fft_plan carries out "forward" FFT. // To carriy out the Inverse FFT, the sign of imaginary part // is inverted after FFT. fft_execute(minimum_phase->inverse_fft); minimum_phase->cepstrum[0][1] *= -1.0; for (int i = 1; i < minimum_phase->fft_size / 2; ++i) { minimum_phase->cepstrum[i][0] *= 2.0; minimum_phase->cepstrum[i][1] *= -2.0; } minimum_phase->cepstrum[minimum_phase->fft_size / 2][1] *= -1.0; for (int i = minimum_phase->fft_size / 2 + 1; i < minimum_phase->fft_size; ++i) { minimum_phase->cepstrum[i][0] = 0.0; minimum_phase->cepstrum[i][1] = 0.0; }
fft_execute(minimum_phase->forward_fft);
// Since x is complex number, calculation of exp(x) is as following. // Note: This FFT library does not keep the aliasing. double tmp; for (int i = 0; i <= minimum_phase->fft_size / 2; ++i) { tmp = exp(minimum_phase->minimum_phase_spectrum[i][0] / minimum_phase->fft_size); minimum_phase->minimum_phase_spectrum[i][0] = tmp * cos(minimum_phase->minimum_phase_spectrum[i][1] / minimum_phase->fft_size); minimum_phase->minimum_phase_spectrum[i][1] = tmp * sin(minimum_phase->minimum_phase_spectrum[i][1] / minimum_phase->fft_size); }}
void DestroyMinimumPhaseAnalysis(MinimumPhaseAnalysis *minimum_phase) { fft_destroy_plan(minimum_phase->forward_fft); fft_destroy_plan(minimum_phase->inverse_fft); delete[] minimum_phase->cepstrum; delete[] minimum_phase->log_spectrum; delete[] minimum_phase->minimum_phase_spectrum;}