怎么知道在谐波乘积频谱中下采样了多少次？

信息处理 fft 频率 信号分析
2022-02-18 22:18:09

我一直在研究谐波乘积谱算法。现在，我读过的所有关于该主题的文献都告诉我要下采样 N 次。如何确定这个 N 值应该是多少？到目前为止，这是我对谐波产品频谱的实现。如果我在某个地方出错了，请随时告诉我。

private int HarmonicProductSpectrum(Complex[] fftData, int n){
    Complex[][] data = new Complex[n][fftData.length/n];
    for(int i = 0; i<n; i++){
        for(int j = 0; j<data[0].length; j++){
            data[i][j] = fftData[j*(i+1)];
        }
    }
    Complex[] result = new Complex[fftData.length/n];//Combines the arrays
    for(int i = 0; i<result.length; i++){
        Complex tmp = new Complex(1,0);
        for(int j = 0; j<n; j++){//multiplies arrays together
            tmp = tmp.times(data[j][i]);
        }
        result[i] = tmp;
    }
    //Calculates Maximum Magnitude of the array
    double max = Double.MIN_VALUE;
    int index = -1;
    for(int i = 0; i<result.length; i++){
        Complex c = result[i];
        double tmp = c.getMagnitude();
        if(tmp>max){
            max = tmp;;
            index = i;
        }
    }
    return index*getFFTBinSize(fftData.length);
}
2个回答

取决于音调源频谱、可能的最低音调和 FFT 长度。

如果 N 太小，算法可能会遗漏一些包含音高频谱能量的重要部分的高次谐波。所以你需要知道在你的特定音高源中有多少泛音可能是重要的。

但是，如果 N 太大，在对频谱进行下采样后，单个低音的多个泛音可能最终会出现在同一个 bin 中，从而混淆结果。

对于具有极其丰富的高次谐波的非常低的音高，这两个约束可能重叠，因此表明 HPS 需要更长的 FFT 窗口，甚至需要完全不同的音高估计方法。

看看第 4 页的这篇论文。它对算法有很好的图形表示。我将尝试在这里编写一些伪代码，您可以根据需要适应您的语言（在我看来就像 Java）。

/* Harmonic Product spectrum PSEUDOCODE
 * @param fftData Discrete Fourier transform of data
 * @param N Number of times we downsample the spectrum to get HPS
 */
int HPS( Complex[] fftData, int N )
{
    // Find magnitude of the FFT
    Real fullSpectrum[] = absOfComplex(fftData);

    // Keep only the positive frequencies (DC to Nyquist)
    Real spectrum[] = discardNegativeFrequencies(fullSpectrum);

    // Make a new array to store HPS
    Real hps[] = copyOf(spectrum);

    // Perfrom HPS:
    // Go through each downsampling factor
    for (int downsamplingFactor = 1; downsamplingFactor <= N; downsamplingFactor++)
    {
        // Go through samples of the downsampled signal and compute HPS at this iteration
        for(int idx = 0; idx < spectrum.length()/downsamplingFactor; idx++)
        {
            hps[idx] *= spectrum[idx * downsamplingFactor];
        }
    }

    return findIndexOfMax(hps);
}

这不是防错的，但它应该让你开始。

其它你可能感兴趣的问题
上一篇使用 Karplus Strong 合成泛音 下一篇复合信号的带宽