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(f0floor∗0.9f0ceil∗1.1)/ln2∗channelsinoctavenumberofchannels = 1 + log(\frac {f0floor*0.9} {f0ceil*1.1}) /ln2 * channelsinoctave numberofchannels=1+log(f0ceil∗1.1f0floor∗0.9)/ln2∗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);// 创建一个 NuttallWindowdouble *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);

卷积

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

// Convolutiondouble 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，后续的频率会逐渐增加，这里的柱状条会逐渐变胖

【Harvest源码分析】GetFilteredSignal函数相关推荐

jQuery源码分析-each函数
本文部分截取自且行且思 jQuery.each方法用于遍历一个数组或对象,并对当前遍历的元素进行处理,在jQuery使用的频率非常大,下面就这个函数做了详细讲解: 复制代码代码 /*! * jQuer ...
x265源码分析 main函数 x265.cpp
图片转载于x265源码流程分析_Dillon2015的博客-CSDN博客_x265编码流程 cliopt.prase main ()函数--解析函数参数并进行编码准备工作:x265.cpp (1)Ge ...
GCC源码分析(十) — 函数节点的gimple高端化
版权声明:本文为CSDN博主「ashimida@」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明. 原文链接:https://blog.csdn.net/lidan1 ...
GCC源码分析(十一) — 函数节点的gimple低端化
版权声明:本文为CSDN博主「ashimida@」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明. 原文链接:https://blog.csdn.net/lidan1 ...
【Linux 内核内存管理】mmap 系统调用源码分析 ④ ( do_mmap 函数执行流程 | do_mmap 函数源码 )
文章目录一.do_mmap 函数执行流程二.do_mmap 函数源码调用 mmap 系统调用 , 先检查 " 偏移 " 是否是 " 内存页大小 " 的 & ...
c++imread 函数_OpenCV学习笔记与源码分析: imread( )函数
引言 imread()函数在opencv使用比较. imread()函数声明: Mat imread(const string& filename, int flags); 这很标准的写法, ...
【Harvest源码分析】GetF0CandidateContour函数
该函数,顾名思义,就是描绘F0 Candidate的轮廓. 在GetFourZeroCrossingIntervals函数中,已经介绍了如何获取ZeroCrossings了.不清楚的可以再看一下我的那 ...
【Harvest源码分析】GetFourZeroCrossingIntervals函数
背景一个完整的正弦波存在如下性质, 波峰间隔,波谷间隔,向上过零间隔,向下过零间隔这四者的值理论上应该一致的. 那么该函数,顾名思义,就是取这四段值的代码数据结构用来保存这四段值得结构 Zer ...
【Harvest源码分析】GetWaveformAndSpectrumSub函数
GetWaveformAndSpectrumSub函数基本就是Matlab版本的decimate实现.用于降采样. Matlab版本的decimate使用在这里在输入信号x的首尾,分别加了140长度 ...

【Harvest源码分析】GetFilteredSignal函数

入参