FIR Interpolation question

"Mark Borgerding" <mark@borgerding.net> wrote in message
news:gk_Tb.78797$DE.72575@fe2.columbus.rr.com...
>
>
> Rick Lyons wrote:
> > On 3 Feb 2004 11:45:20 -0800, mark@borgerding.net (Mark Borgerding)
>
> [snip]
>
> >>This technique works great if you can fit your entire sequence to be
> >>upsampled into a single fft.  For sequences where this is not true,
> >>you could get noise at the buffer boundaries.
>
> [snip]
>
> > Hi Mark,
> >
> >   you idea works for periodic signals (signals
> > whose spetra conatin *NO* leakage).
> > For real-world signals the "frequency-domain zero
> > stuffing to implement time-domain interpolation"
> > method induces errors in the interpolated
> > time-domain sequence.
> >
> > [-Rick-]
> >
>
> Yes, I agree with your observation.
>
> Do you disagree with the above disclaimer?
>
> Let me rephrase: It is a technique that works for 1 block, but not for a
> stream.
>
> The method can be applied to longer streams, BUT one would have to use
> overlapping windows to remove the edge effects.  This added complexity
> makes it less attractive.
>
> If you have a short sequence that you want to up(re)sample,
> frequency-zero-stuffing is a quick way to make sure that you get the
> same spectral components as in the original sampled sequence. *Up to the
> resolution of the original sequence*.
>
> Spectral analysis of any finite length discrete-time signal is limited
> to the frequency resolution dictated by its sampling rate and time
> duration.  The zero-stuffing method achieves accuracy up to that level,
> and no more.
>
> Agreed?

Not exactly....
The sampling rate has nothing to do with spectral resolution - just spectral
extent.
Only the time duration affects the spectral resolution.
Extending the length of the number of samples in frequency by adding a block
of zeros is the dual of extneding the legnth of the number of samples in
time by adding a block of zeros.
The former gives the appearance of adding temporal resolution /
interpolates.
The latter gives the appearance of adding spectral resolution.

Rick touched on the problem with adding zeros in frequency for the purpose
of getting temporal interpolation.  What he said about periodic functions
with no leakage is certainly a sufficient condition - because it's the most
stringent requirement imaginable.
Whether it's a necessary condition can be figured out from this I believe:

The addition of a block of zeros in frequency is equivalent to doing two
operations:
1) Repeat the spectrum for as many times as you need.  This is the same as
(the dual of) inserting zeros into a sequence of temporal samples to
increase the sample rate by 2:1, 3:1, 4:1, etc.
2) Multiply the spectrum by a rectangular window / perfect lowpass filter -
which makes all the samples around fs/2 go to zero.  This is the same as
circularly convolving in time with a sinc or Dirichlet which can cause
temporal aliasing.  The question is a matter of degree.

If the multiplication by the rectangular window causes a sharp discontinuity
in the spectrum at the edges, then the sinc in time, shall we say, has a
pronounced affect.
If the multiplication by the rectangular window doesn't cause a
discontinuity in the spectrum - and perhaps we could have used a less
drastic window to yield "nearly the same" results - then the temporal
aliasing may be acceptable.

Case in point:  Instead of multiplying by a perfect rectangular window in
frequency, we use the response of a lowpass FIR filter.  If the filter is
good enough then the stopband ripple is small and it's as if we stuffed the
spectrum with zeros.  And, the corresponding convolution in time is with a
FIR of known length - so temporal aliasing can be completely avoided by
proper sizing of the arrays as would always be done with circular
convolution.  Note that such a filter can't create sharp discontinuities at
the band edge as is done with a rectangular window.

So, this suggests a compromise approach:
1) Replicate the original spectrum by a factor of 2.
2) Filter with a maximally flat halfband lowpass filter so that there is at
least a double zero at fs/2.
3) Further extend the spectrum by adding a block of zeros around fs/2.
(This is the same as multiplying a replicated spectrum by a window that ends
at fs/4).
Because of the filter in step 2, there's no sharp discontinuity introduced
at fs/4.

