LashaBukhnikashvili's blog

By LashaBukhnikashvili, history, 6 years ago, In English

Hello, I'm working on Audio Algorithm: Pitch-Shifter, unfortunately I can't find Text material about it design, But I found one of the code, implemented on c#: pastebin link.

I would post code here as well, and if you have any idea how I can reduce time or implement new one, please share it, Thanks in advance.

/****************************************************************************
*
* NAME: PitchShift.cs
* VERSION: 1.2
* HOME URL: http://www.dspdimension.com
* KNOWN BUGS: none
*
* SYNOPSIS: Routine for doing pitch shifting while maintaining
* duration using the Short Time Fourier Transform.
*
* DESCRIPTION: The routine takes a pitchShift factor value which is between 0.5
* (one octave down) and 2. (one octave up). A value of exactly 1 does not change
* the pitch. numSampsToProcess tells the routine how many samples in indata[0...
* numSampsToProcess-1] should be pitch shifted and moved to outdata[0 ...
* numSampsToProcess-1]. The two buffers can be identical (ie. it can process the
* data in-place). fftFrameSize defines the FFT frame size used for the
* processing. Typical values are 1024, 2048 and 4096. It may be any value <=
* MAX_FRAME_LENGTH but it MUST be a power of 2. osamp is the STFT
* oversampling factor which also determines the overlap between adjacent STFT
* frames. It should at least be 4 for moderate scaling ratios. A value of 32 is
* recommended for best quality. sampleRate takes the sample rate for the signal 
* in unit Hz, ie. 44100 for 44.1 kHz audio. The data passed to the routine in 
* indata[] should be in the range [-1.0, 1.0), which is also the output range 
* for the data, make sure you scale the data accordingly (for 16bit signed integers
* you would have to divide (and multiply) by 32768). 
*
* COPYRIGHT 1999-2006 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
*
* 						The Wide Open License (WOL)
*
* Permission to use, copy, modify, distribute and sell this software and its
* documentation for any purpose is hereby granted without fee, provided that
* the above copyright notice and this license appear in all source copies. 
* THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
* ANY KIND. See http://www.dspguru.com/wol.htm for more information.
*
*****************************************************************************/

/****************************************************************************
*
* This code was converted to C# by Michael Knight
* madmik3 at gmail dot com. 
* http://sites.google.com/site/mikescoderama/ 
* 
*****************************************************************************/


