Integer upsampling - leave original samples unaltered

Question

I'm looking at the problem of integer upsampling (i.e. interpolation). This can be achieved by first increasing the sample rate by factor $L$ (i.e. inserting $L-1$ zeros between the original samples) and then low-pass filtering the signal (see e.g. Oppenheim/Schafer, "Discrete Time Signal Processing", Sect. 4.6.2).

When I compare the the samples of the original signal with the corresponding samples of the upsampled signal, I observe a slight difference in amplitude. Is there a way to avoid this and to leave the original samples unaltered? That's what I would expect of an interpolation method.

The filter I used is a symmetric FIR filter with an even number of coefficients and for comparison I compensate the delay introduced by that filter.

I know that there are other interpolation methods (linear, cubic, hermite, etc.), but I would like to use this particular method.

I hope I could make myself understood. Thanks in advance for your help!

Answer 1

You need to be careful with the delay compensation. If your symmetric filter has an even number of taps, then you get a non-integer delay. If you use a symmetric filter with an odd number of coefficients then the delay of $\frac{N-1}{2}$ samples is an integer number ($N$ is the number of filter coefficients).

Answer 2

this answer comes only from my direct experience:

first, if you are upsampling and leaving the original samples unaltered, you are necessarily using a windowed sinc() function as the interpolation kernel. even if your method does not ostensibly appear to be the windowing method (say you're using Lagrange or Hermite polynomials), you can take whatever resulting impulse response and point-by-point divide by the sinc(), deal with the 0/0 division as a "removable singularity" (use the limit), and the result will be your effective window. when looking at the interpolation FIR coefficients for the various phases and then look at the specific phase that corresponds to the original sample time, for windowed-sinc you will see 0s for every coefficient except one coef which has the value 1. that corresponds to "copying the original sample over".

my experience is that, for the windowed-sinc kernel, the Kaiser window is best, but using either Parks-McClellan (firpm) or Least Squares (firls) is better, but the output samples that lie on top of the original samples are not unaltered, since they have to be filtered like every other sample. so you have to compute them. but, comparing apples-to-apples, for the same stopband attenuation, the P-McC or LS will have a shorter impulse response and fewer samples are needed in the FIR computation.

if you're upsampling by exactly 2, then it makes sense to save 50% of the computation, and compute only the in-between samples and copy the original sample. the FIR will be a little longer with a windowed-sinc. but if you upsample by 4, windowed-sinc is not worth it. you may save 25% of the computation by copying the original, but the additional costs on the 3 other samples (compared to what you would have with P-McC or LS) exceeds what you save by not needing to compute the sample of the original.

Answer 3

You can design an $L$:th-band (also known as Nyquist) filter, i.e., an FIR filter where the mid-tap is $1/L$ and every $L$:th coefficient starting from the mid-tap is zero. This gives you a filter which is "complementary anti-symmetric" around $\pi/L$ and is suitable for upsampling/downsampling by $L$. (Note that you should multiply the output of the interpolation filter by $L$, so in principle you will get the mid-tap 1 and all other coefficients in the corresponding polyphase branch will be zero, and, hence, you keep the original value.)

The simplest case is $L=2$ which gives you a half-band filter, i.e., a filter with an impulse response looking like $h_0 0 h_2 0 h_4 1/2 h_4 0 h_2 0 h_0$. In addition, you will for half-band filters have $\omega_cT + \omega_sT = \pi$ (passband edge + stopband edge = $f_s/2$) and $\delta_c = \delta_s$ (same passband and stopband ripples). These properties are the complementary anti-symmetry around $\pi/2$, i.e., $|H(e^{j\omega T}) + H(e^{j(\pi-\omega T)})|=1$.

For the general case, the relation is slightly more involved, but there will be a relation between passband and stopband ripples and shifted versions of the frequency response will sum up to one.

To design the filter you can "trick" Remez/Parks-McClellan, at least in the $L=2$ case, by setting symmetric edges and identical ripples and, if needed, change every other coefficient to exactly 0 and 1/2. For the general case it works sometimes better, sometimes worse, so in general you may need to rely on either linear programming design or e.g. fminimax in Matlab.