# Upsampling with time offsets

Suppose I have done 4x oversampling for a continuous time signal, but the successive sampling times have a linearly increasing offset. Specifically, the samples with indices {4k; k=0, 1, 2,...} are correctly sampled at times {4kTs}, but the samples at indices {4k+1; k= 0, 1, 2,..} are sampled at {(4k+1)Ts+e}, the samples at indices {4k+2; k= 0, 1, 2,..} are sampled at {(4k+2)Ts+2e)}, and the samples with indices{4k+3; k= 0, 1, 2,..} are sampled at {(4k+3)Ts+3e)}. How does the spectrum of this oversampled signal look ?

• Can you choose one of the answers if satisfactory? Or ask if any clarifications needed. – DSP Rookie Apr 19 at 6:12
• how do I choose ? And what's the relevance of choosing ? – voy82 Apr 20 at 7:57
• You choose by accepting whichever answer satisfied your doubt, by clicking on the tick mark. If none of the answers satisfied your doubt then you comment on the answers what is lacking. – DSP Rookie Apr 20 at 8:12

Answer: You will see residual images of $$X(f)$$ at multiples $$f_s$$, $$2f_s$$ and $$3f_s$$, and distorted image of $$X(f)$$ at non-zero multiples of $$4f_s$$, when sampling in the manner you explained. Depending on value $$e$$, the size of residual will change. I have explained how in detail below.

Ideally, sampling at $$4f_s$$ would have completely cancelled those images of $$X(f)$$ at multiples of $$f_s$$, $$2f_s$$ and $$3f_s$$ and you would've seen images of $$X(f)$$ only at multiples of $$4f_s$$ and the magnitude scaled by $$4f_s$$.

(Explanation seems long only because I have included lots of pictures to show. Please follow through.)

I would like to give you an intuition on how to visualize sampling at any rate $$f_s$$. You probably have a pretty good idea about that. But then I would like to show pictorially what happens when you sample at $$2f_s$$ and then you can extend the idea to $$4f_s$$.

1. Sampling a bandlimited signal $$x(t)$$ at sampling rate $$f_s$$ :

When you sample $$x(t)$$ at sampling rate $$f_s$$, you are basically, multiplying $$x(t)$$ with a periodic pulse train with period $$T_s = \frac{1}{f_s}$$ in time domain. Hence, in frequency domain you see a convolution of $$X(f)$$ with Fourier representation of that periodic pulse train.

Since the pulse train is periodic, it's Fourier representation will be obtained by computing Fourier Series. The pulse train can be represented by: $$\sum^{\infty}_{k=-\infty} \delta(t - kT_s)$$ and it's Fourier transform as: $$\frac{1}{T_s}\sum^{\infty}_{k=-\infty} \delta(f - kf_s)$$

Notice that the magnitude of Fourier domain representation of the sampling pulse train has scaled up magnitude of $$\frac{1}{T_s} = f_s$$. Assume, frequency representation of $$x(t)$$ as below : Now, as described above, Sampling is nothing but multiplication of $$x(t)$$ with a $$T_s$$ periodic pulse train in time domain and hence a convolution of $$X(f)$$ with Fourier transform of the pulse train in frequency domain. Mathematically, $$x(kT_s) = x(t).\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$$ $$X_{sampled}(f) = X(f) * f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$ 1. Now, lets explore what happens when we sample at $$2f_s$$ or when the sampling pulse train becomes $$\frac{T_s}{2}$$ periodic :

The sampling pulse train becomes as shown below : The sampling pulse train can be represented by a sum of two $$T_s$$ periodic pulse trains:

1. One pulse train as if $$T_s$$ periodic and centered at 0.

2. Second pulse train as if $$T_s$$ periodic but shifted by $$\frac{Ts}{2}$$.

Hence, mathematically it will be as follows: $$\sum^{\infty}_{k=-\infty}\delta(t-k\frac{T_s}{2}) = \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2})$$ Hence, the Fourier Transform of the sampling pulse train will also sum of these 2 trains, as convolution is a linear operation. Also, use the time shift property of Fourier transform to get the result as follows:

$$\mathcal F \{ \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2}) \}$$ $$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.e^{-j\pi \frac{f}{f_s}} \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.(cos(\pi\frac{f}{f_s})-\mathbb i.sin(\pi \frac{f}{f_s})). \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

Notice that the sum is evaluated only at integral multiples of $$f_s$$, because of the $$\delta(f - kf_s)$$. What this means is that $$sin(\pi \frac{f}{f_s})$$ will always be $$0$$, so, no imaginary images of $$X(f)$$ will be seen, and $$cos(\pi \frac{f}{f_s})$$ will be $$1$$ at even multiples of $$f_s$$ and $$-1$$ at odd multiples of $$f_s$$. Pictorially, the fourier transform of pulse train which is used to sample at $$2f_s$$ will look like following : So, the even multiples of $$f_s$$ will be doubled in magnitude to $$2f_s$$ and odd multiples of $$f_s$$ will cancel out each other to cancel the images of $$X(f)$$. This is the reason you see images of $$X(f)$$ only at multiples of $$2f_s$$ when sampling at double the rate, because images of $$X(f)$$ at odd multiples of $$f_s$$ cancel each other out.

