0

I have a 1d signal obtained using a Fourier based resample method (TDIFDZP) for which the resampled points don't necessarily go through the original samples. I want to transform the upsampled signal to go through these points to reduce the error between the resampled and the original signal as shown in the image below. I'm using scipy.sigal.resample implementation in python.

I'm looking for a transformation of the resampled signal so that it goes through the original samples.

enter image description here

Edit: A little of background, I'm exploring different interpolation methods for sub-Nyquist sampling. My original signal's Nyquist rate is ~20Hz but the number of samples available in my application is <6.

Massimo Panda
  • 92
  • 8
  • that sounds like you've picked the wrong resampling method, then. Why did you chose it, then? – Marcus Müller Jun 26 '20 at 14:16
  • I'm comparing different sub-Nyquist sampling methods for an application where the number of available samples is very low (<6). The Nyquist rate for this signal is ~20Hz. I was trying to think outside-the-box when I came up with the idea of transforming the resampled signal so that I would "artificially" introduce higher frequency components at the same time that I eliminate the Gibbs phenomenom. – Massimo Panda Jun 26 '20 at 14:47

1 Answers1

2

To have all the samples in the interpolated waveform land on the original samples implies using an integer interpolation rate (otherwise time domain distortion would need to be introduced, and I can't think of an application where that would make sense to do).

For integer interpolation, a very simple approach is to use the DFT with zero padding since the samples chosen will be exact, but this will have more deviation in between then other optimized interpolation methods further detailed below. This approach is very simple, so if the deviation between samples is of no concern then this approach would be favored and would be as follows:

Compute the FFT of the samples without zero padding.

Scale the result by $N$, where $N$ is the total number of original samples, zero pad the result and then compute the inverse FFT. (Due to reciprocity the fft and ifft operations can be swapped).

In MATLAB or Octave this is:

N = length(test);
interp = 5;
result = ifft(N*fft(sequence), interp*N);

Other optimized interpolation methods using zero-insert and subsequent image filtering would also preserve the values, and with post processing the image filtering can be done with a zero-phase filter using filtfilt() and for an optimized solution in the least squares sense, together with multiband filter coefficients using firls(). The zero-insert and filtering approach is further detailed here: What is the impulse response used in an interpolation filter when upsampling?

If your samples are non-uniform to begin with, first resample to a uniform rate, where this post already addresses that approach: What is an algorithm to re-sample from a variable rate to a fixed rate? But if the new samples do not fall on the same time locations (which will occur if there isn't an integer relationship in the duration between the original samples) then the interpolated samples should not be expected to match the original samples when they occur at new time samples that were not in the original sequence.

Dan Boschen
  • 50,942
  • 2
  • 57
  • 135
  • What requirements for FFT size and interpolation factor need to be met in order for the reconstructed signal to go through the original samples? I understand the approach you described is the one implemented in scipy.signal.resample in Python. Thank you for the answer. I'll look at the alternative methods proposed. – Massimo Panda Jun 26 '20 at 14:50
  • @Pedro No requirement. How many samples do you want to interpolate by? It works in all cases, try it out! You can do this directly in Python using fftpack identical to what I did and confirm for yourself that it is not exactly what is implemented, or rather the way you used it. – Dan Boschen Jun 26 '20 at 14:52
  • I will check and get back to you. It was not working for me using scipy's interpolation. My signal length is 4, 5, or 6 and my reconstructed signal length is 50 or more. – Massimo Panda Jun 26 '20 at 14:54
  • "test" refers to the original sequence whatever it is you may have. "interp" is the integer interpolation rate you desire, but assumes you are interpolating by an integer. If not an integer then the underlying waveform itself wouldn't go through all the samples so wouldn't really make sense to desire that outcome. – Dan Boschen Jun 26 '20 at 14:55
  • Fundamentally to do what you want, you have to be looking for an integer interpolation rate. Doesn't make sense otherwise as you would then imply that you are introducing time distortion to the signal. Just want to make sure that is clear to you. So 5 going to 50 samples, or 6 going to 60 samples for example would make sense to do. Of course you could interpolate the 6 samples to 60 and then truncate the result to less samples, but you can't interpolate by 5.5 and expect all samples to land on the originals. – Dan Boschen Jun 26 '20 at 14:56
  • By time-domain distortion, do you refer to transforming the resampled signal in the time domain (after interpolating)? – Massimo Panda Jun 26 '20 at 15:03
  • Let us continue this discussion in chat. – Dan Boschen Jun 26 '20 at 15:04
  • Thank you again Dan. I understand it's strictly required for the sampling time to be uniform. Because of the way I choose the samples based on an integer index there are some small differences for some samples. – Massimo Panda Jun 26 '20 at 15:34