Discussion: Interpolation Can of Worms

"Erik de Castro Lopo" <nospam@mega-nerd.com> wrote in message 
news:4507CA5A.663BFAEB@mega-nerd.com...
> dbell wrote:
>>
>> To promote some interesting discussion I thought I would raise the
>> following question:
>>
>> Possibly as part of a sequence of processing, to interpolate a sampled
>> sequence by 2 is it correct to keep the original samples and derive the
>> in-between samples, or replace the original samples in the process of
>> interpolating the signal.
>
>
> I would like to contradict what everyone else in this thread
> has said ;-). I prefer the second method or rather I consider the
> first method wrong and here's why.
>
<snip Erik's case>

One application where it is very important to keep the original "samples"
is when the original is a set of calibration values, and the interpolation 
is
used to provide intermediate cal values (for example, measured amplitude
versus frequency tables, to provide amplitudes at frequencies between the
measurement points).

Regards
    Ian 

Ian wrote:
> 
> One application where it is very important to keep the original "samples"
> is when the original is a set of calibration values, and the interpolation
> is used to provide intermediate cal values (for example, measured amplitude
> versus frequency tables, to provide amplitudes at frequencies between the
> measurement points).

Agree completely. My comments were specifically aimed at
interpolation of audio signals :-).

Erik
-- 
+-----------------------------------------------------------+
  Erik de Castro Lopo
+-----------------------------------------------------------+
"C++ is an atrocity, the bletcherous scab of the computing world,
responsible for more buffer overflows, more security breaches, 
more blue screens of death, more mysterious failures than any 
other computer language in the history of the planet Earth."
-- Eric Lee Green

Erik de Castro Lopo wrote:
> dbell wrote:
[...]
> > Possibly as part of a sequence of processing, to interpolate a sampled
> > sequence by 2 is it correct to keep the original samples and derive the
> > in-between samples, or replace the original samples in the process of
> > interpolating the signal.
>
>
> I would like to contradict what everyone else in this thread
> has said ;-). I prefer the second method or rather I consider the
> first method wrong and here's why.


Hi Erik,
You aren't exactly contradicting what I said.  I conceded it might be
desirable to have a  corner frequency other than the original folding
frequency.  But you are certainly a more eloquent advocate for the
idea.



>
[...]
> Consider a source signal with frequency components all the
> way up to a frequency of
>
>      fx = 0.5 * fs - epsilon
>
> where epsilon is orders of magnitude smaller than fs.


If the first approach (interleave half-sample delay) is used with a
filter having a transition band of epsilon radians (or 2*epsilon Hz),
then no aliasing occurs.  There is no spectral energy there to alias.


For a given order anti-alias filter (compared at 2*input fs).  There is
a tradeoff to be decided between approach 1 and approach 2.

1) keeping all the spectral content within epsilon Hz of the original
folding frequency at the cost of allowing some highly attenuated
aliasing in a band that was probably vacant anyway.

2) eliminating all aliasing at the cost of increasing the transition
band and possibly losing information.  Note, the transition band need
not be doubled exactly, but that's another can of worms.


To be honest, there is no right approach 100% of the time.

The half-sample-delay approach will always be faster, but it will also
always have aliasing in the transition band.  The engineer's job is to
decide when this aliasing is ineffectual and allowable.

-- 
Mark Borgerding
3dB Labs, Inc.
Innovate. Develop. Deliver.

I have been giving some thought to this topic since I posted the
original question for discussion.

A second question came to mind:

What are you really doing when you perform the following process?

1) zero pad to double the sample rate (from Fs to 2*Fs), creating an
undesired spectral image
2) run a good lowpass filter to get rid of the spectral image created
3) replace the data at times corresponding to the original data with
the original data (making the appropriate shift for filter delay)

Is there some simple way to describe the effects of what you are doing
that may help you decide how to answer the original question?

I think there is a simple way.

Lets say you have designed a good N-tap anti-imaging filter {h(n),
n=0..(N-1)) with the following characteristics:

A) lowpass filter
B) gain of 2 for unity gain interpolate-by-2 system
C) odd length, N
D) even symmetric coefficients
E) much greater stop band ripple than pass band ripple
F) pass-band frequency (Fpb) < stop-band frequency (Fsb) <= (0.5*Fs)

If you apply this filter in steps 1-3 above, then what you are
effectively doing is:

1) zero padding the data
2) filtering with a half-band filter that is effectively derived from
h(n)

THAT'S ALL!   NOTHING MORE!

How is the half-band filter specified?

a) starting at the center coefficient, zero every second coefficient
that precedes it
b) starting at the center coefficient, zero every second coefficient
that follows it
c) replace the center coefficient with 1.0