I haven't completely characterized this approach.... that's a work in
progress.  There will be some temporal aliasing but not as bad as if one
lets a sharp discontinuity occur.

Now, most FFTs use powers of 2 samples or at least an even number so there
is always a sample exactly at fs/2.  If the sampling theorem is met, then
this sample at fs/2 must be zero.  So, a first order discontinuity at fs/2
shoud not be possible, only a discontinuity in the first derivative.  As a
practical matter this may not be the case.  So, having a maximally flat stop
band filter with multiple zeros at (the new) fs/4 should reduce the temporal
aliasing had by simply adding a block of zeros without first filtering.

Fred

Fred Marshall wrote:
> "Mark Borgerding" <mark@borgerding.net> wrote in message
> news:gk_Tb.78797$DE.72575@fe2.columbus.rr.com...
> 
>>
>>Rick Lyons wrote:
>>
>>>On 3 Feb 2004 11:45:20 -0800, mark@borgerding.net (Mark Borgerding)
>>
>>[snip]
>>
>>
>>>>This technique works great if you can fit your entire sequence to be
>>>>upsampled into a single fft.  For sequences where this is not true,
>>>>you could get noise at the buffer boundaries.
>>
>>[snip]
>>
>>
>>>Hi Mark,
>>>
>>>  you idea works for periodic signals (signals
>>>whose spetra conatin *NO* leakage).
>>>For real-world signals the "frequency-domain zero
>>>stuffing to implement time-domain interpolation"
>>>method induces errors in the interpolated
>>>time-domain sequence.
>>>
>>>[-Rick-]
>>>
>>
>>Yes, I agree with your observation.
>>
>>Do you disagree with the above disclaimer?
>>
>>Let me rephrase: It is a technique that works for 1 block, but not for a
>>stream.
>>
>>The method can be applied to longer streams, BUT one would have to use
>>overlapping windows to remove the edge effects.  This added complexity
>>makes it less attractive.
>>
>>If you have a short sequence that you want to up(re)sample,
>>frequency-zero-stuffing is a quick way to make sure that you get the
>>same spectral components as in the original sampled sequence. *Up to the
>>resolution of the original sequence*.
>>
>>Spectral analysis of any finite length discrete-time signal is limited
>>to the frequency resolution dictated by its sampling rate and time
>>duration.  The zero-stuffing method achieves accuracy up to that level,
>>and no more.
>>
>>Agreed?
> 
> 
> Not exactly....
> The sampling rate has nothing to do with spectral resolution - just spectral
> extent.
> Only the time duration affects the spectral resolution.

If you specify frequency in radians, yes.  Hertz, no.

Mark Borgerding wrote:

> Fred Marshall wrote:

   ...

>> Not exactly....
>> The sampling rate has nothing to do with spectral resolution - just 
>> spectral
>> extent.
>> Only the time duration affects the spectral resolution.
> 
> 
> If you specify frequency in radians, yes.  Hertz, no.

I assume you mean radians per second, but I'm still confused. Specifying
the radian frequency determines the cycle frequency, and vice versa.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message news:<8-2dna0ZScxYNL3dRVn-sQ@centurytel.net>...
> "Mark Borgerding" <mark@borgerding.net> wrote in message
[snip]
> 2) Multiply the spectrum by a rectangular window / perfect lowpass filter -
> which makes all the samples around fs/2 go to zero.  This is the same as
> circularly convolving in time with a sinc or Dirichlet which can cause
> temporal aliasing.  The question is a matter of degree.

Eureka! Capisco! Comprendo! Verstehe!
I get it! 

Discontinuity between the last & first sample in the seqeunce will
cause a sinc ripple error at the edges.  The zero crossings of the
sinc coincide with the original sequence points.  Those zero crossings
are what make it *seem* to have no error, since all the error occurs
between the original points.


Thanks for the explanation.


- Mark