using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public static class PitchShifter
{

    #region Private Static Memebers
    private static int MAX_FRAME_LENGTH = 16000;
    private static float[] gInFIFO = new float[MAX_FRAME_LENGTH];
    private static float[] gOutFIFO = new float[MAX_FRAME_LENGTH];
    private static float[] gFFTworksp = new float[2 * MAX_FRAME_LENGTH];
    private static float[] gLastPhase = new float[MAX_FRAME_LENGTH / 2 + 1];
    private static float[] gSumPhase = new float[MAX_FRAME_LENGTH / 2 + 1];
    private static float[] gOutputAccum = new float[2 * MAX_FRAME_LENGTH];
    private static float[] gAnaFreq = new float[MAX_FRAME_LENGTH];
    private static float[] gAnaMagn = new float[MAX_FRAME_LENGTH];
    private static float[] gSynFreq = new float[MAX_FRAME_LENGTH];
    private static float[] gSynMagn = new float[MAX_FRAME_LENGTH];
    private static long gRover, gInit;

    public static int Progress = 0;
    #endregion

    #region Public Static  Methods
    public static float[] PitchShift(float pitchShift, long numSampsToProcess,
       float sampleRate, float[] indata)
    {
        return PitchShift(pitchShift, numSampsToProcess, (long)2048, (long)10, sampleRate, indata);
    }
    public static float[] PitchShift(float pitchShift, long numSampsToProcess, long fftFrameSize,
        long osamp, float sampleRate, float[] indata)
    {
        int MAX_FRAME_LENGTH = 16000;
        gInFIFO = new float[MAX_FRAME_LENGTH];
        gOutFIFO = new float[MAX_FRAME_LENGTH];
        gFFTworksp = new float[2 * MAX_FRAME_LENGTH];
        gLastPhase = new float[MAX_FRAME_LENGTH / 2 + 1];
        gSumPhase = new float[MAX_FRAME_LENGTH / 2 + 1];
        gOutputAccum = new float[2 * MAX_FRAME_LENGTH];
        gAnaFreq = new float[MAX_FRAME_LENGTH];
        gAnaMagn = new float[MAX_FRAME_LENGTH];
        gSynFreq = new float[MAX_FRAME_LENGTH];
        gSynMagn = new float[MAX_FRAME_LENGTH];
        gRover = 0; gInit=0;
        float magn, phase, tmp, window, real, imag;
        float freqPerBin, expct;
        long i, k, qpd, index, inFifoLatency, stepSize, fftFrameSize2;


        float[] outdata = indata;
        /* set up some handy variables */
        fftFrameSize2 = fftFrameSize / 2;
        stepSize = fftFrameSize / osamp;
        freqPerBin = sampleRate / (float)fftFrameSize;
        expct = 2.0f * Mathf.PI * (float)stepSize / (float)fftFrameSize;
        inFifoLatency = fftFrameSize - stepSize;
        if (gRover == 0) gRover = inFifoLatency;


        /* main processing loop */
        Debug.Log(numSampsToProcess);
        for (i = 0; i < numSampsToProcess; i++)
        {
            Progress = (int)(((i + 1) * 100) / (1.0f * numSampsToProcess));
            /* As long as we have not yet collected enough data just read in */
            gInFIFO[gRover] = indata[i];
            outdata[i] = gOutFIFO[gRover - inFifoLatency];
            gRover++;

            /* now we have enough data for processing */
            if (gRover >= fftFrameSize)
            {
                gRover = inFifoLatency;
                
                /* do windowing and re,im interleave */
                for (k = 0; k < fftFrameSize; k++)
                {
                    window = -.5f * Mathf.Cos(2.0f * Mathf.PI * (float)k / (float)fftFrameSize) + .5f;
                    gFFTworksp[2 * k] = (float)(gInFIFO[k] * window);
                    gFFTworksp[2 * k + 1] = 0.0F;
                }
                
                /* ***************** ANALYSIS ******************* */
                /* do transform */
                ShortTimeFourierTransform(gFFTworksp, fftFrameSize, -1);
                
                /* this is the analysis step */
                for (k = 0; k <= fftFrameSize2; k++)
                {

                    /* de-interlace FFT buffer */
                    real = gFFTworksp[2 * k];
                    imag = gFFTworksp[2 * k + 1];

                    /* compute magnitude and phase */
                    magn = 2.0f * Mathf.Sqrt(real * real + imag * imag);
                    phase = Mathf.Atan2(imag, real);

                    /* compute phase difference */
                    tmp = phase - gLastPhase[k];
                    gLastPhase[k] = (float)phase;

                    /* subtract expected phase difference */
                    tmp -= (float)k * expct;

                    /* map delta phase into +/- Pi interval */
                    qpd = (long)(tmp / Mathf.PI);
                    if (qpd >= 0) qpd += qpd & 1;
                    else qpd -= qpd & 1;
                    tmp -= Mathf.PI * (float)qpd;

                    /* get deviation from bin frequency from the +/- Pi interval */
                    tmp = osamp * tmp / (2.0f * Mathf.PI);

                    /* compute the k-th partials' true frequency */
                    tmp = (float)k * freqPerBin + tmp * freqPerBin;

                    /* store magnitude and true frequency in analysis arrays */
                    gAnaMagn[k] = (float)magn;
                    gAnaFreq[k] = (float)tmp;

                }

                /* ***************** PROCESSING ******************* */
                /* this does the actual pitch shifting */
                for (int zero = 0; zero < fftFrameSize; zero++)
                {
                    gSynMagn[zero] = 0;
                    gSynFreq[zero] = 0;
                }

                for (k = 0; k <= fftFrameSize2; k++)
                {
                    index = (long)(k * pitchShift);
                    if (index <= fftFrameSize2)
                    {
                        gSynMagn[index] += gAnaMagn[k];
                        gSynFreq[index] = gAnaFreq[k] * pitchShift;
                    }
                    else break;
                }


                /* ***************** SYNTHESIS ******************* */
                /* this is the synthesis step */
                for (k = 0; k <= fftFrameSize2; k++)
                {

                    /* get magnitude and true frequency from synthesis arrays */
                    magn = gSynMagn[k];
                    tmp = gSynFreq[k];

                    /* subtract bin mid frequency */
                    tmp -= (float)k * freqPerBin;

                    /* get bin deviation from freq deviation */
                    tmp /= freqPerBin;

                    /* take osamp into account */
                    tmp = 2.0f * Mathf.PI * tmp / osamp;

                    /* add the overlap phase advance back in */
                    tmp += (float)k * expct;

                    /* accumulate delta phase to get bin phase */
                    gSumPhase[k] += (float)tmp;
                    phase = gSumPhase[k];

                    /* get real and imag part and re-interleave */
                    gFFTworksp[2 * k] = (float)(magn * Mathf.Cos(phase));
                    gFFTworksp[2 * k + 1] = (float)(magn * Mathf.Sin(phase));
                }

                /* zero negative frequencies */
                for (k = fftFrameSize + 2; k < 2 * fftFrameSize; k++) gFFTworksp[k] = 0.0F;

                /* do inverse transform */
                ShortTimeFourierTransform(gFFTworksp, fftFrameSize, 1);

                /* do windowing and add to output accumulator */
                for (k = 0; k < fftFrameSize; k++)
                {
                    window = -.5f * Mathf.Cos(2.0f * Mathf.PI * (float)k / (float)fftFrameSize) + .5f;
                    gOutputAccum[k] += (float)(2.0 * window * gFFTworksp[2 * k] / (fftFrameSize2 * osamp));
                }
                for (k = 0; k < stepSize; k++) gOutFIFO[k] = gOutputAccum[k];

                /* shift accumulator */
                //memmove(gOutputAccum, gOutputAccum + stepSize, fftFrameSize * sizeof(float));
                for (k = 0; k < fftFrameSize; k++)
                {
                    gOutputAccum[k] = gOutputAccum[k + stepSize];
                }

                /* move input FIFO */
                for (k = 0; k < inFifoLatency; k++) gInFIFO[k] = gInFIFO[k + stepSize];
            }
        }

        return outdata;
    }
    #endregion

    #region Private Static Methods
    public static void ShortTimeFourierTransform(float[] fftBuffer, long fftFrameSize, long sign)
    {

        float wr, wi, arg, temp;
        float tr, ti, ur, ui;
        long i, bitm, j, le, le2, k;

        for (i = 2; i < 2 * fftFrameSize - 2; i += 2)
        {
            for (bitm = 2, j = 0; bitm < 2 * fftFrameSize; bitm <<= 1)
            {
                if ((i & bitm) != 0) j++;
                j <<= 1;
            }
            if (i < j)
            {
                temp = fftBuffer[i];
                fftBuffer[i] = fftBuffer[j];
                fftBuffer[j] = temp;
                temp = fftBuffer[i + 1];
                fftBuffer[i + 1] = fftBuffer[j + 1];
                fftBuffer[j + 1] = temp;
            }
        }
        long max = (long)(Mathf.Log(fftFrameSize) / Mathf.Log(2.0f) + .5);
        for (k = 0, le = 2; k < max; k++)
        {
            le <<= 1;
            le2 = le >> 1;
            ur = 1.0F;
            ui = 0.0F;
            arg = (float)Mathf.PI / (le2 >> 1);
            wr = (float)Mathf.Cos(arg);
            wi = (float)(sign * Mathf.Sin(arg));
            for (j = 0; j < le2; j += 2)
            {

                for (i = j; i < 2 * fftFrameSize; i += le)
                {
                    tr = fftBuffer[i + le2] * ur - fftBuffer[i + le2 + 1] * ui;
                    ti = fftBuffer[i + le2] * ui + fftBuffer[i + le2 + 1] * ur;
                    fftBuffer[i + le2] = fftBuffer[i] - tr;
                    fftBuffer[i + le2 + 1] = fftBuffer[i + 1] - ti;
                    fftBuffer[i] += tr;
                    fftBuffer[i + 1] += ti;

                }
                tr = ur * wr - ui * wi;
                ui = ur * wi + ui * wr;
                ur = tr;
            }
        }
    }
    #endregion
}

