FIR Interpolation question

This question relates to the Multirate FAQ at the dspGuru site:
http://dspguru.com/info/faqs/multrate/interp.htm

I have a 256 element array that I want to interpolate at a variable
density, between 4 and 10.

Let's take the first case, i.e. the interpolation factor or density, L
= 4.

Thus, I need a 256 * 4 = 1024 tap filter.

My question is this: How do I obtain the coefficients of this filter

"Roger Waters" <rogerwaters@fastmail.fm> wrote in message
news:6975d4b5.0402021625.710dad46@posting.google.com...
> This question relates to the Multirate FAQ at the dspGuru site:
> http://dspguru.com/info/faqs/multrate/interp.htm
>
> I have a 256 element array that I want to interpolate at a variable
> density, between 4 and 10.
>
> Let's take the first case, i.e. the interpolation factor or density, L
> = 4.
>
> Thus, I need a 256 * 4 = 1024 tap filter.
>
> My question is this: How do I obtain the coefficients of this filter

Roger,

Is the "array" a filter?  Or, is it more like "data"?

You can certainly interpolate 256 data points with a filter that's shorter
than 256 taps - in fact, to avoid the transient response, it must be shorter
than the length of the data.

If the array is a filter that you want to interpolate, then yes, you'd end
up with 1024 taps.

If you're multiplying in the frequency domain then you need the data and the
filter to be the same length before computing the Discretee Fourier
Transform or FFT.  That length is at least M+L-1 where M and L are the
lengths of the data and the filter respectively.
Note that M is the data length after being stuffed with zeros between the
original data points.

Below, I'm discussing interpolation by a factor of 2 that can be repeated to
get a factor of 4.

Zero stuff 1 for 1 for a factor of 2
Filter
Done
Repeat for as many factors of 2 are necessary.


You may want to use a halfband filter so that the original samples are
unchanged.  You can find one for a PC at:

ftp://ftp.mission-systems-inc.com/outgoing/Halfband/

You can try various lengths to see what makes sense for your application.
Then, there are fast ways to do it as in the fft/multiply/ifft process
mentioned above.
[Note: a "halfband" filter is a special lowpass filter that limits the
output to fs/4 with is half of fs/2 and where fs means the output sample
rate].

Regarding the transient response, the halfband filter is symmetrical and of
odd length, has the center 3 coefficients nonzero and then has zeros
alternating with nonzero coefficients.  Unless you want to consider zero
coefficients at the ends, which really makes no sense, then the length of
the filter will be 3+4*z where z is the number of zero coefficients on
either side of the center.
A filter with 1 zero on each side will have length 7 and 5 nonzero
coefficients.
A filter with 2 zeros on each side will have length 11 and 7 nonzero
coefficients.
etc.

The transient response will persist for L-1 samples at both ends of the
interpolation.  Another way to say this is that there won't be valid
interpolation points until the filter fully overlaps one end of the data -
that will generate the first valid output.  So, you throw away 2*(L-1) data
points, (L-1) at each end of the output.

If you must have a filter that does 4x interpolation in one step, then you
will start by zero padding the data 3 for 1 to get a factor of 4.
The corresponding filter will be a quarter band filter with cutoff at fs/8
and will have every fourth coefficient of zero with the exception of the
center 7 coefficients - to preserve the original samples - which may or may
not be necessary.

Fred

In article <6975d4b5.0402021625.710dad46@posting.google.com>,
Roger Waters <rogerwaters@fastmail.fm> wrote:
>This question relates to the Multirate FAQ at the dspGuru site:
>http://dspguru.com/info/faqs/multrate/interp.htm
>
>I have a 256 element array that I want to interpolate at a variable
>density, between 4 and 10.
>
>Let's take the first case, i.e. the interpolation factor or density,
>L = 4.
>
>Thus, I need a 256 * 4 = 1024 tap filter.
>
>My question is this: How do I obtain the coefficients of this filter

If the array consists of filter coefficients for a lowpass FIR filter
with an cutoff frequency of just the right value, you might be able to
use the array to interpolate itself for the intermediate points.

There's another thread in this forum which suggests that a windowed
sinc reconstruction using just a few nearby points might be as good
as any common polynomial interpolation scheme.


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

In article <h7idnfJN56Va0ILdRVn-jw@centurytel.net>,
Fred Marshall <fmarshallx@remove_the_x.acm.org> wrote:
>Below, I'm discussing interpolation by a factor of 2 that can be repeated to
>get a factor of 4.
>
>Zero stuff 1 for 1 for a factor of 2
>Filter
>Done
>Repeat for as many factors of 2 are necessary.

Why zero pad?  These points contribute nothing to the result.

Just calculate the (windowed) sinc reconstruction coefficients directly
for each output point.  If the phase is repeated, which will happen every
4th point in 4-to-1 upsampling, then you can save and reuse coefficients
for repeated phases (inter-input-point fractional values).

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.


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

> Let's take the first case, i.e. the interpolation factor or density, L
> = 4.
> 
> Thus, I need a 256 * 4 = 1024 tap filter.
> 
> My question is this: How do I obtain the coefficients of this filter

