Trying to understand CIC interpolation

I am using Matlab's mfilt.cicinterp on a piece of data. Out of curiousity,
I also tried to write my own CIC interpolation filter. I referred to
Hogenhauers original paper (An Economical Class of Digital Filters for
Decimation and Interpolation) as well as Matlab's documentation for
mfilt.cicinterp at
http://www.mathworks.com/access/helpdesk/help/toolbox/filterdesign/index.html?/access/helpdesk/help/toolbox/filterdesign/ref/mfilt.cicinterp.html

Based on that I have the following dilemma:

>> u=(-1).^(1:20); % alternating +1 and -1;
>> Hm=mfilt.cicinterp(2,1,1); % 2x interpolation with one stage CIC 
>> y1=filter(Hm, u); 

This chunk of code gives me y1= [0 -1 -1 1 1 -1 -1 1 1 ... -1 -1 1]; 

Now I try to simulate Matlab's CIC block diagram in following code:

>> function z=CICInt(u,R,N)

>> L=length(u);
>> y=u; % buffer for differentiator stage
>> y1=0*y;
>> dcG=R^(N-1); % dc gain

>> for i=1:N
>>    for j=2:L
>>        y1(j)=y(j)-y(j-1); % differentiator
>>    end
>>    y=y1;
>>    y1=0*y;
>> end
>> z=upsample(y,R); % upsampling
>> z1=0*z; % initialize z1, output of each stage
>> for i=1:N
>>    for j=2:L*R
>>        z1(j)=z(j)+z1(j-1); % delaying integrator
>>    end
>>    z=[0 z1(1:end-1)]; % delay output for next stage and call it input
for next stage
>>    z1=0*z; % re-initialize output buffer
>> end
>> z=z/dcG; % scale output

When I run this on the same input I get [0 0 0 2 2 0 0 2 2 0 0 2 2 0 0 ...
0 2].

Where am I going wrong? I put in a probe point after the upsampling and
the input to the integrator stage is [0 0 2 0 -2 0 2 0 -2 ... 2 0] which
makes sense. 

Probably something immensely trivial that I am missing, but 4 cups of
coffee later, I need help :-)



_____________________________________
Do you know a company who employs DSP engineers?  
Is it already listed at http://dsprelated.com/employers.php ?

On Mon, 21 May 2007 12:40:38 -0500, "somdeb" <somdeb@gmail.com> wrote:

>I am using Matlab's mfilt.cicinterp on a piece of data. Out of curiousity,
>I also tried to write my own CIC interpolation filter. I referred to
>Hogenhauers original paper (An Economical Class of Digital Filters for
>Decimation and Interpolation) as well as Matlab's documentation for
>mfilt.cicinterp at
>http://www.mathworks.com/access/helpdesk/help/toolbox/filterdesign/index.html?/access/helpdesk/help/toolbox/filterdesign/ref/mfilt.cicinterp.html
>

   (snipped)

Hello somdeb,

  I don't know if it will help you 
much, or not, but you might take a look 
at:

http://www.embedded.com/shared/printableArticle.jhtml?articleID=160400592

Good Luck,
[-Rick-]

