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

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

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!

Hello Chad,

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

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

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

Nothing changes provided that 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 . Then you resample the sequence with where 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 . They will be (considering the non-conjugate cycle frequencies here) multiplied by . Suppose one of the original cycle frequencies was 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 changed to .

For that first case, where we find the cycle frequency of 0.5, suppose the physical sampling rate is . Then that cycle frequency in physical Hz is . In the second case, we’ve resampled that signal so that and we find a normalized (unit sampling rate) cycle frequency of , which is a physical cycle frequency of . 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 and you resample the data (again, being careful not to distort it), you still have a signal with symbol rate even though the number of samples per symbol has changed.