<commit_before>//-----------------------------------------------------------------------------
// Copyright 2012-2015 <NAME>. All Rights Reserved.
// Author: mmorise [at] yamanashi.ac.jp (<NAME>)
//
// Spectral envelope estimation on the basis of the idea of CheapTrick.
//-----------------------------------------------------------------------------
#include "./cheaptrick.h"

#include <math.h>

#include "./common.h"
#include "./constantnumbers.h"
#include "./matlabfunctions.h"

namespace {

//-----------------------------------------------------------------------------
// SmoothingWithRecovery() carries out the spectral smoothing and spectral
// recovery on the Cepstrum domain.
//-----------------------------------------------------------------------------
void SmoothingWithRecovery(double current_f0, int fs, int fft_size, double q1,
    ForwardRealFFT *forward_real_fft, InverseRealFFT *inverse_real_fft,
    double *spectral_envelope) {
// We can control q1 as the parameter. 2015/9/22 by M. Morise
// const double q1 = -0.09;  // Please see the reference in CheapTrick.
  double *smoothing_lifter = new double[fft_size];
  double *compensation_lifter = new double[fft_size];

  smoothing_lifter[0] = 1;
  compensation_lifter[0] = (1.0 - 2.0 * q1) + 2.0 * q1;
  double quefrency;
  for (int i = 1; i <= forward_real_fft->fft_size / 2; ++i) {
    quefrency = static_cast<double>(i) / fs;
    smoothing_lifter[i] = sin(world::kPi * current_f0 * quefrency) /
      (world::kPi * current_f0 * quefrency);
    compensation_lifter[i] = (1.0 - 2.0 * q1) + 2.0 * q1 *
      cos(2.0 * world::kPi * quefrency * current_f0);
  }

  for (int i = 0; i <= fft_size / 2; ++i)
    forward_real_fft->waveform[i] = log(forward_real_fft->waveform[i]);
  for (int i = 1; i < fft_size / 2; ++i)
    forward_real_fft->waveform[fft_size - i] = forward_real_fft->waveform[i];
  fft_execute(forward_real_fft->forward_fft);

  for (int i = 0; i <= fft_size / 2; ++i) {
    inverse_real_fft->spectrum[i][0] = forward_real_fft->spectrum[i][0] *
      smoothing_lifter[i] * compensation_lifter[i] / fft_size;
    inverse_real_fft->spectrum[i][1] = 0.0;
  }
  fft_execute(inverse_real_fft->inverse_fft);

  for (int i = 0; i <= fft_size / 2; ++i)
    spectral_envelope[i] = exp(inverse_real_fft->waveform[i]);

  delete[] smoothing_lifter;
  delete[] compensation_lifter;
}

//-----------------------------------------------------------------------------
// GetPowerSpectrum() calculates the power_spectrum with DC correction.
// DC stands for Direct Current. In this case, the component from 0 to F0 Hz
// is corrected.
//-----------------------------------------------------------------------------
void GetPowerSpectrum(int fs, double current_f0, int fft_size,
    ForwardRealFFT *forward_real_fft) {
  int half_window_length = matlab_round(1.5 * fs / current_f0);

  // FFT
  for (int i = half_window_length * 2 + 1; i < fft_size; ++i)
    forward_real_fft->waveform[i] = 0.0;
  fft_execute(forward_real_fft->forward_fft);

  // Calculation of the power spectrum.
  double *power_spectrum = forward_real_fft->waveform;
  for (int i = 0; i <= fft_size / 2; ++i)
    power_spectrum[i] =
      forward_real_fft->spectrum[i][0] * forward_real_fft->spectrum[i][0] +
      forward_real_fft->spectrum[i][1] * forward_real_fft->spectrum[i][1];

  // DC correction
  DCCorrection(power_spectrum, current_f0, fs, fft_size, power_spectrum);
}

//-----------------------------------------------------------------------------
// SetParametersForGetWindowedWaveform()
//-----------------------------------------------------------------------------
void SetParametersForGetWindowedWaveform(int half_window_length, int x_length,
    double temporal_position, int fs, double current_f0, int *base_index,
    int *index, double *window) {
  for (int i = -half_window_length; i <= half_window_length; ++i)
    base_index[i + half_window_length] = i;
  for (int i = 0; i <= half_window_length * 2; ++i)
    index[i] = MyMin(x_length - 1, MyMax(0,
        matlab_round(temporal_position * fs + base_index[i])));

  // Designing of the window function<commit_msg>Updated cheaptrick.cpp to use MyMin<Type> and MyMax<Type>.<commit_after>//-----------------------------------------------------------------------------
// Copyright 2012-2015 <NAME>. All Rights Reserved.
// Author: mmorise [at] yamanashi.ac.jp (<NAME>)
//
// Spectral envelope estimation on the basis of the idea of CheapTrick.
//-----------------------------------------------------------------------------
#include "./cheaptrick.h"

#include <math.h>

#include "./common.h"
#include "./constantnumbers.h"
#include "./matlabfunctions.h"

namespace {

//-----------------------------------------------------------------------------
// SmoothingWithRecovery() carries out the spectral smoothing and spectral
// recovery on the Cepstrum domain.
//-----------------------------------------------------------------------------
void SmoothingWithRecovery(double current_f0, int fs, int fft_size, double q1,
    ForwardRealFFT *forward_real_fft, InverseRealFFT *inverse_real_fft,
    double *spectral_envelope) {
// We can control q1 as the parameter. 2015/9/22 by <NAME>
// const double q1 = -0.09;  // Please see the reference in CheapTrick.
  double *smoothing_lifter = new double[fft_size];
  double *compensation_lifter = new double[fft_size];

  smoothing_lifter[0] = 1;
  compensation_lifter[0] = (1.0 - 2.0 * q1) + 2.0 * q1;
  double quefrency;
  for (int i = 1; i <= forward_real_fft->fft_size / 2; ++i) {
    quefrency = static_cast<double>(i) / fs;
    smoothing_lifter[i] = sin(world::kPi * current_f0 * quefrency) /
      (world::kPi * current_f0 * quefrency);
    compensation_lifter[i] = (1.0 - 2.0 * q1) + 2.0 * q1 *
      cos(2.0 * world::kPi * quefrency * current_f0);
  }

  for (int i = 0; i <= fft_size / 2; ++i)
    forward_real_fft->waveform[i] = log(forward_real_fft->waveform[i]);
  for (int i = 1; i < fft_size / 2; ++i)
    forward_real_fft->waveform[fft_size - i] = forward_real_fft->waveform[i];
  fft_execute(forward_real_fft->forward_fft);

  for (int i = 0; i <= fft_size / 2; ++i) {
    inverse_real_fft->spectrum[i][0] = forward_real_fft->spectrum[i][0] *
      smoothing_lifter[i] * compensation_lifter[i] / fft_size;
    inverse_real_fft->spectrum[i][1] = 0.0;
  }
  fft_execute(inverse_real_fft->inverse_fft);

  for (int i = 0; i <= fft_size / 2; ++i)
    spectral_envelope[i] = exp(inverse_real_fft->waveform[i]);

  delete[] smoothing_lifter;
  delete[] compensation_lifter;
}

//-----------------------------------------------------------------------------
// GetPowerSpectrum() calculates the power_spectrum with DC correction.
// DC stands for Direct Current. In this case, the component from 0 to F0 Hz
// is corrected.
//-----------------------------------------------------------------------------
void GetPowerSpectrum(int fs, double current_f0, int fft_size,
    ForwardRealFFT *forward_real_fft) {
  int half_window_length = matlab_round(1.5 * fs / current_f0);

  // FFT
  for (int i = half_window_length * 2 + 1; i < fft_size; ++i)
    forward_real_fft->waveform[i] = 0.0;
  fft_execute(forward_real_fft->forward_fft);

  // Calculation of the power spectrum.
  double *power_spectrum = forward_real_fft->waveform;
  for (int i = 0; i <= fft_size / 2; ++i)
    power_spectrum[i] =
      forward_real_fft->spectrum[i][0] * forward_real_fft->spectrum[i][0] +
      forward_real_fft->spectrum[i][1] * forward_real_fft->spectrum[i][1];

  // DC correction
  DCCorrection(power_spectrum, current_f0, fs, fft_size, power_spectrum);
}

//-----------------------------------------------------------------------------
// SetParametersForGetWindowedWaveform()
//-----------------------------------------------------------------------------
void SetParametersForGetWindowedWaveform(int half_window_length, int x_length,
    double temporal_position, int fs, double current_f0, int *base_index,
    int *index, double *window) {
  for (int i = -half_window_length; i <= half_window_length; ++i)
    base_index[i + half_window_length] = i;
  for (int i = 0; i <= half_window_length * 2; ++i)
    index[i] = MyMinInt(x_length - 1, MyMaxInt(0,
        matlab_round(temporal_position * fs + base_index[i])));

  // Designing of the window function
{'fraction_non_alphanumeric': 0.26577552986512526, 'fraction_numerical': 0.019990366088631986, 'mean_word_length': 2.1000373273609556, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1503476', 'n_tokens_mistral': 2830, 'n_tokens_neox': 2661, 'n_words': 620}