Changing the Sampling Rate

7.1 Rate Conversion by Discrete-Time Processing

Suppose a continuous-time signal $x_c(t)$ was sampled with period $T$ to give $x[n] = x_c(nT)$, but we now need the same underlying signal sampled at a different period $T_d$:

This is needed constantly in practice: audio moves between 48 kHz (studio) and 44.1 kHz (CD), biomedical recorders sample at 2000 Hz but analysis may only need 250 Hz, and a modem may need fractional rate offsets corrected. The question of this chapter: can we obtain $x_d[n]$ from $x[n]$ using only discrete-time operations?

7.1.1 The analog detour vs the DSP path

One conceptual approach is to go back to the analog world:

1. Reconstruct $x_c(t)$ from $x[n]$ with a D/A converter and reconstruction filter;
2. Resample $x_c(t)$ with an A/D converter at the new period $T_d$.

7.2 Downsampling: Rate Reduction by an Integer Factor M

7.2.1 The compressor in the time domain

The sampling-rate compressor（取樣率壓縮器）, written $\downarrow M$, keeps every $M$-th sample and discards the rest:

The output is exactly what we would have obtained by sampling $x_c(t)$ directly with the longer period $T' = MT$ — the sampling rate drops by $M$, the sample period grows by $M$. A decimator adds a protective lowpass filter in front of the compressor (next subsection):

7.2.2 Spectrum of the downsampled signal

To see what downsampling does in frequency, write $x_d[n]$ as a two-step operation: (1) multiply $x[n]$ by the periodic sampling train $c_M[n] = \frac{1}{M}\sum_{i=0}^{M-1} e^{j2\pi i n/M}$, which equals 1 when $n$ is a multiple of $M$ and 0 otherwise; (2) discard the zeros. Step (1) gives $\tilde{x}[n] = x[n]\,c_M[n]$, whose DTFT is an average of $M$ shifted copies of $X(e^{j\omega})$; step (2) rescales the frequency axis by $M$. The result:

Equivalently, in terms of the original continuous-time spectrum $X_c(j\Omega)$, the downsampled sequence looks exactly like sampling with period $MT$:

7.2.3 Aliasing and the anti-alias prefilter — the decimator

From Fig. 7-3, downsampling is alias-free only when the original spectrum leaves room to be stretched:

If $x[n]$ is not bandlimited to $\pi/M$, we must lowpass-filter first. The cascade of an ideal lowpass filter (gain 1, cutoff $\pi/M$) and the compressor $\downarrow M$ is called a decimator（十中抽一器／降頻器）, and the whole operation decimation（抽取）. The prefilter deliberately throws away the band $|\omega| > \pi/M$ so that the surviving band comes through clean.

7.3 Upsampling: Rate Increase by an Integer Factor L

7.3.1 The expander

To raise the sampling rate we want $x_i[n] = x_c(nT_i)$ with $T_i = T/L$ — i.e. we must create $L-1$ new samples between every pair of existing ones. The first stage is the sampling-rate expander（取樣率擴張器） $\uparrow L$, which inserts $L-1$ zeros between samples:

7.3.2 Spectrum: compression and images

The DTFT of the zero-stuffed sequence follows directly from the definition:

7.3.3 Ideal bandlimited interpolation

The ideal anti-imaging filter has gain $L$ and cutoff $\pi/L$:

Convolving $x_e[n]$ with $h_i[n]$ gives the discrete-time version of Shannon sinc interpolation (the formula on slide 17):

Two properties follow immediately, because $h_i[0]=1$ and $h_i[kL]=0$ for $k \ne 0$:

• Sample preservation: $x_i[kL] = x[k]$ — the original samples pass through unchanged.
• Exact reconstruction: if $x[n]=x_c(nT)$ came from a properly bandlimited $x_c(t)$ (i.e. $X_c(j\Omega)=0$ for $|\Omega| \ge \Omega_N$ and $T \le \pi/\Omega_N$), then $x_i[n] = x_c(nT/L)$ exactly — the interpolator behaves like an ideal C/D converter running $L$ times faster.

7.3.4 Practical interpolators: zero-order hold and linear

The ideal sinc filter is infinitely long and non-causal, so practical systems use short approximations. The two classics:

The ZOH simply repeats each sample $L$ times (a staircase); the linear interpolator draws straight lines between samples (note $h_{\mathrm{lin}} = \tfrac{1}{L}\, h_{\mathrm{zoh}} * h_{\mathrm{zoh}}$, a triangle). Their frequency responses,

only approximate the ideal brick wall: the passband droops slightly and the images around $2\pi i/L$ are attenuated but not eliminated (the linear interpolator at least places spectral zeros exactly at the image centers $\omega = 2\pi i /L$). They work well when the signal is heavily oversampled, i.e. its energy is concentrated at $|\omega| \ll \pi$.

7.4 General Rate Change by a Rational Factor L/M

7.4.1 Cascade and combined filter

A noninteger rate change is achieved by interpolating by $L$ first, then decimating by $M$, giving a new sampling period