Full text and comments »

By LashaBukhnikashvili, history, 8 years ago, In English

Here is given n(n=1000) points on the plane, each described by (x,y). need to place a rectangle of area S on the plane. Sides of a rectangle should be paralel to X and Y axis. required to find maximal number of points that can be contained by a single rectangle.

I need to find maximum O(N^2 * K) solution to solve problem, where K is answer (maximum amount of points among all rectangle with area S),

thanks in advance,

Full text and comments »

By LashaBukhnikashvili, 8 years ago, In English

Consider case when we have disjoint polygons on the plane,and they are represented as obstacles,what is the fastest method to calculate distance between 2 point? that's the problem statement...

Let's take variable N and declare it as the total number of vertices on this plane...

I was thinking about some approaches and found this: PATH BETWEEN THESE POINTS GOES ON VERTICES belongs to POLYGONS and GOES ON THIS POINTS ITSELF,that's all... because of that fact I found first approach : make a graph with N+2 points,with expected N^2 edges and run djixstra, complexity of this approach is O(N^2 log N) with memory O(N^2)

But then I realise that it was not enough fast, I started searching and found this material,where(after 11 page in pdf) is described algorithm by [• John Hershberger and Subhash Suri, 1999]: basic idea which works int [time=O(N log N),memory=O(N log N)] time. Here are mentioned dividing plane into 4 part recursively,untill we won't reach single points on plane,like this:

