SPTK Addendum: Problems with resampling using MATLAB’s resample.m

Previous SPTK Post: Echo Detection Next SPTK Post: The Laplace Transform

In this brief Signal Processing Toolkit note, I warn you about relying on resample.m to increase the sampling rate of your data. It works fine a lot of the time, but when the signal has significant energy near the band edges, it does not.

I use MATLAB a lot for the CSP Blog because I think it produces nice graphics. I can make all kinds of plots with fine control over the details, and I can easily make videos (see here and here for examples) that show how functions or estimates change over time or over selections from some parameter set.

I also use MATLAB for my non-CSP-Blog signal-processing day-job work on some projects, but mostly I use C and gnuradio these days.

Sometimes I need to resample a discrete-time signal. Perhaps the signal was provided to me by someone else, or I generated it, or more likely it is derived from some other signal through something like spectral segmentation or other kinds of filtering. In any case, I’ve got some signal $x(k)$ sampled at some physical sampling rate $f_s$ and I want to resample it so that the sample rate of the new signal is $F_s > f_s$ . One common reason for this is that I am going to subject the signal to a nonlinear operation, which generally increases bandwidth, and also can produce finite-strength additive sine-wave components. I don’t want the frequencies of those sine waves (the cycle frequencies of the signal) to alias, so I have to keep them small by increasing the sampling rate. Physically, the cycle frequencies don’t change! I’m just talking about the smallness relative to the sampling rate.

Resampling in MATLAB

In MATLAB, I have a variable $x$ that is a complex vector having, say, $N$ samples. I want a new variable $y$ that is also a complex vector, but corresponds to a resampling of $x$ such that there are now plausible signal values inserted between each adjacent pair of elements of $x$ . This is also called interpolation, which is the opposite of decimation, where we wish to decrease the sampling rate, not increase it.

MATLAB provides a function called resample.m to do both kinds of sample-rate conversion (interpolation and decimation). Here is the usage message you get when you type ‘help resample’ in MATLAB version R2022b Rev 3:

>> help resample
 resample  Resample uniform or nonuniform data to a new fixed rate.
    Y = resample(X,P,Q) resamples the values, X, of a uniformly sampled
    signal at P/Q times the original sample rate using a polyphase
    antialiasing filter. If X is a matrix, then resample treats each
    column as an independent channel.
 
    In its filtering process, resample assumes the samples at times before
    and after the given samples in X are equal to zero. Thus large
    deviations from zero at the end points of the sequence X can cause
    inaccuracies in Y at its end points.
 
    [Y,Ty] = resample(X,Tx) resamples the values, X, of a signal sampled at
    the instants specified in vector Tx. resample interpolates X linearly
    onto a vector of uniformly spaced instants, Ty, with the same endpoints
    and number of samples as Tx.  Tx may either be a numeric vector
    expressed in seconds or a datetime object.  NaNs and NaTs (for datetime
    objects) are treated as missing data and are ignored.
 
    [Y,Ty] = resample(X,Tx,Fs) uses interpolation and an anti-aliasing
    filter to resample the signal at a uniform sample rate, Fs, expressed
    in hertz (Hz).
 
    [Y,Ty] = resample(X,Tx,Fs,P,Q) interpolates X to an intermediate
    uniform grid with sample rate equal Q*Fs/P and filters the result
    using UPFIRDN to upsample by P and downsample by Q.  Specify P and Q
    so that Q*Fs/P is least twice as large as the highest frequency in the
    input signal.

    [Y,Ty] = resample(X,Tx,...,METHOD) specifies the interpolation method.
    The default is linear interpolation.  Available methods are:
      'linear' - linear interpolation
      'pchip'  - shape-preserving piecewise cubic interpolation
      'spline' - piecewise cubic spline interpolation
 
    Y = resample(...,P,Q,N) uses a weighted sum of 2*N*max(1,Q/P) samples
    of X to compute each sample of Y.  The length of the FIR filter
    resample applies is proportional to N; by increasing N you will get
    better accuracy at the expense of a longer computation time.
    resample uses N = 10 by default. If N = 0, resample performs
    nearest neighbor interpolation: the output Y(n) is
    X(round((n-1)*Q/P)+1), with Y(n) = 0 for round((n-1)*Q/P)+1 >      length(X)).
 
    Y = resample(...,P,Q,N,BTA) uses BTA as the BETA design parameter for
    the Kaiser window used to design the filter.  resample uses BTA = 5 if
    you don't specify a value.
 
    Y = resample(...,P,Q,B) uses B to filter X (after upsampling) if B is a
    vector of filter coefficients.  resample assumes B has odd length and
    linear phase when compensating for the filter's delay; for even length
    filters, the delay is overcompensated by 1/2 sample.  For non-linear
    phase filters consider using UPFIRDN.
 
    [Y,B] = resample(X,P,Q,...) returns in B the coefficients of the filter
    applied to X during the resampling process (after upsampling).
 
    [Y,Ty,B] = resample(X,Tx,...) returns in B the coefficients of the
    filter applied to X during the resampling process (after upsampling).
 
    Y = resample(...,'Dimension',DIM) specifies the dimension, DIM,
    along which to resample an N-D input array. If DIM is not specified,
    resample operates along the first array dimension with size greater
    than 1.

    [Y,B] = resample(TT,...) resamples the data in timetable TT and returns
    a timetable Y. TT must contain double-precision data and must have at
    least two rows. Each variable in TT is treated as an independent
    signal. If TT has an N-D array as a variable, then resample operates
    along the first dimension. In other words, it treats columns as
    channels. The RowTimes in TT may either be a duration vector or a
    datetime object with unique and finite values. If RowTimes in TT are
    not sorted, then resample sorts them in an ascending order. Non-finite
    time values (NaNs and NaTs) are treated as missing data and are
    ignored. You can replace X,Tx with a nonuniformly sampled timetable TT
    in the above syntaxes that use X and Tx as inputs. You can replace X
    with a uniformly sampled timetable TT in all the other syntaxes. Use
    isregular to check if TT is uniformly sampled or not.

