0
$\begingroup$

I am faced with the following code snippet to get an understanding of a F0 estimation algorithm (DIO). Have a quick look below:

static void DesignLowCutFilter(int N, int fft_size, double *low_cut_filter) {
  for (int i = 1; i <= N; ++i)
    low_cut_filter[i - 1] = 0.5 - 0.5 * cos(i * 2.0 * world::kPi / (N + 1));
  for (int i = N; i < fft_size; ++i)
    low_cut_filter[i] = 0.0;
  double sum_of_amplitude = 0.0;
  for (int i = 0; i < N; ++i)
    sum_of_amplitude += low_cut_filter[i];
  for (int i = 0; i < N; ++i)
    low_cut_filter[i] = -low_cut_filter[i] / sum_of_amplitude;
  for (int i = 0; i < (N - 1) / 2; ++i)
    low_cut_filter[fft_size - (N - 1) / 2 + i] = low_cut_filter[i];
  for (int i = 0; i < N; ++i)
    low_cut_filter[i] = low_cut_filter[i + (N - 1) / 2];
  low_cut_filter[0] += 1.0;
}

I struggle to understand how this realises a low-pass filter.

  • The arguments: N seems to be the cut-off frequency in samples, fft_size is the number of bins in our FFT, and low cut filter seems to be a pointer to an array that we should populate with our low-pass filter coefficients
  • The first line seems to be a Hann function with a window length of N.
  • Then a cumulative amplitude is calculated to negate and normalise the filter. I guess this is to ensure all numbers below are below -1.
  • Then we seem to do something which looks like copying the signal to the other end of the spectrum.
  • Then there is another copying part, which I don't get.
  • Then we add one to the first bin.

I realise the main confusion for me is that I would start making a low pass filter in a completely different way and I don't see the connections with the above algorithm.

For an n-th order Butterworth low pass filter, I would make a transfer function:

$$ G^{2}(\omega) = \frac{G_0^2}{1 + (\frac{j \omega}{j \omega_c})^{2n}} $$ I would then multiply this transfer function together with the FFT of my signal, that would be the easiest. I can see that this is not a Butterworth filter, but I don't know what kind of filter the above is.

Could you help me in providing some pointers to find the missing concepts in my understanding?

$\endgroup$
1
  • $\begingroup$ as the name of your function suggests this is not a low pass filter this is a low cut filter. and it's quite opposite of low pass filter. $\endgroup$ – m-sh-shokouhi Aug 6 '20 at 9:48
2
$\begingroup$

This is one tortured piece of code. I strongly recommend to toss it and write whatever you need from scratch.

Arguments

  1. N filter order or some design parameter
  2. fft_size: size of the final vector
  3. low_cut: the actual filter impulse Response padded to FFT size

Code

  1. First loop creates a cosine window of length N and bias it with 0.5
  2. 2nd loop zero pads to FFT length
  3. 3rd loop calcuates the sum of the amplitudes
  4. 4th loop divides by the sum and multiplies with -1.
  5. 5th loop seems to time reverse the filter (pointless since it's symmetrical) AND copies towards the end of FFT buffer AND also writes past the end of the FFT buffer in the process.
  6. 6th loop shift N samples forward by N/2-1 overwriting what was previously in there

What happens exactly here depends on values of N and fft_size and there is some minor low cut behavior in there, but overall it seems like and extremely convoluted, imprecise and dangerous way of doing a low cut

$\endgroup$
1
  • $\begingroup$ I can assure you my intention is to rewrite this code, but first, I need to understand what it does. I believe it doesn't write past the buffers in loop 5 as the last step will be low_culter_filter[fft_size - 1] due to the strict inequality. My question is directed towards what individual steps intend to do, but I realise that this might be a flawed implementation. Could you (1) either point out whats the role of these steps in realising a low-cut/high-pass filter (2) or point me to a similar but less "tortured" implementation? $\endgroup$ – boomkin Aug 6 '20 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.