On May 21, 10:40 pm, "somdeb" <som...@gmail.com> wrote:
> I am using Matlab's mfilt.cicinterp on a piece of data. Out of curiousity,
> I also tried to write my own CIC interpolation filter. I referred to
> Hogenhauers original paper (An Economical Class of Digital Filters for
> Decimation and Interpolation) as well as Matlab's documentation for
> mfilt.cicinterp athttp://www.mathworks.com/access/helpdesk/help/toolbox/filterdesign/in...
>
> Based on that I have the following dilemma:
>
> >> u=(-1).^(1:20); % alternating +1 and -1;
> >> Hm=mfilt.cicinterp(2,1,1); % 2x interpolation with one stage CIC
> >> y1=filter(Hm, u);
>
> This chunk of code gives me y1= [0 -1 -1 1 1 -1 -1 1 1 ... -1 -1 1];
>
> Now I try to simulate Matlab's CIC block diagram in following code:
>
>
>
>
>
> >> function z=CICInt(u,R,N)
> >> L=length(u);
> >> y=u; % buffer for differentiator stage
> >> y1=0*y;
> >> dcG=R^(N-1); % dc gain
> >> for i=1:N
> >>    for j=2:L
> >>        y1(j)=y(j)-y(j-1); % differentiator
> >>    end
> >>    y=y1;
> >>    y1=0*y;
> >> end
> >> z=upsample(y,R); % upsampling
> >> z1=0*z; % initialize z1, output of each stage
> >> for i=1:N
> >>    for j=2:L*R
> >>        z1(j)=z(j)+z1(j-1); % delaying integrator
> >>    end
> >>    z=[0 z1(1:end-1)]; % delay output for next stage and call it input
> for next stage
> >>    z1=0*z; % re-initialize output buffer
> >> end
> >> z=z/dcG; % scale output
>
> When I run this on the same input I get [0 0 0 2 2 0 0 2 2 0 0 2 2 0 0 ...
> 0 2].

Looking at the code snippet which you have provided, it seems a bit
too verbose than
what is required. The output seems to be wrong too... :-)
If I understand correctly, your objective is to
u=(-1).^(1:20);
u1 = [u ; zeros(1,length(u))];
u1 = u1(:).'; % zero inserted
o = conv(u1, ones(1,2)); % op of filtering

for a simple for-loop code (matching with the hardware block diagram),
maybe you can have a look @ the
blog post http://dsplog.blogspot.com/2007/07/example-of-cascaded-integrator-comb.html

>
> Where am I going wrong? I put in a probe point after the upsampling and
> the input to the integrator stage is [0 0 2 0 -2 0 2 0 -2 ... 2 0] which
> makes sense.
>
> Probably something immensely trivial that I am missing, but 4 cups of
> coffee later, I need help :-)

Hope this helps.

Regards,
Krishna
http://dsplog.blogspot.com

>
> _____________________________________
> Do you know a company who employs DSP engineers?
> Is it already listed athttp://dsprelated.com/employers.php?

On May 21, 1:40 pm, "somdeb" <som...@gmail.com> wrote:
> I am using Matlab's mfilt.cicinterp on a piece of data. Out of curiousity,
> I also tried to write my own CIC interpolation filter. I referred to
> Hogenhauers original paper (An Economical Class of Digital Filters for
> Decimation and Interpolation) as well as Matlab's documentation for
> mfilt.cicinterp athttp://www.mathworks.com/access/helpdesk/help/toolbox/filterdesign/in...
>
> Based on that I have the following dilemma:
>
> >> u=(-1).^(1:20); % alternating +1 and -1;
> >> Hm=mfilt.cicinterp(2,1,1); % 2x interpolation with one stage CIC
> >> y1=filter(Hm, u);
>
> This chunk of code gives me y1= [0 -1 -1 1 1 -1 -1 1 1 ... -1 -1 1];
>
> Now I try to simulate Matlab's CIC block diagram in following code:
>
>
>
>
>
>
>
> >> function z=CICInt(u,R,N)
> >> L=length(u);
> >> y=u; % buffer for differentiator stage
> >> y1=0*y;
> >> dcG=R^(N-1); % dc gain
> >> for i=1:N
> >>    for j=2:L
> >>        y1(j)=y(j)-y(j-1); % differentiator
> >>    end
> >>    y=y1;
> >>    y1=0*y;
> >> end
> >> z=upsample(y,R); % upsampling
> >> z1=0*z; % initialize z1, output of each stage
> >> for i=1:N
> >>    for j=2:L*R
> >>        z1(j)=z(j)+z1(j-1); % delaying integrator
> >>    end
> >>    z=[0 z1(1:end-1)]; % delay output for next stage and call it input
> for next stage
> >>    z1=0*z; % re-initialize output buffer
> >> end
> >> z=z/dcG; % scale output
>
> When I run this on the same input I get [0 0 0 2 2 0 0 2 2 0 0 2 2 0 0 ...
> 0 2].
>
> Where am I going wrong? I put in a probe point after the upsampling and
> the input to the integrator stage is [0 0 2 0 -2 0 2 0 -2 ... 2 0] which
> makes sense.
>
> Probably something immensely trivial that I am missing, but 4 cups of
> coffee later, I need help :-)
>
> _____________________________________
> Do you know a company who employs DSP engineers?  
> Is it already listed athttp://dsprelated.com/employers.php?- Hide quoted text -
>
> - Show quoted text -