Now, consider the case which you have explained in your question. When you break your pulse train into 4 pulse trains which are $$T_s$$ periodic individually, but shifted as below:

1. $$\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$$

2. $$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{T_s}{4} - e)$$

3. $$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{2T_s}{4} - 2e)$$

4. $$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{3T_s}{4} - 3e)$$

Conclusion:

When you check their Fourier transforms, you will find that images at multiples of $$f_s$$, $$2f_s$$ and $$3f_s$$ will not get cancelled completely because the negative impulses (both real and imaginary) are shifted by $$e$$, $$2e$$ and $$3e$$ respectively. And, image at multiples of $$4f_s$$ will also not be aligned exactly to give a scaling of $$4f_s$$ but they will be fudged around to give a distorted image of $$X(f)$$ except at $$k=0$$, that is the original image of $$X(f)$$ centered around DC.

Depending upon the value of e, the real and imaginary images of $$X(f)$$ will have residuals at $$f_s$$, $$2f_s$$ and $$3f_s$$, and images at non-zero multiples of $$4f_s$$ will be fudged around.

• Fun-fact: you can still reconstruct the original x(t), because X(f) at DC is still untouched. Meaning, you can still throw 3 in 4 samples and interpolate to get x(t) back. – DSP Rookie Apr 7 at 13:40

The spectrum of the oversampled signal as you describe it would look very odd and is likely not what you are intending to do. Let me explain:

First consider each group separately, instead representing what would be the decimated signal if you were sampling at the higher 4x rate as $$f_s$$: Each of the four spectrums would span from frequency = $$0$$ to $$f_s/4$$ and contain all energy within that frequency range as well as all the aliasing from $$f_s/4$$ to $$f_s$$. Each would be translated by a different phase slope in frequency consistent with the time delay each shift represents. First consider the unit delays only and assume the additional time offset $$e = 0$$ to best understand the effect of the unit delays, then we can add in the effect of non-zero $$e$$.

Each delay of one sample would add a linear phase in frequency extending negatively from $$0$$ to $$2\pi$$ as the the normalized radian frequency goes from $$0$$ to $$2\pi$$ (Meaning as the frequency goes from $$0$$ to $$f_s$$).

Consider the z-transform of an m sample delay:

$$\sum_{n=0}^{N-1}x[n-m]z^{-n} = X(z)z^{-m}$$

And the DFT is simply the z transform with z limited to the unit circle; $$z=e^{j2\pi k/N}$$ as k = $$0$$ to $$N-1$$

So here we see that for each delay m, the DFT would be $$X(k)e^{-mj2\pi k/N}$$, with the phase negatively increasing to $$2\pi$$ for m =1, to $$4\pi$$ for m=2, etc. Thus the higher frequencies that are above $$f_s/4$$ are aliased into the $$0$$ to $$f_s/4$$ spectrum with a different phase slope for each of the four groups.

As far as the additional time offsets, which may be fractional samples, consider the Fourier Transform of a time delay given as

$$\mathscr{F}\{x(t-\tau)\} = e^{-\tau}$$

And we see that the time offsets just add an additional slope to the phase slopes given by the unit delays. The time offset as a fraction of a sample given as $$d$$ would introduce a phase slope according to $$z^{-d}$$ using the method above with limiting z to the unit circle.

Importantly, you cannot recreate the interpolated spectrum simply by cascading the samples in frequency as you are describing. As I explained above above you would only be showing the decimated and folded spectrum extending from $$0$$ to $$f_s/4$$, just with four phases of the same bin next to each other. Combining the interpolated spectrum is more complicated as it involves summing all four spectrums with the appropriate phase slope to compensate for the phase introduced as described above.

If your signal is bandlimited to $$[-f_0/2,+f_0/2]$$, the digitized signal sampled at $$f_s$$ can be considered as a summation of 4 different sampled signals with offset of $$0$$, $$e$$, $$2e$$, $$3e$$. Assuming $$3e \le 1/f_s$$.

