Forums

FIR Interpolation question

Started by Roger Waters February 2, 2004
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