So the basic usage is to select $P$ and $Q$ , the over- and undersampling parameters, such that $F_s = (P/Q)f_s$ . If you want a sampling rate that is ten times the input-signal rate, pick $P=10$ and $Q=1$ . But you could also pick $P=20$ and $Q=2$ .

In this usage message, there is one warning. It concerns the edges of the input $x$ . If the signal starts, relative to zero, too large or ends, relative to zero, too large, the endpoints of the output $y$ might exhibit larger error than interior points in $y$ . There is no warning about the arrangement of the signal’s power as a function of frequency; no warnings about acceptable or problematic functional forms for power spectral densities. But the output is not as desired if the power spectrum has significant values near $f_s/2$ or $-f_s/2.$

Resampling with the Fourier Transform

An alternative to MATLAB’s resample.m is to use a form of zero padding to provide linearly interpolated samples between each pair of samples in the input $x$ . The recipe is to Fourier transform $x$ to obtain $X(j)$ for $j = 0, 1, \ldots, N-1$ , then prepend and append zeros to $X$ to form $Y$ , and finally to inverse Fourier transform $Y$ to obtain the resampled signal $y$ .

Let’s choose $P = 2^k$ , where the integer $k > 1$ . A typical choice might be $P = 2$ or $P=4$ . Then the vector $Y$ is defined by

$\displaystyle Y = [\underbrace{0\ 0\ \ldots 0}_{N(P-1)/2\ zeros} \ \underbrace{X}_{N\ values} \ \underbrace{0\ 0\ \ldots 0}_{N(P-1)/2\ zeros}] \hfill (1)$

so that the vector $Y$ has total length $NP$ ; it is $P$ times as long as $x$ and $X$ .

You can work out the interpolation effect by writing down the inverse discrete Fourier transform of $Y$ –but remember this is an inverse transform of length $NP$ , not $N$ .

To see why the FFT/IFFT method works intuitively, consider what the Fourier transform of a narrowband signal looks like. To do that, you need to grasp what narrowband means in that last sentence. Narrow relative to what? Relative to the sampling rate. Figure 1 shows a typical narrowband signal. The bandwidth of the signal (however you define it, within reason) is much smaller than the sampling rate of unity.

Figure 1. Power spectral density of a generic narrowband signal. The occupied bandwidth of the signal is much smaller than the sampling rate of 1.0. The characteristic appearance is that the energy of the signal is concentrated at midband. In the absence of noise, the signal spectrum would be surrounded by long sequences of zeros. Note the similarity between that structure and the structure of (1).

So if we pretend that our signal $x$ is a narrowband signal, we can imagine its Fourier transform will be an energetic central portion, surrounded by zeros on each side. By using this form of zero padding, we’re essentially forcing the transform to have the appearance it would have if it really were narrowband.