First imagine $$e=0$$. This means, you have perfectly oversampled 4x. That is $$f_s = 4 \times f_0$$. It is like you have sampled the signal 4 times with sample rate of $$f_0$$, each of them with offset of $$k/(4f_0)$$, where $$k \in \{0,1,2,3\}$$, interleaved them with 3 zeros and added them. Assume $$T_s = 1/f_s = (1/(4f_0))$$. This is equivalent to
$$x_1 = x(0), 0 ,0 ,0, x(4T_s), 0, 0, 0..\\ x_2 = 0 ,x(T_s),0,0,0, x(5T_s),0,0,0.. \\ x_3 = 0 ,0,x(2T_s),0,0,0, x(6T_s),0,0,0.. \\ x_4 = 0 ,0,0,x(3T_s),0,0,0, x(7T_s),0,0,0.. \\ x = x_1 +x_2+x_3+x_4$$ Equivalently in fourier domain $$X(e^{j\omega}) = X_1(e^{j4\omega})+X_2(e^{j4\omega})e^{-j\omega}+X_3(e^{j4\omega})e^{-j2\omega}+X_4(e^{j4\omega})e^{-3\omega}$$ which is same as if you had sampled original signal at $$f_s$$.(This is proved later using MATLAB example for 2x oversampled case). If there were offsets of $$e,2e,3e$$ as you mentioned in each of those copies,

$$x_1 = x(0), 0 ,0 ,0, x(4T_s), 0, 0, 0..\\ x_2 = 0 ,x(T_s+e),0,0,0, x(5T_s+e),0,0,0.. \\ x_3 = 0 ,0,x(2T_s+2e),0,0,0, x(6T_s+2e),0,0,0.. \\ x_4 = 0 ,0,0,x(3T_s+3e),0,0,0, x(7T_s+3e),0,0,0.. \\ \hat{x} = x_1 +x_2+x_3+x_4$$

Equivalently in fourier domain $$\hat{X(e^{j\omega})} = X_1(e^{j4\omega})+X_2(e^{j4\omega})e^{-j\omega(1+e)}+X_2(e^{j4\omega})e^{-j2\omega(1+e)}+X_4(e^{j4\omega})e^{-3\omega(1+e)}$$

APPENDIX: To show that if $$x[n]$$ is an $$N$$ sample signal sampled at rate $$f_s$$, it can be shown that it is sum of 2 signals sampled at rate $$f_s/2$$ but with offset of $$1$$ sample and interleaved with zeros and added. That is $$X(e^{j\omega}) = X_1(e^{j2\omega})+X_2(e^{j2\omega})e^{-j\omega}$$.

That is $$e=0$$ in the below code, the spectrum of $$x$$ and $$X_1(e^{j2\omega})+X_2(e^{j2\omega})e^{-j\omega}$$ would be same. If $$e=0.25$$, spectrum of $$x$$ would not be same as $$X_1(e^{j2\omega})+X_2(e^{j2\omega})e^{-j\omega(1+e)}$$.

clc
clear all
close all

x=randn(1,16)+1i*randn(1,16);
xph1 = x(1:2:end);
xph2 = x(2:2:end);
xph1ups=upsample(xph1,2);
xph2ups=upsample(xph2,2);

e=0.25;
F1= fft(xph1ups);
F2= fft(xph2ups);
F=F1+exp(-1i*2*pi/16*(0+e:1:15+e)).*F2;

plot(1:16,abs(fft(x)),'b',1:16,abs(F),'r')

• @voy82 They are not the same. You cannot take out the exponentials as common factor. You can check that with the MATLAB script I uploaded. It is not same as $X_1(e^{j2\omega})(1+e^{-j\omega})$ – jithin Apr 7 at 5:40
• thanks ! I was also thinking along same lines. Since X1, X2, X3, X4 in my question represent the spectrums of Nyquist sampled x(t), are they not the same ? So, we should be able to take them out, so as to get X(e^(j4ω)) x (1+exp^(-jw(1+e)) + exp^(-j2w(1+e))+exp^(-j3w(1+e))) – voy82 Apr 7 at 5:40
• hmm.. i think, using your notations above, we should get : X(ej4ω)+X(ej4ω) e−j4ω/(1+e) + X(ej4ω) e−j8ω/(1+e) + X(ej4ω) e−j12ω/(1+e); where X(ejω) is the spectrum of the Nyquist sampled original continuous signal x(t). Would you agree ? I am thinking if for e=0, this gives X(ej4ω). – voy82 Apr 7 at 6:00
• sorry I did not mean to have the "/" there in the "e−j4ω/(1+e)" like terms. I was basically repeating my first comment but with the 4ω, 8ω, 12ω instead of ω, 2ω, 3ω. But you clearly do not seem to agree with it (that we can take the X(ej4ω) out) ! Looking at your matlab code now. thanks ! – voy82 Apr 7 at 6:13
• Repeating my comment because of formatting issue earlier : @voy82 No the (1+e) terms are not in denominator. Sorry if my MATLAB code appears like that. The (1+e) terms are all in numerator. Also for the term $X_1(e^{j4\omega)}$ , the $4\omega$ cannot be taken outside. So you would not get 4$\omega$,8$\omega$,12$\omega$ – jithin Apr 7 at 6:14