somedeb,

This is certainly a late response to the original post, but WTH,
Krishna popped it up to the top of the list.

The difference in your two outputs appears to be because:

1) The first output of your differentiator is not set properly.
2) The first output of your integrator is not set properly.

Comments:

1) The delay of the output you intoduce by inserting a leading zero is
only required to align your output to that from your 3 line MATLAB
code and can be eliminated in practice.
2) The CIC interpolation filter is really intended to be used on
something that is already oversampled, your test signal is not.
Consequently, the filter gain giving the desired output signal is not
the low freq gain of the filter set in the MATLAB script, because it
is not a low freq.
3) Test your final code rigorously on some other signals.


Dirk

[snip]
> 2) The CIC interpolation filter is really intended to be used on
> something that is already oversampled, your test signal is not.
> Consequently, the filter gain giving the desired output signal is not
> the low freq gain of the filter set in the MATLAB script, because it
> is not a low freq.

In somdeb's code, there is code for upsampling (in between
differentiator code and integrator code).

Regards,
Krishna
http://dsplog.blogspot.com

On Jul 4, 8:17 am, Krishna <krishna.pil...@gmail.com> wrote:
> [snip]
>
> > 2) The CIC interpolation filter is really intended to be used on
> > something that is already oversampled, your test signal is not.
> > Consequently, the filter gain giving the desired output signal is not
> > the low freq gain of the filter set in the MATLAB script, because it
> > is not a low freq.
>
> In somdeb's code, there is code for upsampling (in between
> differentiator code and integrator code).
>
> Regards,
> Krishnahttp://dsplog.blogspot.com


Krishna,

All that code does is stuff zeros between original samples. That is
part of the CIC filter. A CIC interpolator is intended to be applied
to oversampled original signals.

Dirk

On Jul 4, 1:53 pm, dbell <bellda2...@cox.net> wrote:
> On Jul 4, 8:17 am, Krishna <krishna.pil...@gmail.com> wrote:
>
> > [snip]
>
> > > 2) The CIC interpolation filter is really intended to be used on
> > > something that is already oversampled, your test signal is not.
> > > Consequently, the filter gain giving the desired output signal is not
> > > the low freq gain of the filter set in the MATLAB script, because it
> > > is not a low freq.
>
> > In somdeb's code, there is code for upsampling (in between
> > differentiator code and integrator code).
>
> > Regards,
> > Krishnahttp://dsplog.blogspot.com
>
> Krishna,
>
> All that code does is stuff zeros between original samples. That is
> part of the CIC filter. A CIC interpolator is intended to be applied
> to oversampled original signals.
>
> Dirk

Krishna,

Rereading your response again I am now not sure what part of the
comment 2) you were responding to.  The zero padding is taken into
account when the gain is calculated at the end of the script.  If you
interpolate alternating signed ones (magnitude 1 sine wave sampled at
peaks, ... 1,-1,1,-1,1,...) and get pairs of alternating signed ones
(...1,1,-1,-1,1,1,-1,-1,,...) that does not correspond to samples from
an amplitude 1 sign wave.  A CIC interpolator has gain near 1 near 0
Hz, and good attenuation of interpolation-related images for a GROSSLY
oversampled (ie. severely lowpass frequency limited ) original signal.

If I haven't responded to your objection, please restate.

Regards,

Dirk

Dirk,