but after that,I didn't get what is explained, if someone has touch/has understood how that algorithm works,please help me to understand how it really works,I want to implement it,thanks in advance :)

Full text and comments »

By LashaBukhnikashvili, 8 years ago, In English

Tomorrow on Feb/07/2016 13:05UTC+4 these rounds will be overlapped to each other,I'm interested if the authors of this Rounds know that fact, and if it is possible make changes about round start time :\

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

given an array,consist of N=10^5 elements.U need to find LIS after each update operation.Update means that u should change some position(1 number) by the some another number.

P.S. I have not any link about problem,it's just my designed task.

Full text and comments »

By LashaBukhnikashvili, history, 9 years ago, In English

I'm trying to solve this problem and I am using Trie.my code.

Here maximum number nodes of Trie is less than 25*1000000,But for each node I have array of 26 elements,and memory is growing up to 25*25*1000000 in worst case.I saw other codes and they are using 2 links,can't understand it,if you can,help me to understand idea of such Trie with compressed array.

P.S to use vector for each node, instead of array, then I will lost time for Build Trie(for finding some character of some node) and there should be *26 TL,and I don't want this

solved my oversight :)

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

Nice to see 3 Div1 contest in the "Current or upcoming contests" bar ^_^

But Div1 Round without div2? it's something new ^_^ or typo

fixed: 308 CF round is for Div2

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

I know you are tired coz of kquery's problems.I am just interested to solve this problem.without update I had on each vertex of segment tree sorted subarray,but now I can't do this with same idea (it was O(n log^2)). also is any idea which solves this problem in O(n log n) ?

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

Help me to prove,that there are infinitely many prime P,for which 2*P+1 also prime.(In fact,I have no idea is it correct or not?)

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

