Hi Jason,
I think the fundamental issue is that one shouldn't expect that Result = 0.
Suppose we have a continuous time signal s(t).
Then we form
s1(t) = s(t + t_d)
s2(t) = s(t - t_d)
Taking Fourier transforms we have
S1(w) = exp(1j*w*t_d)S(w) (1)
S2(w) = exp(-1j*w*t_d)S(w)
So, we could multiply S1(w) and S2(w) by the inverse exponentials to recover S(w), and therefore s(t), from either.
But, when working with the DFT/iDFT (implemented by fft and ifft) we have to apply a window before sampling. In this case, the window, w(t), is rectangular from -1 <= t <= 1 (sorry for the abuse of notation between frequency (w) and window function (w(t)) ). So we really have s1p(t) = s(t + t_d)w(t) (2)
and similarly for s2(t).
Therefore, the Fourier transform of s1p(t) does not have the form of (1), so we can't expect to recover S(w) by multiplying S1p(w) by the complex exponential.
In fact, we see that if we take the Fourer transform of (2) and then multipily the complex exponential, we get
S1r(w) = exp(-1j*w*t_d)*S1p(w)
with inverse transform
s1r(t) = s1p(t - t_d) = s(t)w(t - t_d)
which is our base signal multiplied by a shifted window (and similarly for s2r(t)).
So we shouldn't expect the s1r(t) - s2r(t) to be zero. It will really be
D(t) = s1r(t) - s2r(t) = s(t) ( w(t - time_delay) - w(t + time_delay) )
And that's before we get to issues with sample rates (as identified by @David Goodmanson in his Answer) and circular convolution, as discussed below. Let's look at the problem in continuous time
I'll leave this line commented and come back to it at the end.
Define the base signal
theta = 2*sym(pi)/lambda*V^2/R*t^2;
Define a numerical function for later use
sfunc = matlabFunction(s);
Define the shifted functions
s1(t) = s(t + time_delay);
s2(t) = s(t - time_delay);
Define the window
w(t) = rectangularPulse(-a,a,t);
Continuous-time Fourier transforms,
syms S1p(omega) S2p(omega)
S1p(omega) = fourier(s1(t)*w(t),t,omega);
S2p(omega) = fourier(s2(t)*w(t),t,omega);
Multiply by the complex exponentials to compensate for the time delay AS IF the s1(t) and s2(t) are not windowed.
S1r(omega) = exp(-1j*omega*time_delay)*S1p(omega);
S2r(omega) = exp(+1j*omega*time_delay)*S2p(omega);
Go back to time domain
s1r(t) = ifourier(S1r(omega),omega,t);
s2r(t) = ifourier(S2r(omega),omega,t);
We actually didn't need to go through all that because we already know what s1r(t) (and s2r(t)) should be, which we can verify, but maybe it was instructive.
simplify(s1r(t) - s(t)*w(t - time_delay))
Compute the difference
D(t) = simplify(s1r - s2r)
We see that D(t) is composed of a complex exponential multiplied by a heaviside sum that's strictly real. So we can plot the magnitude of D by plotting the magnitude of the second term over the time interval of interest
fplot(abs(Dc{2}),[-1,1],'-o'),grid
D(t) is zero, EXCEPT at the endpoints, which will be an important point below.
In order to use discrete-time processing, we need to define the sampling frequency (inspired by @David Goodmanson) First we get the instantaneous frequency as a function of time
frequency(t) = diff(theta,t);
By the nature of the signal, we know the maximum frequency occurs at the end of the window.
maxfrequency = frequency(1);
maxfrequency = maxfrequency/2/sym(pi);
Define an initial guess for the sampling frequency as 20x the max frequency in the signal. Note that we also need the sampling frequency to be large enough to capture the effect of the delay.
Fa = double(20*maxfrequency);
Get an integer number of points that will cover the length of the window (N is odd, which is nice for symmetry)
Define the sampled time vector over the window of interest
Get the actual sampling period and frequency.
Define the corresponding frequency vector
nvec = (-(N-1)/2:(N-1)/2);
Note that time_delay is not an integer multiple of the sampling period, but at least Fa is large enough so that the delay is roughly 7 samples and so will be apparent on the plots below.
Samples of the delayed and windowed signals
s1 = sfunc(ta + time_delay);
s2 = sfunc(ta - time_delay);
In discrete-time, mulitiplication in the frequency domain by exp(-1j*2*pi*f*n0) corresponds to a circular shift of n0 samples in the time domain. Recall from above that fthe signal that we are recovering from this process, D, is very significant at both ends. Hence, the circular shift would cause the points on the left to couple over to the right, and the points on the right to couple over to the left after muliplication by exp(-1j*2*pi*f*n0) in the frequency domain. The way to mitigate this is to zero-pad the signals on both ends in the time domain, then process, then extract the central portion. This way, the circular shift is just shifting zeros around. Here I'll zero pad with 1e3 samples, though that might be excessive for this problem
nvecpad = [-fliplr(npad) + nvec(1),nvec, nvec(end) + npad];
The corresponding frequency vector is
fapad = nvecpad*(Fa/numel(nvecpad));
And the zero-padded signals are
Now we do the FFT stuff. I believe the original code was missing a couple of shift operations, which I've added below.
S1 = fftshift(fft(ifftshift(s1)));
S1 = S1.*exp(-1j*2*pi.*fapad.*time_delay);
s1r = fftshift(ifft(ifftshift(S1)));
Extract the part we care about
s1r = s1r(Npad+1:end-Npad);
S2 = fftshift(fft(ifftshift(s2)));
S2 = S2.*exp(1j*2*pi.*fapad.*time_delay);
s2r = fftshift(ifft(ifftshift(S2)));
s2r = s2r(Npad+1:end-Npad);
The resulting difference looks a lot like the continuous-time plot above based on the symbolic stuff.
plot(ta,abs(Result),'-o'); xlabel("t(s)");
title('After Substraction');
Zoom in at both endpoints and compare magnitude and phase of the discrete-time result to the continuous-time result
plot(h1,ta,abs(Result),'-o')
xlim(h1,[-1.0002,-0.9998])
plot(h1,ta,abs(Df(ta)),'-x')
plot(h2,ta,abs(Result),'-o')
plot(h2,ta,abs(Df(ta)),'-x')
We see that the continuous-time signal amplitude is a step function at both ends, as expected from the heaviside sum in D. The discrete-time processing kind of recovers the same amplitude.
plot(h1,ta,angle(Result),'-o')
xlim(h1,[-1.0002,-0.9998])
plot(h1,ta,angle(Df(ta)),'-x')
plot(h2,ta,angle(Result),'-o')
plot(h2,ta,angle(Df(ta)),'-x')
The phase matches well in the region where the amplitude is large.
Why don't the amplitude and phase match better? I'm pretty sure it's because the time_delay is not an integer multiple of Ts (see Ndelay above) and therefore the multiplication by the complex exponentials in the discrete frequency domain isn't really doing what we we'd like it to do. To see this, you can uncomment the commented time_delay line at the top of the code, which will force time_delay to be exactly 7 samples, which will make the above plots look much better. There will still be some error, which I think is just due to numerical innaccuracy.
Summary: The fundamental approach for recovering s(t) isn't quite right, and if we still want to proceed down the discrete-time processing path, we need choose a high-enough sampling frequency and understand the effect of multiplying by exp(-1j*w*N0) in the discrete-frequency domain.
