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 sampled at some physical sampling rate and I want to resample it so that the sample rate of the new signal is . 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 that is a complex vector having, say, samples. I want a new variable that is also a complex vector, but corresponds to a resampling of such that there are now plausible signal values inserted between each adjacent pair of elements of . 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 and , the over- and undersampling parameters, such that . If you want a sampling rate that is ten times the input-signal rate, pick and . But you could also pick and .

In this usage message, there is one warning. It concerns the edges of the input . If the signal starts, relative to zero, too large or ends, relative to zero, too large, the endpoints of the output might exhibit larger error than interior points in . 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 or

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 . The recipe is to Fourier transform to obtain for , then prepend and append zeros to to form , and finally to inverse Fourier transform to obtain the resampled signal .

Let’s choose , where the integer . A typical choice might be or . Then the vector is defined by

so that the vector has total length ; it is times as long as and .

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

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.

So if we pretend that our signal 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 . The symbol rates are both and the excess bandwidths are , so that each signal has an occupied bandwidth of . The number of samples is . 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 and . 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!

In the next example, all variables are the same as in the first except the carrier frequencies are now instead of . 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 is twice that of the signal at , you can see that the energy near the resampled BPSK signal at must be due to the signal at and vice versa.

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