In the cascade, the interpolation filter (gain $L$, cutoff $\pi/L$) and the decimation filter (gain 1, cutoff $\pi/M$) operate back-to-back at the same (high) rate, so they collapse into one lowpass filter taking the tighter of the two cutoffs:

7.4.2 Worked example: 3/2 rate change

7.5 MATLAB Practice: decimate, downsample, resample, rat

7.5.1 decimate vs downsample（有無抗混疊濾波的差別）

The lecture's running example: a signal sampled at 1000 Hz must be reduced to $1000/8 = 125$ Hz. Two one-line options exist, and they are not equivalent:

7.5.2 resample and rat（任意倍率重新取樣）

Y = resample(X, P, Q) resamples X at $P/Q$ times the original rate (polyphase implementation; length of Y is $\lceil P/Q \cdot \mathrm{length}(X)\rceil$). It applies an anti-aliasing lowpass FIR filter during the conversion and compensates the filter delay. [N, D] = rat(X, tol) returns integers with $|N/D - X| \le tol$ (continued-fraction truncation; default tolerance $10^{-6}\,\lVert X\rVert_1$) — the standard way to get $P$ and $Q$ from a ratio of sampling rates:

org_SR = 2000;                       % original rate (Hz)
new_SR3 = 300;                       % target rate (Hz)
[p, q] = rat(new_SR3/org_SR);        % p = 3, q = 20
new_signal3 = resample(org_signal, p, q);
new_taxis3 = (1:length(new_signal3))'/new_SR3;   % REBUILD the time axis!

7.6 Summary

Exercises

All five problems of ICE503 DSP Homework #7, restated with full worked solutions.

fs = 1000;  t = 0:1/fs:2-1/fs;
x  = sin(2*pi*20*t) + sin(2*pi*100*t);

for M = [5 8]                          % (a) M = 5, (b) M = 8
    x_down = downsample(x, M);         % bare compressor, no filter
    x_deci = decimate(x, M);           % Chebyshev LPF, then compressor
    t_M    = downsample(t, M);

    figure;
    subplot(2,1,1);
    plot(t, x, 'Color', [.75 .75 .75]); hold on;
    stem(t_M, x_down, '.'); xlim([0 0.25]);
    title(sprintf('downsample, M = %d  (fs'' = %g Hz)', M, fs/M));
    subplot(2,1,2);
    plot(t, x, 'Color', [.75 .75 .75]); hold on;
    stem(t_M, x_deci, '.'); xlim([0 0.25]);
    title(sprintf('decimate, M = %d', M));
    xlabel('time (s)');
end

load('EMG_ICA.mat');
org_signal = fdata(:,1);               % channel 1, fs = 2000 Hz
org_SR = 2000;
t_org  = (1:length(org_signal))'/org_SR;

% ---------- (a) 2000 -> 5000 -> 2000 ----------
[p1, q1] = rat(5000/org_SR);           % p1 = 5,  q1 = 2
up_sig   = resample(org_signal, p1, q1);     % now 5000 Hz
[p2, q2] = rat(org_SR/5000);           % p2 = 2,  q2 = 5
back_a   = resample(up_sig, p2, q2);         % back to 2000 Hz
t_a = (1:length(back_a))'/org_SR;

figure; plot(t_org, org_signal); hold on; plot(t_a, back_a);
legend('original 2000 Hz', 'up 5000 \rightarrow down 2000');

% ---------- (b) 2000 -> 300 -> 2000 ----------
[p3, q3] = rat(300/org_SR);            % p3 = 3,  q3 = 20
low_sig  = resample(org_signal, p3, q3);     % now 300 Hz
[p4, q4] = rat(org_SR/300);            % p4 = 20, q4 = 3
back_b   = resample(low_sig, p4, q4);        % back to 2000 Hz
t_b = (1:length(back_b))'/org_SR;

figure; plot(t_org, org_signal); hold on; plot(t_b, back_b);
legend('original 2000 Hz', 'down 300 \rightarrow up 2000');

fs = 100;  N = 100;
n  = 0:N-1;
x  = cos(2*pi*5*(n/fs));               % 1 s of a 5 Hz cosine, N = 100

k  = 0:N-1;
W  = exp(-1j*2*pi/N) .^ (k'*n);        % N-by-N DFT matrix, W(k+1,n+1) = e^{-j2pi kn/N}

X_mat  = W * x.';                      % (a) DFT via matrix
x_mat  = (W' * X_mat) / N;             % (b) IDFT via matrix  (W' = conjugate transpose)
X_fft  = fft(x).';                     % (c) DFT via fft
x_ifft = ifft(X_fft);                  % (d) IDFT via ifft

% (e) stem plots
figure;
subplot(2,2,1); stem(k, abs(X_mat));   title('(a) |X[k]|, DFT matrix');
subplot(2,2,2); stem(n, real(x_mat));  title('(b) x[n], IDFT matrix');
subplot(2,2,3); stem(k, abs(X_fft));   title('(c) |X[k]|, fft');
subplot(2,2,4); stem(n, real(x_ifft)); title('(d) x[n], ifft');

max(abs(X_mat - X_fft))    % ~1e-12  : matrix DFT == fft
max(abs(x_mat.' - x))      % ~1e-14  : perfect reconstruction