"Mark Borgerding" <mark@borgerding.net> wrote in message
news:5d6e06ef.0402040952.15115b65@posting.google.com...
> "Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:<8-2dna0ZScxYNL3dRVn-sQ@centurytel.net>...
> > "Mark Borgerding" <mark@borgerding.net> wrote in message
> [snip]
> > 2) Multiply the spectrum by a rectangular window / perfect lowpass
filter -
> > which makes all the samples around fs/2 go to zero.  This is the same as
> > circularly convolving in time with a sinc or Dirichlet which can cause
> > temporal aliasing.  The question is a matter of degree.
>
> Eureka! Capisco! Comprendo! Verstehe!
> I get it!
>
> Discontinuity between the last & first sample in the seqeunce will
> cause a sinc ripple error at the edges.  The zero crossings of the
> sinc coincide with the original sequence points.  Those zero crossings
> are what make it *seem* to have no error, since all the error occurs
> between the original points.

I don't think this quite describes what happens.  Let's see....

If we start with a continuous temporal function that is periodic and if we
select exactly an integer number of periods to sample in our "window" or
temporal epoch, then we can say this:

There is no aliasing in frequency due to discrete sampling in time as along
as the sample rate within the epoch is >2B where B is the frequency of the
highest sinusoid in the Fourier Series of the periodic waveform.
There *is* aliasing in frequency if the same waveform is sampled at a rate
>2B as above but a non-integer number of cycles of the waveform are in the
window or temporal epoch.  That is because the window has introduced sharp
discontinuities at the edges.
Note that the sample rate chosen is no longer large enough because the
bandwidth of the continuous windowed function is driven by the edges now and
not by the underlying components of the periodic waveform.  Thus the
frequency aliasing.
In effect, the spectrum is convolved with a sinc and becomes infinite.  This
is the Gibbs phenomenon - but in frequency.

So the "sinc ripple error at the edges" that you're talking about is in
frequency, not time - if I understand what you're saying.

The distance between zeros in this continuous frequency sinc is dependent on
the length of the temporal epoch.  So, I don't see how you can say they
coincide with the original function sample points - those sequence points
are in time and the sinc here is in frequency.

Perhaps you're thinking about the classical reconstruction of strictly
bandlimited functions - of which some periodic functions are a subset.  If
you sample a strictly bandlimited (thus infinite) function of bandwidth B at
a sample rate fs>2B then you can conceptually reconstruct in time with
infinite continous sincs that are related to the sample rate - a sinc with
bandwidth fs/2.  If you truncate a Fourier Series (which is a windowing in
frequency, not time) you'll see the Gibb's  phenomenon in time which is the
effect of a time convolution with a sinc.

In order to understand better it's a good idea to say whether things are
sampled or continuous and whether you're addressing a function in time or
frequency.


Fred

In article <YaOdnYnZtZLjf4LdRVn-gQ@centurytel.net>,
Fred Marshall <fmarshallx@remove_the_x.acm.org> wrote:
>> Calculating the reconstruction coefficients directly allows you to do
>> such things as using different asymmetric windows for the left and right
>> tails of the input data from the window used for the middle.
>
>I'm not sure what you mean: "Calculating the reconstruction coefficients
>directly".
>What process did you have in mind?  Isn't designing a reconstruction filter
>pretty direct?

The filter design is the same whether you precalculate or calculate the
coefficients as needed (at each FIR multiply for instance).  The usual
precalculation method assumes a time versus memory efficiency trade-off
which may or may not be important.

>What benefits do you see for asymmetrical interpolation filters?  They
>aren't linear phase and some applications might not prefer the phase
>distortion.

How else do you interpolate around the first and last data points of
a finite sequence?  Zero padding or reflecting the data is the same as
asymmetrically windowing or folding the filter, which might result in
a computational savings.  If the data or some of its derivative are not
continuous at its boundaries, then linear phase might not be the most
important interpolation criteria at these points.



IMHO. YMMV.
-- 
Ron Nicholson   rhn AT nicholson DOT com   http://www.nicholson.com/rhn/ 
#include <canonical.disclaimer>        // only my own opinions, etc.