Examples

I consider two examples. In the first example, both resample.m and CSP-Blog FFT/IFFT work as desired to increase the sampling rate without introducing significant distortion. In the second, only the FFT/IFFT method works well.

The scenario consists of two simulated square-root raised-cosine BPSK signals and a small amount of noise. In the first example, the center frequencies of the two signals are $\pm 0.35$ . The symbol rates are both $1/15$ and the excess bandwidths are $0.5$ , so that each signal has an occupied bandwidth of $0.1$ . The number of samples is $N = 32768$ . A power spectrum estimate of the scene is shown in the upper plot of Figure 2.

The resample.m command is

outdata = resample(indata, P, Q);

where $P=4$ and $Q=1$ . The FFT/IFFT code is

pad_len = (floor((P-Q)/2)*N + N/2);
X = fftshift(fft(indata)).';
Y = [zeros(1, pad_len) X zeros(1, pad_len)];
y = ifft(fftshift(Y));

and the results are shown in the lower two plots in Figure 2. Everything is great!

Figure 2. Illustration of sample-rate conversion (increasing the sampling rate by a factor of four) for an input signal $x$ that contains two BPSK signals with center frequencies not particularly close to the sampling-band edges of -1/2 and 1/2. Note the high correspondence between the outputs of resample.m and the FFT/IFFT method.

In the next example, all variables are the same as in the first except the carrier frequencies are now $\pm 0.45$ instead of $\pm 0.35$ . A power spectrum estimate for this scene is shown in the upper plot of Figure 3.

The two resampling attempts use exactly the same code as before, and the results are shown in the lower two plots of Figure 3. In this case, there is extra energy, from aliased components of the signals, spectrally adjacent to the properly resampled signals. Because the power of the signal at $fc = -0.45$ is twice that of the signal at $fc = 0.45$ , you can see that the energy near the resampled BPSK signal at $0.45$ must be due to the signal at $-0.45$ and vice versa.

Figure 3. Illustration of sample-rate conversion (increasing the sampling rate by a factor of four) for an input signal $x$ that contains two BPSK signals with center frequencies that are close to the sampling-band edges of -1/2 and 1/2. Note the low correspondence between the outputs of resample.m and the FFT/IFFT method. resample.m produces aliases near the signals, and these aliases have significant energy.

Significance of Incorrect Resampling in CSP

The resampled signal with the aliases (Figure 3) will not give the same results when CSP estimators and detectors are applied as the original signal or the properly sampled signal (Figure 2). That is, the introduction of the aliased signal components, having significant energy relative to the original signals, will give rise to additional cycle frequencies.

So just be careful with resample.m. MATLAB has other functions that can be used to resample signals, which may or may not work better for the band-edge case. Let me know in the Comments if you have experience getting around the resample.m problem.

Previous SPTK Post: Echo Detection Next SPTK Post: The Laplace Transform

Author: Chad Spooner

I'm a signal processing researcher specializing in cyclostationary signal processing (CSP) for communication signals. I hope to use this blog to help others with their cyclo-projects and to learn more about how CSP is being used and extended worldwide. View all posts by Chad Spooner

2 thoughts on “SPTK Addendum: Problems with resampling using MATLAB’s resample.m”

Fabio Casagrande Hirono says:

March 30, 2023 at 11:45 am

Thanks for the post Chad, I was not aware of this limitation with in-built resampling function in Matlab. Coincidentally, a recent paper was published on the Journal of the Audio Engineering Society regarding the FFT method for sample rate conversion [1]. I haven’t read it in detail, but it seems to describe your proposed method and suggest a few tweaks – namely, spectral tapering and low-pass filtering to avoid possible ringing at the Nyquist frequency. It doesn’t contain any cyclic frequencies, though!

[1] V. Valimaki and S. Bilbao, “Giant FFTs for Sample-Rate Conversion”, J. Audio Eng. Soc., vol. 71, no. 3, pp. 88–99, (2023 March). DOI: https://doi.org/10.17743/jaes.2022.0061

Loading...

1. Chad Spooner says:
  
  March 31, 2023 at 5:40 pm
  
  Thank you! I don’t take any special action to reduce ringing or other distortions because I haven’t seen the need. But maybe that is just my inability to see clearly!
  Great tip!
  
  Loading...