2
$\begingroup$

I'm implementing (in software) FIR filters. I have done two implementation: convolutional one and based on multiplication in frequency domain with FFT/IFFT. I'm using standard FFT implementation (not my own) and overlap-add processing.

FFT parameters are automatically selected to be like this: $fftSize = 2 * 2^{\lceil\log2(nTaps)\rceil}$, $samplesPerOp = fftSize - nTaps + 1$ and $overlapSize = nTaps - 1$.

Now I'm testing my implementation. I've generated signal as $\cos(F_1*2*pi*n/F_s) + \cos(F_2*2*pi*n/F_s)$, where $F_1 = 10$, $F_2 = 1000$ and $F_s = 44100$. Also, I design two set of taps via windowed-sinc method: one for low-pass, with transition band $110Hz - 900Hz$ and other high-pass with same transition band. I've tried different windows - Hamming, Blackman-Harris and others.

I apply these two filters to signal, with convoluted and FFT implementations. Low-pass filter works perfectly well in both cases: when I select $nTaps$ large enough (via estimation formula from textbook - "Multirate signal processing" by frederic j. harris, for 90dB attenuation it gives me $nTaps = 327$ taps) I got clean 10Hz signal from both implementations (they are differ on edges of signal, but it is expected, of course)

But high-pass filter behaves strange. Convolutional implementation works. It gives almost undesturbed 1000Hz sine wave. And FFT implementation returns a mess. Result of FFT filtering looks like 1000Hz signal, but it is very corrupted - phase jitter, amplitude jitter, etc.

It looks like, I do something wrong, but I could not understand what. Each textbook, tutorial and article about using FFT for FIR filtering doesn't mention something like this.

Oh, I'm using double, not fixed-point arithmetic, so precision loss should not be a problem.

And, not, it is not HW :)

Here are pictures from gnuplot (1/5 of «second»):

Low-Pass, convolution on top, FFT on down: http://lev.serebryakov.spb.ru/dsp/dsp-low-pass-side-by-side.png
High-Pass, convolution on top, FFT on down: http://lev.serebryakov.spb.ru/dsp/dsp-high-pass-side-by-side.png

$\endgroup$
  • $\begingroup$ If you are trying to filter with the FFT via Overlap-Add, or Overlap-Save, then you shouldn't use any window at all. $\endgroup$ – David Oct 10 '14 at 15:30
  • $\begingroup$ It is window for filter coefficients, not for data. $\endgroup$ – blacklion Oct 10 '14 at 16:08
3
$\begingroup$

Ok, I found bug in my code, and it could be useful to post answer here, as I use «standard» FFT library, and somebody could repeat my mistake in future.

I'm using JTransforms library for Java to make FFT. It uses very compact representation of results for FFT of pure-real signal: it returns only half of components and, what is really important here, first two elements of result array are not (re, im) pair but two real numbers for $X_0$ and $X_{N/2}$ (which are real in case of pure-real data).

My mistake was, that I multiply two vectors (fft(taps) and fft(signal)) as vectors of complex numbers and thrash these two first components! In low-pass version it doesn't matter much, but for high-pass it is critical.

So, proper procedure if you use results of JTransform, is to multiple first two elements as real numbers, and only multiple next $N-2$ elements as $N/2-1$ complex numbers.

It works!

$\endgroup$
1
$\begingroup$

Looks like your approach is correct in principle so it's probably a bug in the actual implementation.

Your corrupted output is clearly "blocky" so it looks like a framing error of sorts. Your corrupted high pass output has clearly blocks of constant amplitude that look like the all the same length. This length is most likely related to either the filter length or the FFT size. It's also worth checking whether you have a phase discontinuity at the block boundary as well. The more info you can provide, the easier it is to speculate about the potential problem.

Hard to tell without actually seeing the code.

$\endgroup$

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.