"Ronald H. Nicholson Jr." <rhn@mauve.rahul.net> wrote in message
news:bvs8gs$ab4$1@blue.rahul.net...
> In article <YaOdnYnZtZLjf4LdRVn-gQ@centurytel.net>,
> Fred Marshall <fmarshallx@remove_the_x.acm.org> wrote:
> >> Calculating the reconstruction coefficients directly allows you to do
> >> such things as using different asymmetric windows for the left and
right
> >> tails of the input data from the window used for the middle.
> >
> >I'm not sure what you mean: "Calculating the reconstruction coefficients
> >directly".
> >What process did you have in mind?  Isn't designing a reconstruction
filter
> >pretty direct?
>
> The filter design is the same whether you precalculate or calculate the
> coefficients as needed (at each FIR multiply for instance).  The usual
> precalculation method assumes a time versus memory efficiency trade-off
> which may or may not be important.

OK - possible of course but not an approach that I've ever heard being
considered in practice.  Not that there couldn't be situations where memory
is indeed the driver.  I guess with very long filters and lots of time
......

>
> >What benefits do you see for asymmetrical interpolation filters?  They
> >aren't linear phase and some applications might not prefer the phase
> >distortion.
>
> How else do you interpolate around the first and last data points of
> a finite sequence?  Zero padding or reflecting the data is the same as
> asymmetrically windowing or folding the filter, which might result in
> a computational savings.  If the data or some of its derivative are not
> continuous at its boundaries, then linear phase might not be the most
> important interpolation criteria at these points.

OK.  Thanks.

Fred

Hello All,

First off, let me thank you all for your replies! 

I think I need to explain the problem at hand a little better.  I'm
reading in the data of a recorded .WAV file.  This file is a recording
of music played on a single string of a guitar.

Once I've read the data in, I segment the data into 256 sample arrays,
on each of which I perform autocorrelation.

Therefore I have, let's say "N", number of segments, each 256 samples
long, which I now want to interpolate.

The reason I'm not doing a pure mathematical interpolation is that
I've read that this leads to signal distortion by adding undesired
spectral images to the signal at multiples of the original sampling
rate.

Fred, thank you for the link to your application. However, as you've
mentioned in the README, filter lengths greater than 300 are not
suitable as inputs.

I basically need to code the interpolation in C#, and need a set of
filter coefficients.  More importantly I want to understand what
exactly is going on.

I've read Paul Embree's approach (C Language Algorithms for DSP).  He
says that the pass band should be 0 to (fs*L)/2*L (L = interp. factor)
and the stopband from fs*L/2*L to fs*L/2

My wav file is sampled at 48KHz.  Therefore, for a 4:1 interpolation,
the cutoff frequency (ideally) should be 24KHz.  The stopband should
extend from 24KHz to 96KHz.

The smallest power of 2 greater than 96KHz is 2^17 = 131072 (for the
IFFT length).  Should I scale 0 - 131072 to 0 - 1024 and take an IFFT
and use the coefficients?

Thank you all for your patience.

Roger Waters wrote:


> The reason I'm not doing a pure mathematical interpolation is that
> I've read that this leads to signal distortion by adding undesired
> spectral images to the signal at multiples of the original sampling
> rate.

Whoa! Right there! The samples are valid only if the sampled signal 
contains no frequencies as high as half the sample rate. (It's up to 
whoever does the sampling to ensure that that condition is met.) 
Spectral images at multiples of the sample rate can be filtered out; 
they are irrelevant.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Roger Waters wrote:


 > The reason I'm not doing a pure mathematical interpolation is that
 > I've read that this leads to signal distortion by adding undesired
 > spectral images to the signal at multiples of the original sampling
 > rate.


Whoa! Stop right there! The samples are valid only if the sampled signal
contains no frequencies as high as half the sample rate. (It's up to
whoever does the sampling to ensure that that condition is met.)
Spectral images at multiples of the sample rate can be filtered out;
they are irrelevant.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������