fork download
copy 
  1. #include <bits/stdc++.h>
  2.  
  3. using namespace std;
  4.  
  5. /****************************************************************************
  6. *
  7. * NAME: PitchShift.cs
  8. * VERSION: 1.2
  9. * HOME URL: http://w...content-available-to-author-only...n.com
  10. * KNOWN BUGS: none
  11. *
  12. * SYNOPSIS:
  13. *
  14. * Routine for doing pitch shifting while maintaining duration
  15. * using the Short-Time Fourier Transform.
  16. *
  17. * DESCRIPTION:
  18. *
  19. * The routine takes:
  20. *   1. a pitchShift factor value
  21. *      between 0.5 (one octave down) and 2. (one octave up).
  22. *      A value of exactly 1 does not change the pitch.
  23. *   2. numSampsToProcess
  24. *      tells the routine how many samples in indata[0...numSampsToProcess-1]
  25. *      should be pitch shifted and moved to outdata[0...numSampsToProcess-1].
  26. *      The two buffers can be identical (ie. it can process the data in-place).
  27. *   3. fftFrameSize
  28. *      defines the FFT frame size used for the processing.
  29. *      Typical values are 1024, 2048 and 4096.
  30. *      It may be any value <= MAX_FRAME_LENGTH, but it MUST be a power of 2.
  31. *   4. osamp
  32. *      is the STFT oversampling factor
  33. *      which also determines the overlap between adjacent STFT frames.
  34. *      It should at least be 4 for moderate scaling ratios.
  35. *      A value of 32 is recommended for best quality.
  36. *   5. sampleRate
  37. *      is the sample rate for the signal in unit Hz, ie. 44100 for 44.1 kHz audio.
  38. *
  39. *   The data passed to the routine in indata[] should be in the range [-1.0, 1.0),
  40. *   which is also the output range for the data, make sure the data is scaled accordingly.
  41. *   For 16-bit signed integers, the input data would have to be divided by 32768, and the
  42. *   output data multiplied by the same number.
  43. *
  44. * COPYRIGHT 1999-2006 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
  45. *
  46. * 						The Wide Open License (WOL)
  47. *
  48. * Permission to use, copy, modify, distribute and sell this software and its
  49. * documentation for any purpose is hereby granted without fee, provided that
  50. * the above copyright notice and this license appear in all source copies.
  51. * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
  52. * ANY KIND. See http://w...content-available-to-author-only...u.com/wol.htm for more information.
  53. *
  54. *****************************************************************************/
  55.  
  56. /****************************************************************************
  57. *
  58. * This code was converted to C# by Michael Knight
  59. * madmik3 at gmail dot com.
  60. * http://sites.google.com/site/mikescoderama/
  61. *
  62. *****************************************************************************/
  63.  
  64. /*****************************************************************************
  65. *
  66. * This code is translated to C++17 By CodingKnight at Codeforces.com
  67. * On August 20, 2018
  68. *
  69. * UPDATES:
  70. *
  71. * 1. The MAX_FRAME_LENGTH condition on the fftFrameSize has been removed.
  72. * 2. All data arrays have been declared using vector_t data structure
  73. * 3. A vector_t window object has been declared and precomputed outside
  74. *    the window interleaving loop, as the array is independent of the data
  75. *    samples, and depends on the fftFrameSize only.
  76. * 4. C++ STL complex<ldbl> library functions arg(), abs(), real(), and imag()
  77. *    have been used to deal with complex numbers.
  78. * 5. The Pitch Shifter has been wrapped up inside a C++ class Phase_Shifter_t
  79. *    with public constructor:
  80. *
  81. *          Pitch_Shifter_t( int fftSampleSize = 2048, int osamp = 10 )
  82. *
  83. *    and public data processing function:
  84. *
  85. *          void pitch_shift( ldbl pitchShift, int numSampsToProcess, ldbl sampleRate, vector_t& outdata )
  86. *
  87. *
  88. *****************************************************************************/
  89.  
  90. typedef long double ldbl;
  91.  
  92. struct vector_t: vector< ldbl >
  93. {
  94.     void initialize( int n )
  95.     {
  96.         resize( n ); for( auto &x: *this ) x = 0.0;
  97.     }
  98.  
  99.     void swap( int i, int j )
  100.     {
  101.         std::swap( at( i ), at( j ) );
  102.     }
  103. };
  104.  
  105. class Pitch_Shifter_t: public vector_t
  106. {
  107.     static constexpr ldbl pi = acos( -1.0 ), twice_pi = 2.0 * pi;
  108.  
  109.     vector_t gInFIFO, gOutFIFO, gAnaFreq, gAnaMagn, gSynFreq, gSynMagn;
  110.  
  111.     vector_t gFFTworksp, gOutputAccum, gLastPhase, gSumPhase, window;
  112.  
  113.     int fftFrameSize, fftFrameSize1, fftFrameSize2, gRover, inFifoLatency, k_max, osamp, stepSize;
  114.  
  115.     ldbl osamp_factor;
  116.  
  117.     void Short_Time_Fourier_Transform( ldbl sign )
  118.     {
  119.  
  120. /*
  121.         FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)
  122.         Fills gFFTwroksp[0...2*fftFrameSize-1] with the Fourier transform of the
  123.         time domain data in fftBuffer[0...2*fftFrameSize-1]. The FFT array takes
  124.         and returns the cosine and sine parts in an interleaved manner, ie.
  125.         gFFTworksp[0] = cosPart[0], gFFTworksp[1] = sinPart[0], asf. fftFrameSize
  126.         must be a power of 2. It expects a complex input signal (see footnote 2),
  127.         ie. when working with 'common' audio signals our input signal has to be
  128.         passed as {in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the transform
  129.         of the frequencies of interest is in gFFTworksp[0...fftFrameSize].
  130. */
  131.         for ( int i = 2, k = fftFrameSize1 - 2; i < k; i += 2 )
  132.         {
  133.             int j = 0;
  134.  
  135.             for ( int bit_mask = 2; bit_mask < fftFrameSize1; bit_mask <<= 1, j <<= 1 )
  136.                 if ( i & bit_mask )
  137.                     j++;
  138.  
  139.             if ( i < j )
  140.                 gFFTworksp.swap( i, j ),
  141.                 gFFTworksp.swap( i + 1, j + 1 );
  142.         }
  143.  
  144.         for ( int k = 0, le = 4, le2 = 2, le1 = 1; k < k_max; k++, le1 = le2, le2 = le, le <<= 1 )
  145.         {
  146.             ldbl   ur = 1.0,
  147.                     ui = 0.0,
  148.                     arg = sign * pi / ldbl( le1 );
  149.  
  150.             auto z = polar( ur, arg );
  151.  
  152.             ldbl   wr = real( z ),
  153.                     wi = imag( z );
  154.  
  155.             for ( int j = 0; j < le2; j += 2 )
  156.             {
  157.                 for ( int i = j, l = i + 1; i < fftFrameSize1; i += le, l += le )
  158.                 {
  159.                     int     u = i + le2,
  160.                             v = u + 1;
  161.  
  162.                     ldbl   tr = gFFTworksp[ u ] * ur - gFFTworksp[ v ] * ui,
  163.                             ti = gFFTworksp[ u ] * ui + gFFTworksp[ v ] * ur;
  164.  
  165.                     gFFTworksp[ u ] = gFFTworksp[ i ] - tr,
  166.                     gFFTworksp[ v ] = gFFTworksp[ l ] - ti;
  167.  
  168.                     gFFTworksp[ i ] += tr,
  169.                     gFFTworksp[ l ] += ti;
  170.                 }
  171.  
  172.                 ldbl tr = ur * wr - ui * wi;
  173.  
  174.                 ui = ur * wi + ui * wr,
  175.                 ur = tr;
  176.             }
  177.         }
  178.     }
  179.  
  180. public:
  181.  
  182.     Pitch_Shifter_t( int _fftFrameSize = 2048, int _osamp = 10 ) : fftFrameSize( _fftFrameSize ), osamp( _osamp )
  183.     {
  184.  
  185. /*
  186.         Purpose: doing pitch shifting while maintaining duration using the Short Time Fourier Transform.
  187.         Author: (c)1999-2015 Stephan M. Bernsee <s.bernsee [AT] zynaptiq [DOT] com>
  188. */
  189.         fftFrameSize1 = fftFrameSize * 2,
  190.         fftFrameSize2 = fftFrameSize / 2,
  191.         stepSize = fftFrameSize / osamp,
  192.         inFifoLatency = gRover = fftFrameSize - stepSize,
  193.         k_max = int( log( fftFrameSize ) / log( 2.0 ) + 0.5 ),
  194.         osamp_factor =  2.0 / ( ldbl( fftFrameSize2 ) * ldbl( osamp ) );
  195.  
  196.         initialize( fftFrameSize );
  197.  
  198.         int max_phase_size = fftFrameSize2 + 1;
  199.  
  200.           gInFIFO.initialize( fftFrameSize ),
  201.          gOutFIFO.initialize( fftFrameSize ),
  202.          gAnaFreq.initialize( fftFrameSize ),
  203.          gAnaMagn.initialize( fftFrameSize ),
  204.          gSynFreq.initialize( fftFrameSize ),
  205.          gSynMagn.initialize( fftFrameSize );
  206.  
  207.           gFFTworksp.initialize( fftFrameSize1 ),
  208.         gOutputAccum.initialize( fftFrameSize1 );
  209.  
  210.         gLastPhase.initialize( max_phase_size ),
  211.          gSumPhase.initialize( max_phase_size );
  212.  
  213.         window.initialize( fftFrameSize );
  214.  
  215.         for( int k = 0; k < fftFrameSize; k++ )
  216.             window[ k ] = 0.5 - .5 * cos( twice_pi *  ldbl( k ) / ldbl( fftFrameSize ) );
  217.     }
  218.  
  219.     void pitch_shift( ldbl pitchShift, int numSampsToProcess, ldbl sampleRate, vector_t& outdata )
  220.     {
  221.         /* set up some handy variables */
  222.  
  223.         ldbl   freqPerBin = sampleRate / ldbl( fftFrameSize ),
  224.                 expct = twice_pi * ldbl( stepSize ) / ldbl( fftFrameSize );
  225.  
  226.         /* main processing loop */
  227.  
  228.         for ( int i = 0; i < numSampsToProcess; i++ )
  229.         {
  230.             /* As long as we have not yet collected enough data just read in */
  231.  
  232.             gInFIFO[ gRover ] = at( i ),
  233.             outdata[ i ] = gOutFIFO[ gRover - inFifoLatency ],
  234.             gRover++;
  235.  
  236.             /* now we have enough data for processing */
  237.  
  238.             if ( gRover >= fftFrameSize )
  239.             {
  240.                 gRover = inFifoLatency;
  241.  
  242.                 /* do windowing and re,im interleave */
  243.  
  244.                 for ( int k = 0, l = 0, m = 1; k < fftFrameSize; k++, l += 2, m += 2 )
  245.                     gFFTworksp[ l ] = gInFIFO[ k ] * window[ k ],
  246.                     gFFTworksp[ m ] = 0.0;
  247.  
  248.                 /* ***************** ANALYSIS ******************* */
  249.                 /* do transform */
  250.  
  251.                 Short_Time_Fourier_Transform( -1.0 );
  252.  
  253.                 /* this is the analysis step */
  254.  
  255.                 for ( int k = 0, l = 0, m = 1; k <= fftFrameSize2; k++, l += 2, m += 2 )
  256.                 {
  257.  
  258.                     /* de-interlace FFT buffer */
  259.  
  260.                     complex<ldbl> z( gFFTworksp[ l ], gFFTworksp[ m ] );
  261.  
  262.                     /* compute magnitude and phase */
  263.                     ldbl   magn  = 2.0 * abs( z ),
  264.                             phase = arg( z );
  265.  
  266.                     /* compute phase difference */
  267.                     ldbl tmp = phase - gLastPhase[ k ];
  268.  
  269.                     gLastPhase[ k ] = phase;
  270.  
  271.                     /* subtract expected phase difference */
  272.                     tmp -= ldbl( k ) * expct;
  273.  
  274.                     /* map delta phase into +/- Pi interval */
  275.                     int qpd = int( tmp / pi );
  276.  
  277.                     if ( qpd >= 0 )
  278.                         qpd += qpd & 1;
  279.                     else
  280.                         qpd -= qpd & 1;
  281.  
  282.                     tmp -= pi * qpd;
  283.  
  284.                     /* get deviation from bin frequency from the +/- Pi interval */
  285.                     tmp = osamp * tmp / twice_pi;
  286.  
  287.                     /* compute the k-th partials' true frequency */
  288.                     tmp = ldbl( k ) * freqPerBin + tmp * freqPerBin;
  289.  
  290.                     /* store magnitude and true frequency in analysis arrays */
  291.                     gAnaMagn[ k ] = magn,
  292.                     gAnaFreq[ k ] = tmp;
  293.                 }
  294.  
  295.                 /* ***************** PROCESSING ******************* */
  296.                 /* this does the actual pitch shifting */
  297.  
  298.                 for ( int zero = 0; zero < fftFrameSize; zero++ )
  299.                     gSynMagn[ zero ] = 0,
  300.                     gSynFreq[ zero ] = 0;
  301.  
  302.                 for ( int k = 0; k <= fftFrameSize2; k++ )
  303.                 {
  304.                     int index = int( ldbl( k ) * pitchShift );
  305.  
  306.                     if ( index <= fftFrameSize2 )
  307.                         gSynMagn[ index ] += gAnaMagn[ k ],
  308.                         gSynFreq[ index ]  = gAnaFreq[ k ] * pitchShift;
  309.                     else
  310.                         break;
  311.                 }
  312.  
  313.                 /* ***************** SYNTHESIS ******************* */
  314.                 /* this is the synthesis step */
  315.                 for ( int k = 0, l = 0, m = 1; k <= fftFrameSize2; k++, l += 2, m += 2 )
  316.                 {
  317.  
  318.                     /* get true frequency from synthesis arrays */
  319.                     ldbl tmp  = gSynFreq[ k ];
  320.  
  321.                     /* subtract bin mid frequency */
  322.                     tmp -= ldbl( k ) * freqPerBin;
  323.  
  324.                     /* get bin deviation from freq deviation */
  325.                     tmp /= freqPerBin;
  326.  
  327.                     /* take osamp into account */
  328.                     tmp = twice_pi * tmp / osamp;
  329.  
  330.                     /* add the overlap phase advance back in */
  331.                     tmp += ldbl( k ) * expct;
  332.  
  333.                     /* accumulate delta phase to get bin phase */
  334.                     gSumPhase[ k ] += tmp;
  335.  
  336.                     auto z = polar( gSynMagn[ k ], gSumPhase[ k ] );
  337.  
  338.                     /* get real and imag part and re-interleave */
  339.                     gFFTworksp[ l ] = real( z ),
  340.                     gFFTworksp[ m ] = imag( z );
  341.                 }
  342.  
  343.                 /* zero negative frequencies */
  344.                 for ( int k = fftFrameSize + 2; k < fftFrameSize1; k++ )
  345.                     gFFTworksp[ k ] = 0.0;
  346.  
  347.                 /* do inverse transform */
  348.                 Short_Time_Fourier_Transform( 1.0 );
  349.  
  350.                 /* do windowing and add to output accumulator */
  351.                 for ( int k = 0, l = 0; k < fftFrameSize; k++, l += 2 )
  352.                     gOutputAccum[ k ] += osamp_factor * window[ k ] * gFFTworksp[ l ];
  353.  
  354.                 for ( int k = 0; k < stepSize; k++ )
  355.                     gOutFIFO[ k ] = gOutputAccum[ k ];
  356.  
  357.                 /* shift accumulator */
  358.                 //memmove(gOutputAccum, gOutputAccum + stepSize, fftFrameSize * sizeof(ldbl));
  359.                 for ( int k = 0, l = stepSize; k < fftFrameSize; k++, l++ )
  360.                     gOutputAccum[ k ] = gOutputAccum[ l ];
  361.  
  362.                 /* move input FIFO */
  363.                 for ( int k = 0, l = stepSize; k < inFifoLatency; k++, l++ )
  364.                     gInFIFO[ k ] = gInFIFO[ l ];
  365.             }
  366.         }
  367.     }
  368. };
  369.  
  370. typedef long long ll;
  371.  
  372. typedef chrono::high_resolution_clock Clock_t;
  373.  
  374. typedef Clock_t::time_point      time_point_t;
  375.  
  376. time_point_t Now() { return Clock_t::now(); }
  377.  
  378. struct Timer_t
  379. {
  380.     time_point_t start_time;
  381.  
  382.     Timer_t() : start_time( Now() ) {}
  383.  
  384.     typedef long double ldbl;
  385.  
  386.     typedef chrono::nanoseconds nanosec_t;
  387.  
  388.     ll elapsed_time()
  389.     {
  390.         nanosec_t timer = chrono::duration_cast< nanosec_t >( Now() - start_time );
  391.  
  392.         return timer.count() / 1000000LL;
  393.     }
  394. };
  395.  
  396. struct random_t: mt19937_64
  397. {
  398.     random_t() : mt19937_64( Clock_t::now().time_since_epoch().count() ) {}
  399.  
  400.     int random( int MIN, int MAX )
  401.     {
  402.         ll x_min = MIN, x_max = MAX, interval = x_max - x_min + 1, number = (*this)();
  403.  
  404.         if ( number < 0 )
  405.             number += LLONG_MAX;
  406.  
  407.         return x_min + ( number % interval );
  408.     }
  409.  
  410. } mt;
  411.  
  412. int main()
  413. {
  414.     ios_base::sync_with_stdio( false ), cin.tie( nullptr ), cout.tie( nullptr );
  415.  
  416.     Pitch_Shifter_t x; vector_t y;
  417.  
  418.     int n = x.size(), r_max = INT_MAX / 2, r_min = -r_max;
  419.  
  420.     ldbl pitch_shift, sample_rate, range = r_max;
  421.  
  422.     ldbl x_avg = 0.0, y_avg = 0.0, x_min = FLT_MAX, y_min = FLT_MAX, x_max = -FLT_MAX, y_max = -FLT_MAX;
  423.  
  424.     for( int i = 0; i < n; i++ )
  425.         x[ i ] = ldbl( mt.random( r_min, r_max ) ) / range,
  426.         x_min  = min( x_min, x[ i ] ),
  427.         x_max  = max( x_max, x[ i ] ),
  428.         x_avg += x[ i ] / ldbl( n );
  429.  
  430.     y.initialize( n );
  431.  
  432.     cin >> pitch_shift >> sample_rate;
  433.  
  434.     Timer_t timer;
  435.  
  436.     x.pitch_shift( pitch_shift, n, sample_rate, y );
  437.  
  438.     ll t = timer.elapsed_time();
  439.  
  440.     for( int i = 0; i < n; i++ )
  441.         y_min  = min( y_min, y[ i ] ),
  442.         y_max  = max( y_max, y[ i ] );
  443.  
  444.     ldbl scaling_factor = 2.0 / ( y_max - y_min ), g = y_min;
  445.  
  446.     y_min = FLT_MAX, y_max = -FLT_MAX;
  447.  
  448.     for( int i = 0; i < n; i++ )
  449.         y[ i ] = scaling_factor * ( y[ i ] - g ) - 1.0,
  450.         y_min  = min( y_min, y[ i ] ),
  451.         y_max  = max( y_max, y[ i ] ),
  452.         y_avg += y[ i ] / ldbl( n );
  453.  
  454.     cout << "x_avg = " << x_avg << ", x_min = " << x_min << ", x_max = " << x_max << endl,
  455.     cout << "y_avg = " << y_avg << ", y_min = " << y_min << ", y_max = " << y_max << endl,
  456.     cout << "scaling factor = " << scaling_factor << endl,
  457.     cout << "elapsed time = " << t << " msec.";
  458. }
  459.  
Success #stdin #stdout 0.02s 4404KB
comments (?)
stdin
copy 
0.5 44100
stdout
copy 
x_avg = 0.028723, x_min = -0.999108, x_max = 0.998987
y_avg = -0.765273, y_min = -1, y_max = 1
scaling factor = 0.275593
elapsed time = 12 msec.
https://ideone.com/sltOKf
Pitch Shifter
language:
C++14 (gcc 8.3)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!