Help me to prove: C(2*n,n) | lcm(1,2,3,...,2*n) for all natural n

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

Need advice about MIT lectures,few days ago I watched 1 of this lecture(about static trees: RMQ,LCA,LA in O(1)),it was great explanation and I've learned,I am interested: If I will watch all this lectures,(how) it will help me?

Full text and comments »

By LashaBukhnikashvili, 9 years ago, In English

I want to prove that for every 200 numbers,it can be choose 100 number which sum will divisible by 100.

I also thinks that,it will be correct: for every n numbers,we can choose x number which sum will divisible by x,if x<=n/2+n%2

can someone help me to prove them(of course,if 2th will be correct)? (sorry for my bad english)

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

I am trying to solve this problem,I need your helps ;)

Is any idea to solve it,without 2d Trees? or what kind of 2d Trees should I use,that dont use O(n*m) memory,and update,query as much fast as it possible(like a logNlogM)

P.S I was reading material about fenwik tree,but it needs O(n*m) memory(sorry for my poor english )

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

Where can I find editorial of the IOIs problems? thanks in advance

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

how prove that sqrt(2)+sqrt(3)+sqrt(5)+sqrt(7)+sqrt(11) is irrational number? thanks in advance

Full text and comments »

  • Vote: I like it
  • +31
  • Vote: I do not like it

By LashaBukhnikashvili, 10 years ago, In English

I need problems about Queries on Tree.

which data structures need problems same as above?

Can you guys suggest problems?

Any help will be appreciated very much :)

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

here is the link of a problem,I wrote O(n log n + m log^3 n) solution,but I have wa16, here is my code.

please help to find my mistake,also I am interested log n and log^2 n solutions.

Any help will be appreciated very much :)

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

I'm trying to solve this porlbem,where you are interested the number of distinct elements in the subsequences.

I remember it was discussed,on codeforces but can't find it(not working "search by tag"),also I was looking in the Internet,but without result.please help me and discuss again or give me the links where already discussed. thank you in advance.

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

what is the complexity of Bitwise operations( &, |, ^, ~, <<, >> ) ?

thank you in advance.

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

during the contest i tried to check is there or not: such m numbers in array that sum of this numbers will be x.

i wrote code,but i dont know why it is not correct,so need your help:

 int d[10001][101]={0};
 d[0][0]=1;
 for (i=1;i<=n;i++){
     cin>>x;
     for (j=10000;j>=0;j--)
          for (k=0;k<n;k++) 
              if (d[j][k]) 
                  d[j+x][k+1]=1;
 };

d[i][j]=1 if we can find in array j numbers where sum of this numbers will be i,otherwise d[i][j]=0;

P.S using this,i get WA6

UPDATE: i didn't know problem good. but this idea is correct,more formally problems same as above can solve in O(s*n*n) where s is the sum of numbers(all numbers positive).

Full text and comments »

By LashaBukhnikashvili, 10 years ago, In English

Can someone help me to deduce this?

For any numbers x,y,z:

UPDATE: problem is wrong see coment.

Full text and comments »

By LashaBukhnikashvili, 11 years ago, In English

The USACO 2013 January contest is available January 11 through January 14.Good Luck to everyone!!!

We can discuss problems here,**ONLY** after the contest is complete.

UPD: Link to Contest Page.You can appear in the maximum of any of the 3 hour window during the contest.Results will be posted within a day or two after the end conest.

UPD: Results

Full text and comments »

By LashaBukhnikashvili, 11 years ago, In English

Hi to everybody....

Today will be another round of COCI.

Those wishing to participate can enter/register here.

To LogIn choose COCI 2012/2013 — Contest #3.

Good luck to all in advance...

UPD: Here are results.

Full text and comments »

By LashaBukhnikashvili, 11 years ago, In English

The USACO 2012 November contest is available November 16 through November 19.Good Luck to everyone!!!

UPD: Contest complete.Results will be posted within a day or two.

UPD: Results here.

Full text and comments »