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数组。
n u m b e r o f c h a n n e l s = 1 + l o g ( f 0 f l o o r ∗ 0.9 f 0 c e i l ∗ 1.1 ) / l n 2 ∗ c h a n n e l s i n o c t a v e numberofchannels = 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);
// 创建一个 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卷积后的频谱
有意思的是,这个图很像潮水从左边冲向右边,原始信号频率(假设为固体)在其中冲刷出来的痕迹
卷积后的时域谱
具体的变化过程,我在这里有记录。