I was trying to comment on the statement:
2) The CIC interpolation filter is really intended to be used on
> > something that is already oversampled, your test signal is not.
I interpreted this as - "The input signal is not oversampled - hence
the code is wrong". Maybe I made the wrong interpretation :)
In the previous post, I tried to point out that the code has an
upsampling stage inbetween differentiator and integrator stage
(Refer the architecture shown in Figure9 of
http://www.us.design-reuse.com/articles/article10028.html).
Hence your statement on 'upsampling not being done' was not accurate.

Additionally, I do not understand the comment:
> If you
> interpolate alternating signed ones (magnitude 1 sine wave sampled attem
> peaks, ... 1,-1,1,-1,1,...) and get pairs of alternating signed ones
> (...1,1,-1,-1,1,1,-1,-1,,...) that does not correspond to samples from
> an amplitude 1 sign wave.

Repeating the same sample twice is a form of filtering with almost
sinc() shaped transfer
function. Agree?

Regards,
Krishna
http://dsplog.blogspot.com

On Jul 4, 10:40 pm, Krishna <krishna.pil...@gmail.com> wrote:
> Dirk,
>
> I was trying to comment on the statement:
> 2) The CIC interpolation filter is really intended to be used on> > something that is already oversampled, your test signal is not.
>
> I interpreted this as - "The input signal is not oversampled - hence
> the code is wrong".

Krishna,

I have no idea why you would make that interpretation. The comment had
nothing to do with the code, just the test signal.

Remember the CIC filter is really a pretty poor filter as a general
lowpass filter.  When interpolating, it greatly attenuates images over
a very small region/regions, so the original lowpass bandwidth whose
images (resulting from the CIC filter) you want filtered out must be
very narrow bandwidth to start with.  Normally the input signal would
be extremely lowpass bandlimited (compared to its sample rate), either
because it was originally bandlimited, or it was interpolated up to
be, by interpolation methods more appropriate for a signal that has
not been severely bandlimited relative to its original sample rate.  A
CIC interpolation filter is often used as the final filter(s) in a
chain of interpolation filters that might be [FIR filter,
interpolating halfband filters, CIC filters], such as part of a chip
that does 'Digital Up Conversion(DUC)'.

The test signal is not in a range that you would normally use a CIC
filter on for interpolation.  You would probably only use it because
you know what the output should be for that particular input, as one
of your test signals for the code. A further shortcoming of the test
signal is explained below.

Maybe I made the wrong interpretation :)
> In the previous post, I tried to point out that the code has an
> upsampling stage inbetween differentiator and integrator stage
> (Refer the architecture shown in Figure9 ofhttp://www.us.design-reuse.com/articles/article10028.html).
> Hence your statement on 'upsampling not being done' was not accurate.
>

I never made the statement that "upsampling not being done".  What
post are you refering to?

> Additionally, I do not understand the comment:
>
> > If you
> > interpolate alternating signed ones (magnitude 1 sine wave sampled attem
> > peaks, ... 1,-1,1,-1,1,...) and get pairs of alternating signed ones
> > (...1,1,-1,-1,1,1,-1,-1,,...) that does not correspond to samples from
> > an amplitude 1 sign wave.
>

First, excuse my proofreading, that should be "sine wave" not "sign
wave". Also, assuming the test signal is intending to represent a
sinusoid, my implication that it had to come from an amplitude 1
sinusoid is incorrect.  In fact, the sinusoidal signal amplitude of
the signal that the test signal could have come from cannot be
determined from the samples, because the signal frequency is half the
sample rate. With the appropriate phase it could be samples of a
sinusoid with signal frequency of half the input sample rate and ANY
amplitude >= 1.0.  So you cannot use the test signal to determine a
non-zero system gain at the test signal frequency. On the other hand,
the output from the test signal (assuming a sinusoid) does have an
unambiguous amplitude of sqrt(2).

> Repeating the same sample twice is a form of filtering with almost
> sinc() shaped transfer
> function. Agree?

Repeating the same signal twice is adequate for interpolation only if
your input is already sufficiently oversampled. If you are looking for
a flat requency response, a sinc() shaped transfer function is pretty
poor.

Regards,

Dirk

>
> Regards,
> Krishnahttp://dsplog.blogspot.com

Thanks Dirk. I understand your perspective now. Your comments were on
the appropriateness of using
CIC filter for the example test signal and I interpreted them as
comments on the code - yes, it finally
seeped in :)

Thanks,
Krishna
http://dsplog.blogspot.com