1

我目前正在使用 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 We 将默认为 0 valuesDownSample[1]。

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

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

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

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

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

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

4

1 回答 1

4

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

首先，我不称之为下采样，而是划分 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 大小。这给了我每箱 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 1705693250000 60 1558419875000 20 555159950 179.99 264981525 7.5 1317353 8.57 1317353

产生输出的输入音频信号是： 60 Hz 方波（基波加奇次谐波） 181.5 Hz 正弦波 302.5 Hz 正弦波 423.5 Hz 正弦波

我正在模拟 60 + 60.5 Hz 的拍频，而基波在 60.5 处完全缺失。这是一个艰难的案例，它奏效了:)

于 2013-11-01T10:46:45.600 回答