java - 在java中实现谐波乘积谱算法

Question

我目前正在使用 Java 编写吉他调音器程序，并且我正在尝试实现谐波乘积频谱算法以确定基频。

目前，我采用了一种方法，将我的频谱下采样 f。

我现在正尝试将所有不同的下采样频谱相乘。我正在用 Java 编码并使用数组。因此，我有用于已下采样的索引的数组和用于我的频谱的不同值的数组。我现在正在尝试将所有不同的数组设置为相同的大小并组织它们的值，以便它们对应于正确的下采样索引。我在大小和值方面有很多问题....

此外，我正在尝试根据我在纸上的一个例子来实现这个算法......因此我只能用 4 个下采样频谱来实现这个算法，但我怀疑当我实际使用一个真实的声音信号。

这是下采样方法的代码:

import org.jfree.ui.RefineryUtilities;

public class SousEchantillonnage {

public static double[] valuesDownSample(double[] tab, int factor){

    int N = tab.length;

    double[] values = new double[N];

    for (int i = 0; i < N; i++){
        values[i] = tab[i];
    }

    int lengthDownSample = N + (facteur - 1) * (N - 1);

    double[] valuesDownSample = new double[lengthDownSample];
    for (int i = 0; i < N; i++){
        valuesDownSample[i] = values[i];
    }
    for (int i = N; i < lengthDownSample; i ++){
        valuesDownSample[i] = 0;
    }

    return valuesDownSample;
}

public static double[] indexDownSample(double[] tab, int factor){

    int N = tab.length;

    double[] indexes = new double[N];

    for (int i = 0; i < N; i++){
        indexes[i] = i;
    }

    int lengthDownSample = N + (factor - 1) * (N - 1);

    double[] indexDownSample = new double [lengthDownSample];
    for (int i = 0; i < lengthDownSample; i++){
        indexDownSample[i] = i / factor;
    }

    return indexDownSample;
}

这个方法似乎有效。

到目前为止，这是我的 HPS 算法方法:

public static double[] hps(double[] tab){

    int N = tab.length;

    int factor2 = 2;
    int factor3 = 3;
    int factor4 = 4;
    int lengthDownSample2 = N/2 + (factor2 - 1) * (N/2 - 1);
    int lengthDownSample3 = N/2 + (factor3 - 1) * (N/2 - 1);
    int lengthDownSample4 = N/2 + (factor4 - 1) * (N/2 - 1);

            // Gives us the spectrogram of the signal tab
    double[] spectrogram = new double[N];                       
    spectrogramme = FFT.calculFFT(tab);

            // We only need the first values of the spectrogram. The other half is the same.
    double[] spectrogramCut = new double[N/2];          
    for (int i = 0; i < N/2; i++){
        spectrogramCut[i] = spectrogram[i];
    }

            // We create the array that contains the values of spectrogramCut that we downsample by a factor 2
       double[] valuesSpect2 = new double [sizeDownSamp2];  
       valuesSpect2 = SousEchantillonnage.valuesDownSample(spectrogramCut, factor2);

          // We create an array of the indexes of spectrogramCut that we downsample by a factor 2
      double[] indexSpect2 = new double[sizeDownSamp2];
      indexSpect2 = SousEchantillonnage.indexDownSample(spectrogramCut, factor2);

         // We create the array that contains the values of spectrogramCut that we downsample by a factor 3
       double[] valuesSpect3 = new double [sizeDownSamp3];  
       valuesSpect3 = SousEchantillonnage.valuesDownSample(spectrogramCut, factor3);

         // We create an array of the indexes of spectrogramCut that we downsample by a factor 3
    double[] indexSpect3 = new double[sizeDownSamp3];
    indexSpect3 = SousEchantillonnage.indexDownSample(spectrogramCut, factor3);;

       // We create the array that contains the values of spectrogramCut that we            downsample by a factor 4
   double[] valuesSpect4 = new double [sizeDownSamp4];  
   valuesSpect4 = SousEchantillonnage.valuesDownSample(spectrogramCut, factor4);

       // We create an array of the indexes of spectrogramCut that we downsample by a factor 4
    double[] indexSpect4 = new double[sizeDownSamp4];
    indexSpect4 = SousEchantillonnage.indexDownSample(spectrogramCut, factor4);

        int sizeIndex = N/2 + 5 * (N/2 - 1); // size of the array that contains all the       indexes of the downsamples

        // We create this array
    double[] indexDowSamp = new double[sizeIndex];
    indexDowSamp[0] = 0;
    indexDowSamp[1] = 0.25;
    indexDowSamp[2] = 0.333;
    indexDowSamp[3] = 0.5;
    indexDowSamp[4] = 0.666;
    indexDowSamp[5] = 0.75;

    int q = sizeIndex / 6;      // quantity of packets of 6 we can do
    int r = sizeIndex%6;        // what we are left with.

    for (int i = 6; i < q * 6; i += 6){
        for (int j = 0; j < 6; j++){
        indexDowSamp[i + j] = indexDowSamp[i + j - 6] + 1;
        }
    }
    for (int i = 0; i < r; i++){
        indexDowSamp[q * 6 + i] = indexDowSamp[q * 6 + i - 6] + 1;
    }

我被困在这里了。我想做一个将两个不同长度的数组相乘的方法。

基本上，如您所见，当我对频谱图进行下采样时，我得到两个数组:

• 一个具有缩减采样的索引
• 另一个具有下采样后的值。

我想做的是创建一个与 indexDownSample: valuesDownSample 大小相同的数组。例如，我们有 indexDownSample[0] = 0。我想为 valuesDownSample[0] 使用 valuesSpectCut[0] *valuesSpect2[0]*valuesSpect3[0]*valuesSpect4[0] 的乘积，因为所有这些数组具有对应于索引 0 的值(indexSpectCut[0] = 0，indexSpect2[0] = 0 = indexSpect3[0] = indexSpect4[0])。

对于 indexDownSample[1]=0.25，我们注意到只有 indexSpect4[1]= indexDownSample[1] = 0.25然后，我们将为 valuesDownSample[1] 默认设置 0。

然后我们继续这样，直到我们填满数组。

如果一切顺利，我们应该在最后:

• 包含产品不同值的 valuesDownSample
• 包含不同下采样索引的 indexDownSample。

我只需要找到最大峰值即可找到我的基频。

我唯一的问题是我不知道如何做乘法......

如果有人有想法，我将不胜感激!

Answer 1

好吧，我有一个实际的答案，这就是我解决问题的方法:

首先，我不称其为下采样，而是将 FFT 的频率分开。为此，我只使用 2 个数组，如下所示:

float findex[1000]; //index of frequencies for harmonic product spectrum
float mindex[1000]; //index of magnitudes for harmonic product spectrum
unsigned int max_findex;    //number of elements in findex[] and mindex[]

这是我的初始化: max_findex = 0; 对于 (i = 0; i < 1000; i++) { mindex[i] = 0.0;

我正在使用 VisualDSP++ 在 Analog Devices Sharc EZ 开发板上进行开发，所以我的代码是用 C 编写的。

我在 12kHz 的低采样率下使用 131072 的大 FFT 大小。这为我提供了每个 bin 0.09Hz 的相对较高的精度。由于高次谐波会更准确，所以我首先从最高分开始:

//first run with max division (highest accuracy)
for (i = 0; i < new_peak; i++) {
    findex[max_findex] = f(peak[i].index) / 9;
    mindex[max_findex] = peak[i].magnitude;
    if (max_findex < 999) max_findex++;
    else xmitUARTmessage("ERROR max_findex\r\n", 100);
}

我的所有 FFT 峰值都在这样的结构中:peak[].index 是 FFT 中的 bin 编号，peak[].magnitude 是峰值的幅度。 f() 函数返回 bin 的频率。

接下来我会转到第 8 部分，然后是第 7 部分，依此类推，最后是第 1 部分(实际频率排在最后)。对于每个被分割的峰，我查看我的阵列，看看此时是否已经有一个频率。我使用 +/- 0.2，我可能会收紧它，但你必须调整它以适应你的 FFT 的精度。

char found;

for (u2 = 8; u2 > 0; u2--) {
    for (i = 0; i < new_peak; i++) {
        tempf = f(peak[i].index) / u2;
        found = 0;
        for (u = 0; u < max_findex; u++) {
            //try to find existing frequency
            if (tempf < findex[u] + 0.2 && tempf > findex[u] - 0.2) {
                //found existing frequency
                mindex[u] *= peak[i].magnitude;
                found = 1;
                break;
            }
        } //for u
        if (!found) {
            //no frequency was found, add new one
            findex[max_findex] = tempf;
            mindex[max_findex] = peak[i].magnitude;
            if (max_findex < 999) max_findex++;
            else xmitUARTmessage("ERROR max_findex\r\n", 100);
        }
    } //for i
} //for u2

就是这样。现在我只需打印出我的值并在 Excel 中按大小对它们进行排序...

for (i = 0; i < max_findex; i++) {
    sprintf(tempstr, "%.2f,%.2f\r\n", findex[i], mindex[i]);
    xmitUARTmessage(tempstr, 100);
}

这是一些输出(我只显示前 6 个):60.53 170569325000060 155841987500020 555159950179.99 2649815257.5 13173538.57 1317353

创建输出的输入音频信号是:60 Hz 方波(基波加奇次谐波)181.5 赫兹正弦波302.5 赫兹正弦波423.5赫兹正弦波

我正在模拟 60 + 60.5 Hz 的拍频，基本频率在 60.5 时完全缺失。这是一个棘手的案例，但它奏效了:)