If N is of the form 3+K*4 (K an integer >=0) then the first and last
point will normally be non-zero.
If N is not of the form 3+K*4 then the first and last point will be
zero (so you could shorten the filter length by 2).

What does the half-band filter frequency response look like?

1) not quite the same shape as if you designed a half-band filter
directly
2) assuming E) above, and that the original filter stop band is at
least as wide as the original filter pass band, the ripple in the
half-band pass band will be about half the ripple in the original
filter pass band
3) the half-band stop-band ripple will be the same as the half-band
pass-band ripple
4) the transition bandwidth will be approximately (Fs - 2*Fpb)
5) all the rest of the normal half-band filter characteristics

What does this imply?

1) as stated before, you are effectively re-sampling using a half-band
filter with all of the advantages and disadvantages that go along with
it
2) any great initial performance of the original filter has probably
been wasted
3) you can probably get away with an original filter having less
stop-band attenuation and less taps
4) you should probably design the half-band filter directly instead of
designing the non-half-band filter, so you can control the design and
specs

FURTHER DISCUSSION ON THIS AND THE ORIGINAL QUESTION ENCOURAGED

Dirk Bell
DSP Consultant


dbell wrote:
> To promote some interesting discussion I thought I would raise the
> following question:
>
> Possibly as part of a sequence of processing, to interpolate a sampled
> sequence by 2 is it correct to keep the original samples and derive the
> in-between samples, or replace the original samples in the process of
> interpolating the signal.
> 
> Dirk Bell
> DSP Consultant

dbell wrote:
> 
> I have been giving some thought to this topic since I posted the
> original question for discussion.
> 
> A second question came to mind:
> 
> What are you really doing when you perform the following process?
> 
> 1) zero pad to double the sample rate (from Fs to 2*Fs), creating an
> undesired spectral image
> 2) run a good lowpass filter to get rid of the spectral image created
> 3) replace the data at times corresponding to the original data with
> the original data (making the appropriate shift for filter delay)
> 
> Is there some simple way to describe the effects of what you are doing
> that may help you decide how to answer the original question?
> 
> I think there is a simple way.
> 
> Lets say you have designed a good N-tap anti-imaging filter {h(n),
> n=0..(N-1)) with the following characteristics:
> 
> A) lowpass filter
> B) gain of 2 for unity gain interpolate-by-2 system
> C) odd length, N
> D) even symmetric coefficients
> E) much greater stop band ripple than pass band ripple
> F) pass-band frequency (Fpb) < stop-band frequency (Fsb) <= (0.5*Fs)
> 
> If you apply this filter in steps 1-3 above, then what you are
> effectively doing is:
> 
> 1) zero padding the data
> 2) filtering with a half-band filter that is effectively derived from
> h(n)
> 
> THAT'S ALL!   NOTHING MORE!

Yes, agree completely.

> How is the half-band filter specified?
> 
> a) starting at the center coefficient, zero every second coefficient
> that precedes it
> b) starting at the center coefficient, zero every second coefficient
> that follows it
> c) replace the center coefficient with 1.0
> 
> If N is of the form 3+K*4 (K an integer >=0) then the first and last
> point will normally be non-zero.
> If N is not of the form 3+K*4 then the first and last point will be
> zero (so you could shorten the filter length by 2).
> 
> What does the half-band filter frequency response look like?
> 
> 1) not quite the same shape as if you designed a half-band filter
> directly
> 2) assuming E) above, and that the original filter stop band is at
> least as wide as the original filter pass band, the ripple in the
> half-band pass band will be about half the ripple in the original
> filter pass band
> 3) the half-band stop-band ripple will be the same as the half-band
> pass-band ripple
> 4) the transition bandwidth will be approximately (Fs - 2*Fpb)
> 5) all the rest of the normal half-band filter characteristics
> 
> What does this imply?
> 
> 1) as stated before, you are effectively re-sampling using a half-band
> filter with all of the advantages and disadvantages that go along with
> it
> 2) any great initial performance of the original filter has probably
> been wasted

If you have a true half band filter of the form:

    h[n] = {..., g[-2], 0, g[-1], 0, g[0], 1, g[1], 0, g[2], 0, g[3], ...}

then an efficient polyphase implementation is possible.

Erik
-- 
+-----------------------------------------------------------+
  Erik de Castro Lopo
+-----------------------------------------------------------+
"You can have quality software, or you can have pointer arithmetic;
but you cannot have both at the same time." -- Bertrand Meyer, 1989

dbell wrote:
> To promote some interesting discussion I thought I would raise the
> following question:
>
> Possibly as part of a sequence of processing, to interpolate a sampled
> sequence by 2 is it correct to keep the original samples and derive the
> in-between samples, or replace the original samples in the process of
> interpolating the signal.
>
> Dirk Bell
> DSP Consultant

