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