フレーム間の位相変化を使用して、FFTビンから正確な周波数を抽出する

https://stackoverflow.com/questions/4633203

08-10-2019
|

質問

私はこの素晴らしい記事を見てきました： http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

素晴らしく、それは非常に硬くて重いです。この素材は本当に私を伸ばしています。

特定のビンの正確な周波数を計算するStefanのコードモジュールから数学を抽出しました。しかし、私は最後の計算を理解していません。誰かが最後に数学的な構造について説明できますか？

コードを掘り下げる前に、シーンを設定させてください。

fftframesize = 1024を設定したので、512+1ビンを扱っています
例として、Bin [1]の理想的な周波数は、フレーム内の単一の波に適合します。 40kHzのサンプルレートでは、ToneFrame = 1024/40K秒= 1/40sなので、Bin [1]は理想的には40Hzの信号を収集します。
OSAMP（オーバーサンプル）= 4を設定すると、256のステップで入力信号に沿って進行します。したがって、最初の分析では、バイトゼロから1023、次に256〜1279などを調べます。

...

void calcBins( 
              long fftFrameSize, 
              long osamp, 
              float sampleRate, 
              float * floats, 
              BIN * bins
              )
{
    /* initialize our static arrays */
    static float gFFTworksp[2*MAX_FRAME_LENGTH];
    static float gLastPhase[MAX_FRAME_LENGTH/2+1];

    static long gInit = 0;
    if (! gInit) 
    {
        memset(gFFTworksp, 0, 2*MAX_FRAME_LENGTH*sizeof(float));
        memset(gLastPhase, 0, (MAX_FRAME_LENGTH/2+1)*sizeof(float));
        gInit = 1;
    }

    /* do windowing and re,im interleave */
    for (long k = 0; k < fftFrameSize; k++) 
    {
        double window = -.5*cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5;
        gFFTworksp[2*k] = floats[k] * window;
        printf("sinValue: %f", gFFTworksp[2*k]);
        gFFTworksp[2*k+1] = 0.;
    }

    /* do transform */
    smbFft(gFFTworksp, fftFrameSize, -1);

    printf("\n");

    /* this is the analysis step */
    for (long k = 0; k <= fftFrameSize/2; k++) 
    {
        /* de-interlace FFT buffer */
        double real = gFFTworksp[2*k];
        double imag = gFFTworksp[2*k+1];

        /* compute magnitude and phase */
        double magn = 2.*sqrt(real*real + imag*imag);
        double phase = atan2(imag,real);

        /* compute phase difference */
        double phaseDiff = phase - gLastPhase[k];
        gLastPhase[k] = phase;

        /* subtract expected phase difference */
        double binPhaseOffset = M_TWOPI * (double)k / (double)osamp;
        double deltaPhase = phaseDiff - binPhaseOffset;

        /* map delta phase into [-Pi, Pi) interval */
        // better, but obfuscatory...
        //    deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);

        while (deltaPhase >= M_PI)
            deltaPhase -= M_TWOPI;
        while (deltaPhase < -M_PI)
            deltaPhase += M_TWOPI;

（編集:)今、私が取得しないビット：

        // Get deviation from bin frequency from the +/- Pi interval 
        // Compute the k-th partials' true frequency    

        // Start with bin's ideal frequency
        double bin0Freq = (double)sampleRate / (double)fftFrameSize;
        bins[k].idealFreq = (double)k * bin0Freq;

        // Add deltaFreq
        double sampleTime = 1. / (double)sampleRate;
        double samplesInStep = (double)fftFrameSize / (double)osamp;
        double stepTime = sampleTime * samplesInStep;
        double deltaTime = stepTime;        

        // Definition of frequency is rate of change of phase, i.e. f = dϕ/dt
        // double deltaPhaseUnit = deltaPhase / M_TWOPI; // range [-.5, .5)
        double freqAdjust = (1. / M_TWOPI) * deltaPhase / deltaTime; 

        // Actual freq <-- WHY ???
        bins[k].freq = bins[k].idealFreq + freqAdjust;
    }
}

顔を見つめているように見えますが、はっきりと見ることができません。誰かがこのプロセスをゼロから、段階的に説明してもらえますか？

解決 4

最後に、私はこれを理解しました。本当に私はそれをゼロから導き出さなければなりませんでした。私はそれを導き出すためのいくつかの簡単な方法があると知っていました、私の（通常の）間違いは、私自身の常識を使用するのではなく、他の人の論理に従うことを試みることでした。

このパズルがかかります 2つのキー ロックを解除します。

最初のキー 過剰サンプリングがビン相でどのように回転を導入するかを理解することです。
2番目のキー ここでグラフ3.3と3.4から来ています： http://www.dspdimension.com/admin/pitch-shifting-using-the-ft/

...

