fft - How to create a distortion compensation filter

I am creating a channelizer that consists of a complex mixer, CIC decimator, and a FIR compensation/decimation filter. The final FIR filter may, if it matters, be implemented as more than one filter.

My question is, how do I design a filter such that it compensates for the very non-flat frequency response of the CIC filter? Do you create the frequency response that you want by calculating the reciprocal of the CIC's response and then run it through an inverse FFT to get the impulse response?

As you can see, though my particular problem pertains to CIC filters, my question is really about how you create any kind of distortion compensation filter.

Thanks for your time.

Answer

The following addresses the CIC compensation and is not a general "distortion" technique. But it is a straight-forward method to "compensate" non-varying "distortion". If the frequency response is known the inverse of the frequency can be used to compensate. Examples like the CIC filter, where a poor filter might be used because of reduced complexity, are compensated later the signal chain. In this example the frequency response is know and the inverse can be used. Note, with multi-rate filters you only want to use the "usable" spectrum after decimation.

Generally, to compensate a CIC filter the inverse of the CIC filters response can be used to generate the compensation filter. The CIC has a response of (see reference [r2] or [r3])

$$ H(\omega) = \left| \frac{sin(\omega D/2)}{sin(\omega M/2)} \right|^N $$

Where D is the differentiate delay, M is the decimation rate, and N is the filter order (number of cascaded filters). The inverse can be specified as

$$ H(\omega) = \left| \frac{sin(\omega M/2)}{sin(\omega D/2)} \right|^N $$

Once we have the frequency response of the compensation filter, we can simply choose the length of FIR filter that we desire. The length of the FIR is application specific. Obviously the longer the FIR filter the better compensation.

The following are plots of this straight forward compensation.

The following is the Python code to create the frequency responses and plots.

import numpy as np
from numpy import sin, abs, pi
import pylab

D = 1; M = 7; N = 3

Hfunc = lambda w : abs( (sin((w*M)/2)) / (sin((w*D)/2.)) )**N
HfuncC = lambda w : abs( (sin((w*D)/2.)) / (sin((w*M)/2.)) )**N

w = np.arange(1024) * pi/1024


G = (M*D)**N
H = np.array(map(Hfunc, w))
Hc = np.array(map(HfuncC, w))
# only use the inverse (compensation) roughly to the first null.
Hc[int(1024*pi/M/2):] = 1e-8
plot(w, 20*log10(H/G))
plot(w, 20*log10(Hc*G))
grid('on')

See [r1] for other approaches and $sinc^{-1}$ approximation.

[r1] Altera, "Understanding CIC compensation filters"

[r2] R. Lyons, "Understanding Digital Signal Processing", 2nd ed., Prentice Hall, Upper Saddle River, New Jersey, 2004

[r3] R. Lyons, "Understanding Cascaded Integrator Comb Filters"