You dont need 1024 taps, you can use any number of taps you like, with
a tradeoff between performance and complexity. When you insert zeros
at the new sampling instants (in this case three zeros between each
sample), you will get a squeezing of the signals periodic spectrum (in
this case the spectrum of the signal originally can occupy anywhere in
the range [0,pi) and after inserting zeros, will be repeated 4 times
in the range [0,pi). The filter you need removes the three unwanted
images, which will occur above the frequence pi/4.

To conclude, you need a low pass filter with cutoff frequency pi/4.
You can design it any way you like, even using IIR. The ideal
interpolation filter is a brickwall filter. Its impulse response is
IIR (sinc shaped). If you truncate this to 21 taps for example, you
will still have a lowpass filter with the right cutoff, but the
spectral response will no longer be ideal (it will in fact be sinc
shaped). I like to apply a Kaiser window to this to improve spectral
performance (more stopband attenuation). The advantage ofthis method
(windowed sinc pulse) is that the filter design is straight forward
and the response is linear phase. Furthermore, the Matlab function
FIR1 will generate the coefficients for you with no effort, so you can
fiddle with the number of coefficients and see which filter best suits
your needs.

i did a software for coefficients calculation:
http://www.winfilter.150m.com

Roger Waters wrote:

> This question relates to the Multirate FAQ at the dspGuru site:
> http://dspguru.com/info/faqs/multrate/interp.htm
> 
> I have a 256 element array that I want to interpolate at a variable
> density, between 4 and 10.
> 
> Let's take the first case, i.e. the interpolation factor or density, L
> = 4.
> 
> Thus, I need a 256 * 4 = 1024 tap filter.
> 
> My question is this: How do I obtain the coefficients of this filter

"Ronald H. Nicholson Jr." <rhn@mauve.rahul.net> wrote in message
news:bvnkkk$67c$2@blue.rahul.net...
> In article <h7idnfJN56Va0ILdRVn-jw@centurytel.net>,
> Fred Marshall <fmarshallx@remove_the_x.acm.org> wrote:
> >Below, I'm discussing interpolation by a factor of 2 that can be repeated
to
> >get a factor of 4.
> >
> >Zero stuff 1 for 1 for a factor of 2
> >Filter
> >Done
> >Repeat for as many factors of 2 are necessary.
>
> Why zero pad?  These points contribute nothing to the result.
>
> Just calculate the (windowed) sinc reconstruction coefficients directly
> for each output point.  If the phase is repeated, which will happen every
> 4th point in 4-to-1 upsampling, then you can save and reuse coefficients
> for repeated phases (inter-input-point fractional values).

Ronald,

This is a conceptual framework within which one can comprehend the spectral
operations as well as the temporal.  It's easy to understand.  As to
efficient implementation, I figure that's another step.  You'll notice I
didn't suggest multiplying by zero or doing additon operations with zeros -
and I didn't warn against doing those things either.

I think it's a lot easier to go from description to implementation in the
process I described.

>
> 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?
What benefits do you see for asymmetrical interpolation filters?  They
aren't linear phase and some applications might not prefer the phase
distortion.

Fred

rogerwaters@fastmail.fm (Roger Waters) wrote in message news:<6975d4b5.0402021625.710dad46@posting.google.com>...
> This question relates to the Multirate FAQ at the dspGuru site:
> http://dspguru.com/info/faqs/multrate/interp.htm
> 
> I have a 256 element array that I want to interpolate at a variable
> density, between 4 and 10.
> 
> Let's take the first case, i.e. the interpolation factor or density, L
> = 4.
> 
> Thus, I need a 256 * 4 = 1024 tap filter.
> 
> My question is this: How do I obtain the coefficients of this filter

You can also upsample/resample in the frequency domain by doing a
larger ifft than fft.

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.

0. Do the N point fft of the input sequence.  In your case N=256

1. zero pad the frequency response (in the middle) to get the output
length you want, in your case 1024.  The Nyquist bin (Nfft/2), if
present, will need to be split into two bins.
2. inverse fft the zero padded frequency response to get the resampled
sequence.


In your example:
Fx = fft(x)
Fxpad = [ Fx(1:127) Fx(128)/2 zeros(1,767) Fx(128)/2 Fx(129:end) ]
x_up4 = ifft( Fxpad )

The x_up4 may have to be multiplied by the resampling factor (i.e. 4)
depending on your fft implementation.

- Mark Borgerding

On 3 Feb 2004 11:45:20 -0800, mark@borgerding.net (Mark Borgerding)
wrote:

>rogerwaters@fastmail.fm (Roger Waters) wrote in message news:<6975d4b5.0402021625.710dad46@posting.google.com>...
>> This question relates to the Multirate FAQ at the dspGuru site:
>> http://dspguru.com/info/faqs/multrate/interp.htm
>> 
>> I have a 256 element array that I want to interpolate at a variable
>> density, between 4 and 10.
>> 
>> Let's take the first case, i.e. the interpolation factor or density, L
>> = 4.
>> 
>> Thus, I need a 256 * 4 = 1024 tap filter.
>> 
>> My question is this: How do I obtain the coefficients of this filter
>
>You can also upsample/resample in the frequency domain by doing a
>larger ifft than fft.
>
>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.
>
>0. Do the N point fft of the input sequence.  In your case N=256
>
>1. zero pad the frequency response (in the middle) to get the output
>length you want, in your case 1024.  The Nyquist bin (Nfft/2), if
>present, will need to be split into two bins.
>2. inverse fft the zero padded frequency response to get the resampled
>sequence.
>
>
>In your example:
>Fx = fft(x)
>Fxpad = [ Fx(1:127) Fx(128)/2 zeros(1,767) Fx(128)/2 Fx(129:end) ]
>x_up4 = ifft( Fxpad )
>
>The x_up4 may have to be multiplied by the resampling factor (i.e. 4)
>depending on your fft implementation.
>
>- Mark Borgerding

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-]

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?


Thanks for your reply,
- Mark Borgerding