Hello Dirk,

The 1st part of your question is really answered by whether or not you
want an interpolation (literally in between the samples) or an
approximation.

Since in sampled systems, we already know the proper interpolation
function we can work with approximations of it. The trade off between
interpolation and approximation is an approximation can yield a smaller
maximium error. Although there are many reasons to go with
interpolation, often one hates to see their original values get changed
as an approximation would do. And thus an interpolation "feels" better
in many circumstances. Also Rune's point should be taken that one
should start with a simple method and build upon it. Because if your
final method has a problem, you can back up to a simpler version.
Currently I'm working with using interpolation methods to solve
differential equations.

Clay

in article 1158072374.880336.237110@e3g2000cwe.googlegroups.com, dbell at
bellda2005@cox.net wrote on 09/12/2006 10:46:

> To promote some interesting discussion I thought I would raise the
> following question:
> 
> Possibly as part of a sequence of processing, to interpolate a sampled
> sequence by 2 is it correct to keep the original samples and derive the
> in-between samples, or replace the original samples in the process of
> interpolating the signal.

it could go either way.  depends on the specs (how long the FIR is and how
sharp the cut-off at pi/2 has to be).  the difference is, as you mentioned
Dirk, that of half-band filter vs. not.  it can equivalently be thought of
as windowed sinc() function:

    h[n] = w[n] * sinc(n/2)

probably the best guess for w[n] would be some Kaiser.  pick some cutoff
spec, try it for Kaiser windowed sinc(n/2), then use the same spec to
design, using Parks-McClellan, a filter with -6 dB cutoff at or around
Nyquist/2.  of the two, which one has fewer non-zero values of h[n]?  i've
looked into this in the past and have decided that, for most audio
application, the latter (P-McC) is computationally less costly.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

robert bristow-johnson wrote:

> in article 1158072374.880336.237110@e3g2000cwe.googlegroups.com, dbell at
> bellda2005@cox.net wrote on 09/12/2006 10:46:
>
> > To promote some interesting discussion I thought I would raise the
> > following question:
> >
> > Possibly as part of a sequence of processing, to interpolate a sampled
> > sequence by 2 is it correct to keep the original samples and derive the
> > in-between samples, or replace the original samples in the process of
> > interpolating the signal.
>
> it could go either way.  depends on the specs (how long the FIR is and how
> sharp the cut-off at pi/2 has to be).  the difference is, as you mentioned
> Dirk, that of half-band filter vs. not.  it can equivalently be thought of
> as windowed sinc() function:
>
>     h[n] = w[n] * sinc(n/2)
>
> probably the best guess for w[n] would be some Kaiser.  pick some cutoff
> spec, try it for Kaiser windowed sinc(n/2), then use the same spec to
> design, using Parks-McClellan, a filter with -6 dB cutoff at or around
> Nyquist/2.  of the two, which one has fewer non-zero values of h[n]?  i've
> looked into this in the past and have decided that, for most audio
> application, the latter (P-McC) is computationally less costly.

