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

Sometimes MATLAB’s resample.m gives results that can be trouble for subsequent CSP.

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 = 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 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 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 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.

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

1. Fabio Casagrande Hirono says:

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 . 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!

 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

1. Chad Spooner says:

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!

2. Anat says:

Thank you for your blog! I’m new to CSP and am learning a lot, both from your posts and the comments section.

My question doesn’t relate to the resampling problem you’ve raised, rather, I’d like to ask –
What would happen to the cyclic autocorrelation function if the resampling ratio P/Q > 1 and is not an integer?
How will the cycle frequencies be effected by this ratio?
I didn’t manage to write the formula for CAF for this case.

Thanks

1. Chad Spooner says:

Welcome to the CSP Blog Anat! And thanks for the comment.

Thank you for your blog! I’m new to CSP and am learning a lot, both from your posts and the comments section.

You are welcome, and I too find the Comments on the CSP Blog helpful.

What would happen to the cyclic autocorrelation function if the resampling ratio P/Q > 1 and is not an integer?
How will the cycle frequencies be effected by this ratio?

Nothing changes provided that $R = P/Q$ isn’t so small that you are distorting the signal. That is, if you continue to obey the sampling theorem, the cyclic autocorrelation doesn’t change. This assumes that you are considering the physical sampling rate for the samples you are processing. Things look a little different if you focus on normalized cycle frequencies.

Suppose you have some samples, and the signal in those samples is comfortably oversampled. You analyze the signal, and find cycle frequencies in some set $A$. Then you resample the sequence with $R = P/Q$ where $R$ is not an integer. Then you view the resulting sequence as corresponding to a sampling rate of one. The cycle frequencies you find for that sequence will not be those in $A$. They will be (considering the non-conjugate cycle frequencies here) multiplied by $R$. Suppose one of the original cycle frequencies was $0.5$ and $R = 10/11$, so that the sampling rate is slightly reduced. Then, in the new sequence, pretending that the sampling rate is one, you’ll see that $0.5$ changed to $(11/10)0.5$.

For that first case, where we find the cycle frequency of 0.5, suppose the physical sampling rate is $F_1$. Then that cycle frequency in physical Hz is $0.5F_1$. In the second case, we’ve resampled that signal so that $F_2 = (10/11)F_1$ and we find a normalized (unit sampling rate) cycle frequency of $(11/10)0.5$, which is a physical cycle frequency of $F_2(11/10)0.5 = (10/11)F_1(11/10)0.5 = 0.5F_1$. So the physical cycle frequency is the same, as it needs to be. This is analogous to saying that if you have a sampled digital signal with symbol rate $f_{sym}$ and you resample the data (again, being careful not to distort it), you still have a signal with symbol rate $f_{sym}$ even though the number of samples per symbol has changed.