Recovering Time Shift Using DFT of Translated Square Pulse?

72 Views Asked by At

As an exercise, I attempted to manually translate a pulse $n_0$ steps to the right and recover the translation using the time-shift property. The problem I'm encountering is that the phase unwrapping of the fft of the translated pulse is not a straight line like it should be, despite that the theoretical and manually translated pulses appear to be exactly the same. I've written a very simple script to illustrate this. Can someone explain what I'm doing wrong? This is important because in reality, I won't have $n_0$.

[EDIT]

Here is the math showing how to recover the translation.

$$y(t) = x(t-t_0) = x(t)*h(t)$$ $$Y(\omega) = X(\omega)H(\omega)$$ By the time shift property: $$H(\omega) = e^{-j\omega t_0}$$ $$e^{-j\omega t_0} = Y(\omega)/X(\omega)$$ $$t_0 = -ang(H(\omega))/w$$

from numpy import *
from scipy import *

N = 1000
n0 = 200

pulse =  zeros(N)
pulset = zeros(N)
pulse[100:200] = 1
pulset[100+n0:200+n0] = 1 # t0 = 100

f = arange(0, 1, 1./N)
w = 2*pi*f

X = fft(pulse)
Y1 = fft(pulset)
Y2 = exp(-1j*w*n0)*X # theoretical transformation

plot(Y1) # Y1 and Y2 look the same
plot(Y2)
print(sum([abs(Y1[i] - Y2[i]) for i in range(0, len(Y1))])) # 1.12701544921e-10

figure()

t1 = -unwrap(angle(Y1/X))/w
t2 = -unwrap(angle(Y2/X))/w

plot(t1)
plot(t2)

figure()

plot(unwrap(angle(Y1/X)))
plot(unwrap(angle(Y2/X))) 

[EDIT2] I posted pictures of the super imposed FFTS of Y1 and Y2, unwrapped phase and calculated time shift.

[EDIT3] Same pictures (minus the FFTs) but with multipling by conj(X) instead of dividing by X

[EDIT4]

So after looking at MATLAB code (in comments) I looked into the phase unwrapping and indeed after I changed the discont parameter to 5, the phase unwrapped more correctly meaning it was unwrapping where it should not have (second last picture) but looking at the calculated time shift (last picture) there were still spikes. I then looked at the super imposed wrapped phase (third last picture - I zoomed in) and indeed they were different, unlike in MATLAB where they are exactly the same. It is worth noting that I am using 32 bit python and I ran the MATLAB code in 64 bit MATLAB although I'm not sure if that has anything to do with it.

Super imposed FFTs

Super imposed unwrapped phases

Super imposed calculated time shifts

Super imposed unwrapped phase - multiplied by conj(X) instead of division by X

Super imposed calculated time shift - multiplied by conj(X)

Superimposed wrapped phase of Y1 and Y2. notice they are not the same

Super imposed Unwrapped phase after increasing discont parameter to 5

Super imposed calculated time shift