【Harvest源码分析】GetFilteredSignal函数

static void GetFilteredSignal(double boundary_f0, int fft_size, double fs,
const fft_complex *y_spectrum, int y_length, double *filtered_signal)

入参

• boundary_f0 边界频率，来自boundary_f0_list。它是一个长度为number_of_channels(185)的double数组。

$$numberofchannels=1+log(\frac{f0floor\ast 0.9}{f0ceil\ast 1.1})/ln2\ast channelsinoctave$$

• fs 降采样后的频率，值为7350
• fft_size 计算fft的长度
• y_spectrum 频率
• y_length 信号长度
• filtered_signal 输出信号

构建NuttallWindow

有关NuttallWindow，请参考这里

  int filter_length_half = matlab_round(fs / boundary_f0 * 2.0);
  // 创建一个 NuttallWindow
  double *band_pass_filter = new double[fft_size];
  NuttallWindow(filter_length_half * 2 + 1, band_pass_filter);

这里NuttallWindow本身的时域谱如下，

它与HanningWindow有何不同呢？
对比一下看一下，它比hanningWindow要“瘦”。相比于hanning Window能量更集中于中部

构建band pass filter

每一个boundary_f0使用NuttallWindow做窗。不够的部分填0.
但是filter_length_half是跟boundary_f0反比关系，所以band pass filter会越来越短。

  // int filter_length_half = matlab_round(fs / boundary_f0 * 2.0);
  for (int i = -filter_length_half; i <= filter_length_half; ++i)
    band_pass_filter[i + filter_length_half] *=  cos(2 * world::kPi * boundary_f0 * i / fs);
  // 做fft后面半部分置零
  for (int i = filter_length_half * 2 + 1; i < fft_size; ++i)
    band_pass_filter[i] = 0.0;

通过图片可以看到band pass filter的变化，是越来越收窄的

做FFT

  fft_complex *band_pass_filter_spectrum = new fft_complex[fft_size];
  fft_plan forwardFFT = fft_plan_dft_r2c_1d(fft_size, band_pass_filter,
    band_pass_filter_spectrum, FFT_ESTIMATE);
  fft_execute(forwardFFT);

卷积

卷积计算公式可以参考这里
这里只提一句话。 频谱的乘积就是时域的卷积。

// Convolution
  double tmp = y_spectrum[0][0] * band_pass_filter_spectrum[0][0] -
    y_spectrum[0][1] * band_pass_filter_spectrum[0][1];
  band_pass_filter_spectrum[0][1] =
    y_spectrum[0][0] * band_pass_filter_spectrum[0][1] +
    y_spectrum[0][1] * band_pass_filter_spectrum[0][0];
  band_pass_filter_spectrum[0][0] = tmp;
  for (int i = 1; i <= fft_size / 2; ++i) {
    tmp = y_spectrum[i][0] * band_pass_filter_spectrum[i][0] -
      y_spectrum[i][1] * band_pass_filter_spectrum[i][1];
    band_pass_filter_spectrum[i][1] =
      y_spectrum[i][0] * band_pass_filter_spectrum[i][1] +
      y_spectrum[i][1] * band_pass_filter_spectrum[i][0];
    band_pass_filter_spectrum[i][0] = tmp;
    band_pass_filter_spectrum[fft_size - i - 1][0] =
      band_pass_filter_spectrum[i][0];
    band_pass_filter_spectrum[fft_size - i - 1][1] =
      band_pass_filter_spectrum[i][1];
  }

逆向FFT

fft_plan inverseFFT = fft_plan_dft_c2r_1d(fft_size,
    band_pass_filter_spectrum, filtered_signal, FFT_ESTIMATE);
  fft_execute(inverseFFT);

补偿

  // Compensation of the delay.
  int index_bias = filter_length_half + 1;
  for (int i = 0; i < y_length; ++i)
    filtered_signal[i] = filtered_signal[i + index_bias];

---

流程可视化

原始信号频谱

band_pass_filter频谱

这里取的是第一个filter，值应该是36.6HZ，后续的频率会逐渐增加，这里的柱状条会逐渐变胖 😃
观测窗的长度变短，会导致频率分辨率下降。这就导致在FFT后能量就越往周围扩散。但又因为能量总和是一致的，这就导致会变“矮胖”。 如果在源码总的窗长度固定的化，这里的频谱将是一个保持身材右移的过程。

这个是330HZ，变胖了 ，也变矮了。

这个是606HZ

为了明显的看出这个变化，我做了个动态图来看一下变化。从又高又瘦变成又矮又平。

卷积后的频谱

下图是原始信号的频谱：

下图为不同频率band pass filter卷积后的频谱
有意思的是，这个图很像潮水从左边冲向右边，原始信号频率（假设为固体）在其中冲刷出来的痕迹

卷积后的时域谱

具体的变化过程，我在这里有记录。