10
$\begingroup$

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!

$\endgroup$
3
  • $\begingroup$ Well. How about you multiply your delay by 1024. It seems to me that your delay works well.when the delay is close an integer number of samples $\endgroup$
    – Ben
    Sep 3, 2019 at 0:52
  • $\begingroup$ @Ben This doesn't work because you can still end up with a fractional delay. $\endgroup$
    – Axemaster
    Sep 3, 2019 at 20:34
  • $\begingroup$ No... I was not suggesting a solution I was pointing out that your solution worked when the delay was close to an integer number of samples and didn't work when it wasn't. $\endgroup$
    – Ben
    Sep 3, 2019 at 20:46

6 Answers 6

4
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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. Fractional Delay Impulse Response with a Tail Artifact

$\endgroup$
0
$\begingroup$

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

$\endgroup$
0
$\begingroup$

Stumbled across this and thought I'd share the logic behind the way to do this for non-integer delays.

To time delay a sequence you convolve it with an "impulse" at the desired delay.

However, since it's discrete time the bandwidth of the system is inherently limited to pi. Therefore this "impulse" is actually a sinc function. Therefore h[n] = sinc(nT - t0). If the delay is an integer number of samples the sequence is h[n]= 0, 0,... 1, 0, 0,... If the delay is a fractional number of samples then the sequence has lots of non-zero elements.

Note that in the frequency domain this sinc function is simply H(f) = e^(k f).

A sinc function is technically infinite therefore h[n] will exist for n < 0 although only 100 or so samples are actually significant.

If we want to do the convolution in the frequency domain the convolution is circular. Therefore the samples of h[n] for n < 0 need to be wrapped to N-1.

So we could generate a sinc sequence at the desired time delay (by sampling the sinc function), wrap the elements for n < 0 and then take the FFT. This is kind of a pain because we don't typically generate sequences for negative indexes. A somewhat easier way to do this is to generate the sinc function offset by N/2 and then circularly shift the sequence back by N/2. The circular shift back is equivalent to swapping the lower and upper halves of the sequence.

Yet another way to do it is to generate the sinc function offset by N/2, take the FFT and then shift the sinc function back in the frequency domain. Due to the duality of the time shift and frequency shift properties of the Fourier Transform this is equivalent to swapping the lower and upper halves of the FFT.

Or we can just do it all in the frequency domain. To generate the sinc offset by N/2 we let H(f) = e^(-2 pi (f-1/2) delta) where delta is the delay in samples and f is the normalized frequency.

Then swap the lower and upper halves to effectively shift the time sequence back by N/2. In Matlab you can use fftshift() to do this (or ifftshift() as they are equivalent in this case).

The algorithm is thus:

  1. Take the input sequence, x[n], and zero-pad it (typically to a power of 2). I.e., if x[n] has 512 samples choose N = 1024 and zero-pad x[n] to 1024 samples.

  2. Generate the sinc sequence, h[n], in the frequency domain using H(f) = e^(-2 pi (f-1/2) delta) and then swapping the left and right halves.

  3. Take the FFT of x[n] -> X(f).

  4. Y(f) = X(f) * H(f).

  5. Then y[n] = real(IFFT(Y(f))).

The actual algorithm is almost trivial. The logic behind it is somewhat complex.

$\endgroup$
0
$\begingroup$

An alternative method:

  1. Upsample by a factor 3. Use pad & filter, or Polyphase, or FFT O&A / O&S - this is more of an implementation issue.
  2. Use a cubic spline to calculated the needed interpolated values

The spline interpolation does not perform well for frequencies that at close to Nyquist, thus the reason for the initial upsampling - it provides more room between the highest frequency and the sampling frequency.

You can achieve better performance if required by up sampling by a larger factor or using a filter with tighter requirements (stopband attenuation, passband ripple, transition width etc).

This is also useful when you need multiple versions of the signal with different time delays e.g. Broad band beamforming.

$\endgroup$
-1
$\begingroup$

You are performing fft's on a REAL signal (all imaginary values are zero). Replace fft and ifft with rfft and irfft and your original code will work.

$\endgroup$
1
  • 2
    $\begingroup$ Hello Jim and welcome to Signal Processing SE. I am not sure this is the correct answer. FFT and IFFT routines (should) give the same results as their respective real counterparts with the only difference that the latter are faster (they compute half the spectrum and then just mirror it). Have you tried your suggestion and reached correct results just by changing the said functions? If so, would you care to share that with the OP? $\endgroup$
    – ZaellixA
    Jan 29, 2023 at 12:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.