Simulating the impulse response of a multirate digital filter

Question

I have designed a CIC filter in Verilog and I would like to test its impulse response, but I am not sure where to start from.

My input sampling frequency is 51200Hz, the filter is a third order decimating CIC with downsampling ratio = 256, so the output sampling frequency is 200Hz. The differential delay was set to one.

The input word is one bit wide, while the output will be 16 bits. I have already calculated all the internal registers dimensions and I am satisfied with that, my current Verilog design does not use true registers but integer numbers so that is not a problem anyway.

How am I supposed to compute its impulse response?

First of all there's the problem that the sequence '1000...' is not an impulse, because 1 corresponds to +1 while 0 corresponds to -1. I have solved that doing two simulations, one with '10101010...' and the other with '11010101010...' as input sequences. The impulse response is then half of the sum of the two responses. So far so good.

Since I want the frequency magnitude response referred with the output sampling frequency I figured that my input impulse should be 256 high frequency samples wide, is that correct? I have done such a simulation and computed the specter in Matlab, the result is qualitatively satisfying but I have a lot of noise, especially in the low frequency part of the specter. I have also tried some windowing without much luck. I simulated my sistem for about half a second, that's 25000 samples.

I must note that my simulation output gives out samples at an higher non constant rate, I've tried resampling my signal with no luck.

Given my limited experience with digital signal processing I am out of ideas. Is there something fundamentally wrong in my approach? Is it ok but there's something I should poke to get better results?

Answer 1

By definition, the Impulse Response of a Filter is the response the filter has to an Impulse: https://en.wikipedia.org/wiki/Impulse_response

Since your Cascaded Integrator Comb Filter is an optimized FIR filter, it's response will be Finite in length (as opposed to IIR - Infinite Impulse Response) and you can calculate this response by sending a single maximum value followed by zeros through it.

You mentioned [1, 0, 0, 0] is not an impulse response since 0 represents -1. I believe the value that represents 0 would be half of the unsigned maximum value in this situation. If your values are Unsigned integers from 0 to 255 with 0 representing the largest negative value (-1) and 255 representing the highest positive value (+1), you should put a maximum value (255) followed by zero values (128). After the filter, you will probably shift/scale/truncate back to the previous bit-width to account for bit growth/gain of the filter. This will give you the values of your Filter's Impulse Response.

You have a right idea for calculating the Impulse Response when you tried input '10101010...' and '11010101010...' (these are really like putting in [1 -1 1 -1 ...] and [1 1 -1 1 -1 1 ...]). Since FIR filters (CIC in this case) are linear systems, you can use the Superposition rule and add the inputs and then filter. You will get the same answer as filtering the signals separately and then adding. That is why your experiment gives you twice the filter impulse response. If you add after you filter, you will still get twice the filter impulse response.

I proved this to myself with this example (MATLAB):

If a is your filter:

a = [1 -1 1 1 1 -1]

and b and c are your inputs:

b = [1 -1 1 -1 1 -1 1 -1 1 -1 1]
c = [1 1 -1 1 -1 1 -1 1 -1 1 -1]

filter(a,1,b) + filter(a,1,c) = [2 -2 2 2 2 -2 0 0 0 0 0]

which is two times your filter impulse response.

However, you still get the filter impulse response by inputting a single max input:

filter(a,1,[1 0 0 0 0 0 0 0 0 0 0 ]) = [1 -1 1 1 1 -1 0 0 0 0 0]

For a very good treatment of CIC filter analysis, see: http://dspguru.com/sites/dspguru/files/cic.pdf

To calculate the Frequency Magnitude Response of your filter, you can take the FFT of the coefficients. To plot it, you should do something like this in MATLAB:

impResp = [1 -1 1 1 1 -1];

freqDomain = fft(impResp);
fSamp = 100;
frequencies = -fSamp/2:fSamp/(numel(impResp)-1):fSamp/2;

figure
plot(frequencies,db(abs(fftshift(freqDomain))))
xlabel('frequencies')
ylabel('db')

This filter does not have a very interesting Impulse response, though.

Answer 2

You should probably start with floating point filter representation in MATLAB. You can specify it as CIC coefficients vector, $cic = [C0, C1, C2]$ for third order one. Next you need a filtering function, something like $y=cicFilter(cic, x)$, you can find the ready one or write it (if you've already done some HDL code for this filter it won't be difficult to implement proper MATLAB code). So you can specify test vector $x = [0,0,0,...,1,0,0,0...]^T$ and use to find impulse response of your filter (there are coefficients $cic$ to be tested, nothing more). By doing $FFT$ of result $y$ you'll obtain the frequency response evidently. For it's magnitude part it will be a smooth curve, the filter magnitude spectrum with no noise, since no noise is presented at the filter input.

After such simple verification you can start to test your HDL implementation. I'm curious about one bit representation of the input signal since filter is designed to work with signed numbers. So maybe you should represent your input as 2 bits signed number instead of one bit unsigned. Then it's possible to feed your filter with test vector the same as we've used in MATLAB ideal model. This test vector are signed integers truncated to 2 bits signed representation. Then capture the output data, store it and analyze in MATLAB. Its frequency magnitude response must be similar to the ideal one. Plot it in Log scale to see fixed point degradation with respect to floating point ideal representation. Also by varying bit width of the input, output and processing operations at each stage you can look different performance with such a plot.

Hope it helps.