0
$\begingroup$

I'm attempting to interpret the following function. To my understanding, this aims to implement an FIR low-pass filter with a Nuttall window. Then, it filters the signal by a simple application of the convolution equivalence theorem.

static void GetFilteredSignal(int half_average_length, int fft_size, const fft_complex *y_spectrum, int y_length, double *filtered_signal) { double *low_pass_filter = new double[fft_size]; // Nuttall window is used as a low-pass filter. // Cutoff frequency depends on the window length. NuttallWindow(half_average_length * 4, low_pass_filter); for (int i = half_average_length * 4; i < fft_size; ++i) low_pass_filter[i] = 0.0; fft_complex *low_pass_filter_spectrum = new fft_complex[fft_size]; fft_plan forwardFFT = fft_plan_dft_r2c_1d(fft_size, low_pass_filter, low_pass_filter_spectrum, FFT_ESTIMATE); fft_execute(forwardFFT); // Convolution double tmp = y_spectrum[0][0] * low_pass_filter_spectrum[0][0] - y_spectrum[0][1] * low_pass_filter_spectrum[0][1]; low_pass_filter_spectrum[0][1] = y_spectrum[0][0] * low_pass_filter_spectrum[0][1] + y_spectrum[0][1] * low_pass_filter_spectrum[0][0]; low_pass_filter_spectrum[0][0] = tmp; for (int i = 1; i <= fft_size / 2; ++i) { tmp = y_spectrum[i][0] * low_pass_filter_spectrum[i][0] - y_spectrum[i][1] * low_pass_filter_spectrum[i][1]; low_pass_filter_spectrum[i][1] = y_spectrum[i][0] * low_pass_filter_spectrum[i][1] + y_spectrum[i][1] * low_pass_filter_spectrum[i][0]; low_pass_filter_spectrum[i][0] = tmp; low_pass_filter_spectrum[fft_size - i - 1][0] = low_pass_filter_spectrum[i][0]; low_pass_filter_spectrum[fft_size - i - 1][1] = low_pass_filter_spectrum[i][1]; } fft_plan inverseFFT = fft_plan_dft_c2r_1d(fft_size, low_pass_filter_spectrum, filtered_signal, FFT_ESTIMATE); fft_execute(inverseFFT); // Compensation of the delay. int index_bias = half_average_length * 2; for (int i = 0; i < y_length; ++i) filtered_signal[i] = filtered_signal[i + index_bias]; fft_destroy_plan(inverseFFT); fft_destroy_plan(forwardFFT); delete[] low_pass_filter_spectrum; delete[] low_pass_filter; }

There are two major points that I don't understand in this algorithm:

  1. Half average length. What is it? It is used to calculate the filter order and the comment says it influences the cutoff frequency. I realise that the window length limits resolution when applying to a signal. My confusion is that all DSP guides I look up usually assume some filter order N as a design parameter, rather than selecting it based on a desired cut-off frequency. Is there a relationship between N and then cutoff frequency then in this case/generally?
  2. Compensation of the delay. I'm not familiar with the details of filter design, but after reading up, I found that it is essential to shift a time-domain signal after FIR filtering. Filtering usually has frequency-dependent group delay, but FIR filters have the advantage of a constant group delay being $ \frac{N-1}{2} $, where N is the number of taps or order of the filter. Assuming $ \text{half_average_length} \cdot 4 = N $, it follows that $ \text{index_bias} = (\text{half_average_length} \cdot 4 - 1) /2 $ . I realise this wouldn't make sense as it's not a whole number, but then there is something wrong with my interpretation of $ \frac{N-1}{2} $. What is wrong with my interpretation of this formula?
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

Half average length. What is it?

Primarily something the author made up. They design a lowpass more or less as a weighted moving average filter. They simply use the window itself as the impulse response of the filter. That's NOT the same as designing an ideal low pass filter at the cutoff frequency you want and then windowing it down to the filter order. This seems to work reasonably well, you get frequency response of a nuttall window (which isn't bad) and the cut off frequency is roughly the sample rate divided by the window length.

Ir your sample rate is 44.1kHz and your window length is, say, N=32, the cut off will be around 1320 Hz or thereabouts.

Compensation of the delay. .... What is wrong with my interpretation of this formula?

Your interpretation is correct. If the window length is even, you can't correct for the group delay by simply shifting samples. The code only "approximately" compensates for group delay and is off by half a sample.

There are a few more things wrong about the code:

  1. If someone calls it with half_average > fft_size/4 , bad things will happen
  2. Multiplication in the frequency domain is circular convolution (not linear convolution) in the time domain. The input must be properly zero-padded before FFT to avoid time domain aliasing.
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you, this gives me new pointers to look at. I don't understand n your answer, how you calculate the cutoff as you say $ f_{cutoff} = \frac{f_{s}}{N} $, but you don't state what N, you have used with 44.1 kHz. $\endgroup$ Commented Aug 10, 2020 at 12:40
  • 1
    $\begingroup$ Sorry, fixed. I was using N=32 as an example. $\endgroup$ Commented Aug 10, 2020 at 14:06

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.