Some while back, Fred posted a link to a program that computes
equiripple halfband filters
(http://groups.google.ch/group/comp.dsp/msg/5476b541eb7c55ce?dmode=source)
- that should beat the pants off any filter for factor 2 upsampling
(w.r.t. computational efficiency). BTW, we've had this (or similar)
discussion before:

http://groups.google.ch/group/comp.dsp/browse_frm/thread/8de91e30641d10e8/4e69c0ef1c9f4edf?#4e69c0ef1c9f4edf

My 0.02$ is that the appropriate filter (either halfband or Erik's
proper downsampler) depends on the application. For audio sample rate
conversion, halfband filters usually are critical. However, for
sidechain upsampling or metering, halfband filters are very useful.

Regards,
Andor

Andor wrote:

> robert bristow-johnson wrote:
> 
>> in article 1158072374.880336.237110@e3g2000cwe.googlegroups.com, dbell at
>> bellda2005@cox.net wrote on 09/12/2006 10:46:
>>
>> > To promote some interesting discussion I thought I would raise the
>> > following question:
>> >
>> > Possibly as part of a sequence of processing, to interpolate a sampled
>> > sequence by 2 is it correct to keep the original samples and derive the
>> > in-between samples, or replace the original samples in the process of
>> > interpolating the signal.
>>
>> it could go either way.  depends on the specs (how long the FIR is and
>> how
>> sharp the cut-off at pi/2 has to be).  the difference is, as you
>> mentioned
>> Dirk, that of half-band filter vs. not.  it can equivalently be thought
>> of as windowed sinc() function:
>>
>>     h[n] = w[n] * sinc(n/2)
>>
>> probably the best guess for w[n] would be some Kaiser.  pick some cutoff
>> spec, try it for Kaiser windowed sinc(n/2), then use the same spec to
>> design, using Parks-McClellan, a filter with -6 dB cutoff at or around
>> Nyquist/2.  of the two, which one has fewer non-zero values of h[n]? 
>> i've looked into this in the past and have decided that, for most audio
>> application, the latter (P-McC) is computationally less costly.
> 
> Some while back, Fred posted a link to a program that computes
> equiripple halfband filters
> (http://groups.google.ch/group/comp.dsp/msg/5476b541eb7c55ce?dmode=source)
> - that should beat the pants off any filter for factor 2 upsampling
> (w.r.t. computational efficiency). BTW, we've had this (or similar)
> discussion before:
> 

I have thought that for interpolation, equiriple is not necessarily the best
possible solution. Especially in the case of audio. I believe aliasing is
more audible for frequencies close to 2*Nyquist than close to Nyquist
frequency, thus it would be better to use some form of weighting based on
subjective sound quality. Just a thought...

-- 
Jani Huhtanen
Tampere University of Technology, Pori

Erik de Castro Lopo wrote:

> dbell wrote:
> 
>>To promote some interesting discussion I thought I would raise the
>>following question:
>>
>>Possibly as part of a sequence of processing, to interpolate a sampled
>>sequence by 2 is it correct to keep the original samples and derive the
>>in-between samples, or replace the original samples in the process of
>>interpolating the signal.
> 
> 
> 
> I would like to contradict what everyone else in this thread
> has said ;-). I prefer the second method or rather I consider the
> first method wrong and here's why.
> 
> So, for the upsample by 2 case lets use:
> 
>     fs -- source sample rate
> 	fd -- desitnation sample rate (fs * 2)
> 
> Consider a source signal with frequency components all the
> way up to a frequency of
> 
>      fx = 0.5 * fs - epsilon
> 
> where epsilon is orders of magnitude smaller than fs.
> 
> Now upsample the source signal by inserting a zero valued sample
> between each of the original samples and filter the resulting
> signal (at sample rate of fd) using a half band FIR anti aliasing
> filter of the form:
> 
>     h[n] = {..., g[-2], 0, g[-1], 0, g[0], 1, g[1], 0, g[2], 0, g[3], ...}
> 
> Notice however than the frequency response of this filter has
> a passband magnitude response of 2, a -3dB point occuring at
> 0.5 of the source sample rate. This looks something like
> (ASCII art):
> 
>     +---------------------     -3dB point
>     |                     \  /
>     |                      \
>     |                      |\
>     |                      | \
>     |                      |  \
>     |                      |   \
>     |                      |    \
>     |                      |     \
>     |                      |      \
>     |                      |       \
>     +----------------------+----------------------+-
>                      fs/2                        fd/2
> 
> 
> So whats the problem here? Well remember that source frequency at a
> frequency of fx? Well, that frequency gets aliased to 
> 
>     fx1 = 0.5 * fs + epsilon
> 
> and the attenuation of the anti aliasing filter at fx1 is only about
> 3dB. Now normally, fs/2 is 22.05kHz or more so that the original 
> fx and the alias frequency fx1 are both out of hearing range. However,
> they do exist and at least one of them shouldn't be there. These 
> alias frequencies should be obvious when upsample from say 8kHz to
> 16kHz.
> 
> As far as I am concerned the filter should look like this:
> 
>     +---------------  
>     |               \ 
>     |                \
>     |                 \
>     |                  \
>     |                   \
>     |                    \     target stopband attenuation
>     |                     \  /  
>     |                      |
>     |                      |\
>     |                      | \
>     +----------------------+----------------------+-
>                      fs/2                        fd/2
> 
> 
> The anti-aliasing filters in Secret Rabbit Code
> 
>     http://www.mega-nerd.com/SRC/
> 
> are designed like this.
> 
> Erik
> 
> 
> 
> Erik


Yes, but half-band filters are generally used with other filters.  If 
you look at the part that is only 6 db down (it is 6dB at Fs/4 for a 
halfband, not 3dB),  it gets folded onto the top part of the pass band 
for this filter.  i.e. Fs/4+delta folds onto Fs/4-delta.  That aliasing 
region is in the stopband of the next filter in a multi-rate filter 
bank.  You generally don't want to use a half-band decimator (or 
interpolator) as the final filter (first in interpolation) because of 
the aliasing.  The advantage of using the halfband filters of course is 
the greatly reduced computational power needed to implement them.  You 
can typically get away with 15 taps or less in the halfband decimators, 
and 8 of those taps are zeros.