for (int k = 0; k <= fftFrameSize/2; k++) 
{
    // compute magnitude and phase 
    bins[k].mag = 2.*sqrt(fftBins[k].real*fftBins[k].real + fftBins[k].imag*fftBins[k].imag);
    bins[k].phase = atan2(fftBins[k].imag, fftBins[k].real);

    // Compute phase difference Δϕ fo bin[k]
    double deltaPhase;
    {
        double measuredPhaseDiff = bins[k].phase - gLastPhase[k];
        gLastPhase[k] = bins[k].phase;

        // Subtract expected phase difference <-- FIRST KEY
        // Think of a single wave in a 1024 float frame, with osamp = 4
        //   if the first sample catches it at phase = 0, the next will 
        //   catch it at pi/2 ie 1/4 * 2pi
        double binPhaseExpectedDiscrepancy = M_TWOPI * (double)k / (double)osamp;
        deltaPhase = measuredPhaseDiff - binPhaseExpectedDiscrepancy;

        // Wrap delta phase into [-Pi, Pi) interval 
        deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);
    }

    // say sampleRate = 40K samps/sec, fftFrameSize = 1024 samps in FFT giving bin[0] thru bin[512]
    // then bin[1] holds one whole wave in the frame, ie 44 waves in 1s ie 44Hz ie sampleRate / fftFrameSize
    double bin0Freq = (double)sampleRate / (double)fftFrameSize;
    bins[k].idealFreq = (double)k * bin0Freq;

    // Consider Δϕ for bin[k] between hops.
    // write as 2π / m.
    // so after m hops, Δϕ = 2π, ie 1 extra cycle has occurred   <-- SECOND KEY
    double m = M_TWOPI / deltaPhase;

    // so, m hops should have bin[k].idealFreq * t_mHops cycles.  plus this extra 1.
    // 
    // bin[k].idealFreq * t_mHops + 1 cycles in t_mHops seconds 
    //   => bins[k].actualFreq = bin[k].idealFreq + 1 / t_mHops
    double tFrame = fftFrameSize / sampleRate;
    double tHop = tFrame / osamp;
    double t_mHops = m * tHop;

    bins[k].freq = bins[k].idealFreq + 1. / t_mHops;
}

他のヒント

基本原則は非常に単純です。特定のコンポーネントがビン周波数と正確に一致する場合、その位相は1つのFTから次のFTに変化しません。ただし、周波数がビン周波数と正確に対応していない場合、連続したFTS間に位相変化があります。周波数deltaは次のようです。

delta_freq = delta_phase / delta_time

そして、コンポーネントの頻度の洗練された推定値は次のとおりです。

freq_est = bin_freq + delta_freq

このアルゴリズムを実装しましたパフォーマンス私自身。タイムオフセットで別のFFTを服用すると、位相がオフセットに応じて変化すると予想されます。つまり、256のサンプルを採取した2つのFFTは、信号に存在するすべての周波数に対して256サンプルの位相差を持つ必要があります（これは、信号自体がそれ自体を想定しています安定していますが、これは256サンプルのような短期間の良い仮定です）。

これで、FFTから得られる実際の位相値はサンプルではなく位相角度であるため、周波数によって異なります。次のコードでは、Phasestep値はビンごとに必要な変換係数です。つまり、ビンXに対応する周波数の場合、位相シフトはx * phasestepになります。 Binセンターの周波数Xの場合Xは整数（ビン数）になりますが、実際に検出された周波数の場合は、任意の実数である可能性があります。

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

補正は、ビン内の信号にビン中心周波数があると仮定し、そのために予想される位相シフトを計算することにより機能します。この予想されるシフトは、実際のシフトから差し引かれ、エラーが残ります。残り（Modulo 2 Pi）が採取され（-PIからPI範囲）、最終周波数はBin Center +補正で計算されます。

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

デルタ補正はいずれにせよ最大0.5 * fft_n / fft_stepビンになる可能性があるため、多くの隣接するビンが同じ周波数に補正されることが多いことに注意してください。不正確さのために必要なものと不正確さ）。

私はこれが役立つことを願っています:)

これは、位相ボコーダー法で使用される周波数推定手法です。

時間内に（固定周波数と固定振幅）正弦波の単一のポイントを見ると、周波数に比例した量だけで位相が前進します。または、逆を行うことができます。副鼻腔の位相が時間単位でどれだけ変化するかを測定すると、そのsinusoidの頻度を計算できます。

位相ボコーダーは、2つのFFTウィンドウを参照して2つのFFTSを使用して位相を推定します。2つのFFTSのオフセットは、時間内の2位測定間の距離です。そこから、そのFFTビンの周波数推定値があります（FFTビンは、正弦波成分またはそのビン内に適合するその他の十分に狭帯域信号を分離するためのほぼフィルターです）。

この方法が機能するためには、使用中のFFTビン近くのスペクトルはかなり静止している必要があります。たとえば、周波数が変化しないなどです。これは、位相ボコーダーが必要とする仮定です。

多分これは助けになるでしょう。 FFTビンは、それぞれがビンの周波数で回転する小さなクロックまたはローターを指定すると考えてください。安定した信号の場合、ローターの（理論的）次の位置は、取得できないビットの数学を使用して予測できます。この「必要のある」（理想的な）位置に対して、いくつかの有用なものを計算できます。（1）隣接するフレームのビンの位相との違いは、 位相ボコーダー ビン周波数をより適切に推定するため、または（2）より一般的に 位相偏差, 、これは、音声のメモ開始またはその他のイベントの肯定的な指標です。

ビン周波数で正確に落ちる信号周波数は、2πの整数倍数によってビン位相を前進させます。ビン周波数に対応するビン相は、FFTの周期的な性質のために2πの倍数であるため、この場合には相変化はありません。あなたが言及する記事もこれを説明しています。

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow