How to do FFT fractional time delay (SOLVED)

Question

I am trying to time-shift a signal using the FFT, however I have encountered some strange effects that depend on the size of the time-shift. I need to be able to shift the time by an arbitrary amount - i.e. floating-point time. I am using the method from this thread:

Here is my code:

import numpy as np
import matplotlib.pyplot as plt

f1 = 1.8
f2 = 2.6
#try tDelay = .02002 and tDelay = .0205
tDelay = .0205 #seconds
samples = 1024 #number of samples in the time interval
tstart = 0.0
tend = 1.0

# create a waveform to use for the time shifting
samplePeriod = (tend - tstart) / (samples)
print("\nThe sampling period is %f seconds" % samplePeriod)
print("The time delay is %f seconds" % tDelay)
tDelayInSamples = tDelay / samplePeriod
print("The time delay in samples is %f samples" % tDelayInSamples)
timeList = np.linspace(tstart, tend, samples)
waveform = np.sin(2 * np.pi * f1 * timeList) + np.sin(2 * np.pi * f2 * timeList)

# do the time shifting
fftOut = np.fft.fft(waveform)
N = fftOut.shape[0]
k = np.linspace(0, N-1, N)
phaseShiftFunction = np.exp((-2*np.pi*1j*k*tDelayInSamples)/(N))
fftWithDelay = np.multiply(fftOut, phaseShiftFunction)
waveform2 = np.fft.ifft(fftWithDelay)

plots = 1
plt.subplot(plots, 1, 1)
plt.plot(waveform)
plt.plot(waveform2)
plt.show()

If you run the above code using tDelay like .02002 and .0205, you will see that the recovered signal is drastically different. I don't understand why this is the case, or what to do about it.

Example of good time-shift (tDelay = .0205):

Example of bad time-shift (tDelay = .02007):

EDIT:

After considerable effort I believe I have solved this problem. There are 2 key insights here:

1. When you do a fractional time shift, you have to do an fftshift on the array of phase shift values, to make it symmetric. (This is not needed for integer-sample time shifts). See this thread.

2. After multiplying the FFTed data with the phase shift values, if you immediately do the IFFT, you will find that the data is correctly time-shifted, but the waveform is rotated at some angle in the complex plane. This angle depends on the sample fraction corresponding to the time shift. For example, a time shift corresponding to N+.5 samples will be entirely imaginary (i.e. rotated by $\pi/2$). To rotate back to the real axis, add an additional phase shift (before doing IFFT): $e^{\pi i D}$ (where D is the time-shift in samples).

So just to be super clear, the process is:

1. Take the FFT
2. Construct the phase shift $e^{\frac{-2\pi i k D}{N} + \pi i D}$ (where k=sample number, D=time shift in samples, N=sample length of the FFT)
3. Do FFTshift on the phase shift coefficients calculated in #2
4. Now multiply the FFT data with the phase shift coefficients
5. Do the IFFT

Example Code:

import numpy as np
import matplotlib.pyplot as plt

f1 = 12.8
f2 = 22.6
samples = 1024
tDelay = .00938
tstart = 0.0
tend = 1.0

# 0. Example waveform to demonstrate the time shift
timeList = np.linspace(tstart, tend, samples)
waveform = np.sin(2 * np.pi * f1 * timeList) + 1*np.sin(2 * np.pi * f2 * timeList)

# 1. Take the FFT
fftData = np.fft.fft(waveform)

# 2. Construct the phase shift
samplePeriod = (tend - tstart) / (samples)
tDelayInSamples = tDelay / samplePeriod
N = fftData.shape[0]
k = np.linspace(0, N-1, N)
timeDelayPhaseShift = np.exp(((-2*np.pi*1j*k*tDelayInSamples)/(N)) + (tDelayInSamples*np.pi*1j))

# 3. Do the fftshift on the phase shift coefficients
timeDelayPhaseShift = np.fft.fftshift(timeDelayPhaseShift)

# 4. Multiply the fft data with the coefficients to apply the time shift
fftWithDelay = np.multiply(fftData, timeDelayPhaseShift)

# 5. Do the IFFT
shiftedWaveform = np.fft.ifft(fftWithDelay)

print("\nThe sampling period is %f seconds" % samplePeriod)
print("The time delay is %f seconds" % tDelay)
print("The time delay in samples is %f samples" % tDelayInSamples)
print("The correction phase shift is %f pi" % (tDelayInSamples))

plots = 1
plt.subplot(plots, 1, 1)
plt.plot(waveform)
plt.plot(shiftedWaveform)
plt.show()

I welcome any additional insight someone might have for better ways to do this, or if I have overlooked anything!

Answer 1

In either Overlap-add or Overlap-save, the FFT is doing the Discrete-Fourier Transform that periodically extends your input data. The DFT does only circular convolution, so you need to make this tool that does circular convolution into a tool that does linear convolution. And Overlap-add and Overlap-save are an adaptation of a circular convolver to do the task of linear convolution. And because the computational savings of the FFT, it turns out to be worth it for very long FIR filters.

So what I meant to say is that; if your purpose is to delay a signal by a delay amount that is expressed in precision of fraction of a sample, then you need to convolve with an impulse response representing the delay and use Overlap-add or Overlap-save to do it. that's the only glitch-free way to implement a fractional delay over a streaming or very long input in blocks.

Answer 2

It’s because .0205 is reeeeaaally close to 30 samples, but .02002 is reeeeaaallly close to 20.5 samples. This scheme only works for integer sample delays, and is not well suited for arbitrary delays. Unfortunately, I haven’t ever tried to do fractional delays, so I don’t know that I can give you any advice on how to do it any better. However, that is why you are seeing what you are seeing.

Edit I'm saving my previous post for posterity, but I see now what you were doing. You are performing a circular shift, in that the end of the waveform wraps around to the beginning. If you were to try this using a real time/causal filter, you might get some nasty results. I used your code to generate the corresponding impulse response. As you can see, the tail does not 'fade' to zero, indicating that you would get artifacts as a result.

Answer 3

The phase shift formula you originally used does not guarantee conjugate symmetry property of the fractionally delayed signal (only works fine for integer delays), as you might have already noticed. A good solution for the problem is described here