Phase gradient in fft

I recently found out that zero padding can affect the phase gradient of the Fourier transformed data. It seems that depending on the position of the signal, the phase gradient can be positive or negative as shown in the attached figure. I am not talking about the additional phase gradient due to a shift in time signal. Could someone please explain why that is so?

The following code is used

t = 0:1/40:10;  t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty = fft(y);  fty(abs(fty)<1e-6) = 0;
theta = unwrap(angle(fty));  dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
subplot(2,2,1)
plot(t,y);  title('Time signal') 
subplot(2,2,3)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

t = 0:1/40:20;  t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty = fft(y);  fty(abs(fty)<1e-6) = 0;
theta = unwrap(angle(fty));  dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
subplot(2,2,2)
plot(t,y);  title('Time signal') 
subplot(2,2,4)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

[Paul]

Can you post the code that produced those plots?

[Jiayun Liu]

I will post the code as soon as I have access to them. The title of the bottom plots is supposed to be unwrapped phase. Sorry for that

Answer 1

[Paul]

Hi Jiayun,

First of all, I don't see anything in the code that does "zero-padding." That term typically refers to augmenting a sequence with zeros prior to taking the FFT. But, the code in the Question is just taking more sample points of the underlying continuous-time signal, which isn't zero padding.

I pretty sure all you're seeing is the effect of

a) taking a finite number of samples of an infinitely long signal (i.e., a windowing effect), so we'd expect that the more samples taken the closer the DFT will be to samples of the CTFT, which is what I think you're expecting for the phase plots.

b) the code sets a lot of samples of the DFT to zero (fty(abs(fty)<1e-6) = 0;) which also loses all phase information for those points.

Start with an exact signal and its Fourier transform

syms t
t0 = -7;
x(t) = sin(2*sym(pi)*10*(t+t0))*exp(-(t+t0)^2);
xfunc = matlabFunction(x(t));
fplot(x(t),[0 10])
xlabel('t'); ylabel('x(t)')

syms w f
X(w) = (fourier(x(t),t,w));
X(f) = X(2*sym(pi)*f);
Xfunc = matlabFunction(X(f));

Plot the phase of the CTFT, both wrapped and unwrapped

fval = 0:.01:20;
figure
plot(fval,angle(Xfunc(fval)))
xlabel('f (Hz)'); ylabel('angle(X(f)) (rad)')

figure
plot(fval,unwrap(angle(Xfunc(fval))))
xlabel('f (Hz)'); ylabel('angle(X(f)) (rad)')

That funny business of zero phase at the edges is because X(f) numerically evaluates to zero at points near the ends, e.g.,

Xfunc(fval([1 2 end-1 end-2]))

ans = 1×4

0 0 0 0

Settting that issue aside, we see (exact?) linear phase as function of frequency.

Now, sample the signal from t = 0 to t =10 a sampling frequency of Ts = 1/40 and compute its DTFT with lots of frequency points

dt = 1/40;
df = 1/dt;
tval = 0:dt:10;
xval = xfunc(tval);

[h,w] = freqz(xval,1,8192);

figure
plot(w/pi*20,angle(h))
xlabel('f (Hz)'); ylabel('angle(XDTFT(f)) (rad)')

If you look carefully, it's apparent that the character of the phase is different for frequencies in small window 8 < f < 12, and frequencies outside that window. We can see this more easily after unwrapping, which reveals change in slope in the center of the graph.

figure
plot(w/pi*20,unwrap(angle(h)))
xlabel('f (Hz)'); ylabel('angle(XDTFT(f)) (rad)')
grid
xlim([5 15])

Instead of the DTFT, let's take the N-point DFT, where N is the number of samples.

N = numel(tval);
Xval1 = fft(xval,N);
fval1 = (0:(N-1))/N*df;

Plot the angle of the DFT samples on the phase plot of the DTFT

plot(w/pi*20,angle(h),fval1,angle(dt*Xval1),'o')
xlabel('f (Hz)'); ylabel('angle(XDTFT(f) XDFT(f)) (rad)')
xlim([0 20])

plot(w/pi*20,angle(h),fval1,angle(dt*Xval1),'o')
xlabel('f (Hz)'); ylabel('angle(XDTFT(f) XDFT(f)) (rad)')
xlim([8 12])

As must be the case, the angle of the DFT samples are samples of the the angle of the DTFT, but the effect of the time-domain windowing and sampling on the phase of the DTFT makes the frequency-domain samples of the DFT take on an unexpected shape.

I suggest going through the exact same code as above, but with

t = 0:dt:20;

to see how taking more samples of the continuous-time signal impacts the analysis. It should yield a DTFT that closer to the CTFT, and therefore the DFT samples may be closer to what may be expected.

[Jiayun Liu]

Hi Paul,

Regarding both comments,

a) From what I understand, zero padding in matlab is achieved through concatenating a vector with zeros in it with your original signal. In my code I did it through increasing the vector t. Since the Gaussian will attenuate the signal, Is this not equivalent to adding zeros? Nevertheless, If I just add zeros to the back of the signal through concatenating, I get the same result.

b) fty(abs(fty)<1e-6) = 0; in the code is just to make the phase in frequency domain look 'nicer' by removing those Fourier components with small amplitude. This is mainly used to remove signal caused by spectra leakage which could be due to windowing like you mentioned. Again, even if this line is removed, I observed the same effect.

We evaluate the time signal differently. You used symbolic function to get the exact signal while I use a sine.*Gaussian (kinda like sampling your signal). You can see that phase gradient that we evaluate are different. Mine is positive while yours is negative. I believe that negative phase gradient should be correct since the fft kernel is exp(-iωt). So why don't I get negative gradient?

I can't evaluate the signal using symbolic function since I am dealing with experimental signals with the gaussian peak at different positions for different parameters.

[Paul]

Hi Jiayun,

a) zero-padding means concatenationg zeros onto the end of the samples. This definition isn't just a Matlab convention, it's a standard terminology as far as I know (example reference). Increasing the length of the t vector is actually taking more samples of the signal. In this particular case, those sample are quite small and so the result is approximately the same as zero-padding, but it's not really zero-padding.

b) Yes, it does remove those small amplitude components, as long as we also realize that it sets the phase of those samples to zero as well. Forcing the phase to zero at those samples will elminate the linear phase response through those samples that I think you're looking for.

c) I'm not really using the symbolic function. I converted x(t) to an anonymous function that is functionally equivalent to how you evaluate y at the samples in the t-vector. Consequently, the result that I get in the variable Xval1 is essentially the same as what you get in fty (before comparing to 1e-6). So the phase response I showed above for angle(Xval1) in the red circles is the same as yours, and if you apply unwrap() to angle(Xval1) we'll get the same result.

syms t
t0 = -7;
x(t) = sin(2*sym(pi)*10*(t+t0))*exp(-(t+t0)^2);
xfunc = matlabFunction(x(t));

dt = 1/40;
df = 1/dt;
tval = 0:dt:10;
xval = xfunc(tval);

N = numel(tval);
Xval1 = fft(xval,N);
fval1 = (0:(N-1))/N*df;

Here's the same plot I generated above, but with the red circles connected:

figure
plot(fval1,angle(dt*Xval1),'-ro')
xlabel('f (Hz)'); ylabel('angle(XDFT(f)) (rad)')
xlim([0 20])

Here is the same plot but with the unwrapped phase; it looks the same as yours.

figure
plot(fval1,unwrap(angle(dt*Xval1)),'-ro')
xlabel('f (Hz)'); ylabel('angle(XDFT(f)) (rad)')
xlim([0 20])

All of these results make sense to me. The DFT samples (as computed by fft()) are frequency-domain samples of the DTFT of the the finite number of samples of x(t). Maybe the issue is with how unwrap() is unwrapping the phase to yield a visually unexpected result?

[Jiayun Liu]

I really do not know what the issue is. Yes, we can get the same plot but should the gradient of the phase be positive in this case? The fft kernel is not

syms t
t0 = -7;
x(t) = sin(2*sym(pi)*10*(t+t0))*exp(-(t+t0)^2);
xfunc = matlabFunction(x(t));
dt = 1/40;
df = 1/dt;
tval = 0:dt:10; tval = [tval,zeros(1,801)];
xval = xfunc(tval);
N = numel(tval);
Xval1 = fft(xval,N);
fval1 = (0:(N-1))/N*df;

figure
plot(fval1,unwrap(angle(dt*Xval1)),'-ro')
xlabel('f (Hz)'); ylabel('angle(XDFT(f)) (rad)')
xlim([0 20])

By adding zeros to tval in your code, we see that the phase gradient is negative now. How do we explain this change from positive to negative phase gradient by simply zero padding?

[Paul]

Adding zeros to the end of tval isn't zero-padding. Instead, that actually adds 801 copies of x(0) to the end of the sequence. In this case, x(0) is quite small, so it's almost the same as zero-padding, but that's just for this case.

Let's start with the continuous-time signal x(t). In this problem x(t) has infinite duration to both negative and positive time. The Continuous Time Fourier Transform (CTFT) of x(t) was derived (I called it X(f)) and its phase plotted in the Answer in this thread. Let's assume for the sake of discussion that the phase of X(f) is linear with constant negative slope.

Then, we take samples of x(t) from t = 0 to t = 10. Call that discrete-time signal x[n], which is assumed to be zero for any samples that correspond to t < 0 or t > 10. In so doing, we lose some information about the original signal x(t). Nevertheless, we can compute the Discrete Time Fourier Transform (DTFT) of x[n], which is called h. The DTFT of x[n] is a continuous function of freuqency, but of course we can only compute it for a finite number of points, which I chose 8192 above. The phase of the DTFT very nearly approximates the phase of the CTFT, but as shown above the phase of the DTFT has a slight change in slope for 9 < f < 11 Hz. However, we can still perhaps think of the DTFT as linear phase with a negative slope.

Now, what happens when we use fft() to compute the N-point DFT of x[n], which I called Xval1? Well, by definiton, the samples of Xval1 must be samples of h(f). The specific samples of h(f) are just h(fval1), even though we computed h(fval1) via fft() and named the result Xval1. Now for the value of N used to compute Xval1, the phase of the samples happened to "migrate" across the "bottom" of the phase of curve of h(f) for the lower frequencies and across the "top" of the phase curve of h(f) for higher frequencies. Is there some mystical reason why it came out that way? I do not know. But the end result is that it looks like the phase has positive slope when we "connect the dots" and unwrap, but it only looks like that because of the where the frequency-domain samples lie on the underlying phase curve of the DTFT.

If we increase the density of the frequency-domain samples, then the phase curve looks more like what expect. How can we increase the density? By zero-padding. In the code below, I've zero-padded to increase the density in the frequency domain by a factor of four. In so doing, we see that every fourth zero-padded sample lies exactly on top of every sample without zero-padding. By zero-padding, we fill in samples in the phase between +pi and -pi and the result will look as we expect after unwrapping. So, in the first case the unwrapping was "fooled" because we just didn't have enough points in the frequency domain.

syms t
t0 = -7;
x(t) = sin(2*sym(pi)*10*(t+t0))*exp(-(t+t0)^2);
xfunc = matlabFunction(x(t));
dt = 1/40;
df = 1/dt;
tval = 0:dt:10;  %tval = [tval,zeros(1,801)];
xval = xfunc(tval);

Take the N-point DFT

N = numel(tval);
Xval1 = fft(xval,N);
fval1 = (0:(N-1))/N*df;

Take the 4N-point DFT

Xval2 = fft(xval,4*N);
fval2 = (0:(4*N-1))/(4*N)*df;

Compare the wrapped phase of Xval1 and Xval2. In the first plot, we see that Xval1 marches across the bottom and then the top of Xval2. In the zoomed-in plot, we see that every point in Xval1 coincides with every fourth point in Xval2.

figure
plot(fval2,(angle(dt*Xval2)),'-bx',fval1,(angle(dt*Xval1)),'-ro')
xlim([0 20])

figure
plot(fval2,(angle(dt*Xval2)),'-bx',fval1,(angle(dt*Xval1)),'-ro')
xlim([7 10])

If we increase the mulitplier on N even more, the DFT samples will approach the DTFT.

In short, after windowing and sampling x(t) to form x[n], the transform that applies is the DTFT. The DFT is just frequency-domain samples of the DTFT, so the look and feel of the DFT will depend on where the DFT samples the DTFT in the frequency domain, with those frequency samples being proportional to: (0:(numel(DFT)-1))/numel(DFT).

Perhaps another way to think of this is that the slope of the phase of the DTFT for frequencies not near 10 Hz is roughly 2*pi rad / 0.1 Hz, i.e., the phase phase change of 2*pi occurs over every 0.1 Hz in frequency. But for the first case, we have df = 40 and N = 401, so every DFT sample is spaced by 0.1 Hz, which is what we see with the red circles in the plot above. So the 401-point DFT frequency spacing isn't small enough for unwrap to be able to "see" the true (i.e., DTFT) slope of the phase curve.

[Jiayun Liu]

Yes I agree with what you are saying. So basically you think there is a true phase gradient that we can get if we use continuous Fourier series and the discrete case approaches the continuous case as the sample points increases. What you said could be true but it does not really explain what I observe.

If we look at the question again, I am really saying that the phase gradient changes depending on the location of the signal (either on the left/right half of the whole signal). Zero padding is one way of shifting the signal. Since we always see that zero padding just increases the resolution without adding more information, I thought this way of asking will make more sense. Regardless, we can shift the signal without changing the sample size.

t = 0:1/40:10;  t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty = fft(y);
dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
figure
subplot(2,2,1)
plot(t,y);  title('Time signal') 
subplot(2,2,3)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

t0 = -3;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty = fft(y);
dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
subplot(2,2,2)
plot(t,y);  title('Time signal') 
subplot(2,2,4)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

As is apparant, the phase gradient changes by simply shifting the time signal. So the phase gradient changes due to DFT approaching CFT can't be used to explain this.

t = 0:1/40:10; 

t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
y = [y zeros(1,800)];
dt = t(2)-t(1); t = 0:1/40:30;
fty = fft(y);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
figure
subplot(2,2,1)
plot(t,y);  title('Time signal') 
subplot(2,2,3)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

t = 0:1/40:10; 
t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
y = [zeros(1,400) y zeros(1,400)];
dt = t(2)-t(1); t = 0:1/40:30;
fty = fft(y);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
subplot(2,2,2)
plot(t,y);  title('Time signal') 
subplot(2,2,4)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

Furthermore we can see that this really doesn't depend on sample size as the y values in both graphs above as the same size. The only difference is their position.

Apologies for using zero padding to ask this question. I should have asked it this way from the beginning. I really do not know what I was thinking...

[Paul]

The first example in the comment above is actually a little tricky, because we have to deal with the effect of shifting in the time domain followed by windowing and sampling. However, the second example is much easier to explain because it's a simple time shift, so let's focus on that one.

Here is the first signal

t = 0:1/40:10; 
t0 = -7;
y1 = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
y1 = [y1 zeros(1,800)];
dt = t(2)-t(1); t = 0:1/40:30;
fty1 = fft(y1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta1 = unwrap(angle(fty1));
%figure
%subplot(2,2,1)
%plot(t,y);  title('Time signal') 
%subplot(2,2,3)
figure
plot(ftx(1:floor(end/2)),theta1(1:floor(end/2)));  title('Phase of Fourier components')

And the second signal

t = 0:1/40:10; 
t0 = -7;
y2 = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
y2 = [zeros(1,400) y2 zeros(1,400)];
dt = t(2)-t(1); t = 0:1/40:30;
fty2 = fft(y2);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta2 = unwrap(angle(fty2));
%subplot(2,2,2)
%plot(t,y);  title('Time signal') 
%subplot(2,2,4)
figure
plot(ftx(1:floor(end/2)),theta2(1:floor(end/2)));  title('Phase of Fourier components')

Because of the way that y2 is constructed, the DFT of y2 can also be computed from that of y1 by

fty1shift = fty1.*exp(-1j*2*pi*ftx*dt*400);  % time shift property
theta1shift = unwrap(angle(fty1shift));
figure
plot(ftx(1:floor(end/2)),theta1shift(1:floor(end/2)));  title('Phase of Fourier components')

We see that the unwrapped phase of y2 in the frequency domain has a positive slope because of the effect of the shift of y2 relative to y1 in the time domain, and however unwrap works. We can verify that theta1shift is the same as theta2

figure
plot(ftx,theta2 - theta1shift)

If these results are still troubling, I guess you'd have to look the wrapped phase in all cases and then see what unwrap is doing.

[Jiayun Liu]

I don't get what you are try to say. By multiplying the frequency domain signal with the exponent is the same as shifting the signal in time domain and doing fft. So we would expect the transformed signal to be exactly the same.

I don't get how this is related to the sign of the phase gradient depending on the position of the time signal. Put it another way, the first signal has positive phase gradient, the second has negative phase gradient. By shifting the second signal to the position of the first signal, I get the same phase gradient. But why do they have different phase gradient in the first place?

[Paul]

I think the answer you're looking for lies in the effect sampling in the frequency domain, which I tried to explain before.

Let's consider an 11-point signal with DFT with linear phase

n = 0:10;
w = n/numel(n)*2*pi;  % associated discrete frequency vector, rad/sample
X = exp(-1j*5*n);     % DFT  of the signal, linear phase
x = ifft(X);          % the signal itself, which isn't really important

Now, if we ignore for the moment that n only takes on integer values, it would appear that the phase of X could be computed as

figure
plot((0:.01:10)/11*2*pi,angle(exp(-1j*5*(0:.01:10))))

And if we unwrap, it surely looks like a negative slope, which is what we might expect based on exp(-1j*5*n).

figure
plot((0:.01:10)/11*2*pi,unwrap(angle(exp(-1j*5*(0:.01:10)))))

But, n only takes on integer values, i.e., the DFT only takes on samples of the blue curve.

figure
plot((0:.01:10)/11*2*pi,angle(exp(-1j*5*(0:.01:10))),(0:10)/11*2*pi,angle(X),'-o')

But looking at just the red curve, it looks like it has a positive slope, becasue the frequency domain samples aren't spaced close enough to capture the rapidly, negatively changing phase. As a result, the DFT samples actually look like a positive slope after unwrapping

figure
plot((0:.01:10)/11*2*pi,unwrap(angle(exp(-1j*5*(0:.01:10)))),(0:10)/11*2*pi,unwrap(angle(X)),'-o')

Now, if we zero-pad, we get more frequency samples and the get the expected slope

figure
plot((0:.01:10)/11*2*pi,angle(exp(-1j*5*(0:.01:10))),(0:100)/101*2*pi,angle(fft(x,101)),'-o')

So, the slope of the phase of the DFT samples can be deceiving if the frequency domain sampling isn't small enough.

Now suppose we apply a circular shift to the signal itself and take its DFT (note that your example above was a circular shift)

x1 = circshift(x,5);
X1 = fft(x1);

The DFT X1 is just the DFT of X multiplied by what looks like another linear phase term with negative slope

norm(X1 - X.*exp(-1j*w*5))

ans = 5.6574e-15

What is the effect of that linear phase term on the phase of the DFT of X1?

figure
plot(w,(angle(X)),'-x',w,(angle(X1)),'-o')
grid

Here we see the red samples acutally look like they have a negative slope, and they do after unwrapping

figure
plot(w,unwrap(angle(X)),'-x',w,unwrap(angle(X1)),'-o')
grid

So the time shift can make the phase look like an apparent change in the sign of the slope of the phase plot. But again, that's just a result of the spacing of the frequency domain samples relative to the phase of the DTFT (not the CTFT).

[Jiayun Liu]

If I understand you correctly, you are saying this change in phase gradient is due to insufficient sampling points. This would mean that increasing the sampling points through zero padding would do the job. If that is the case, I would like to ask do we always have to zero pad behind the signal or we can zero pad infront of the signal i.e. [y,zeros(1,300)] vs [zeros(1,300),y]. As I have shown previously, depending on where we add the zero pad, the phase gradient changes while having the same number of sampling points.

Additionally, this doesn't explain the tricky situation that you mentioned previously. In that example, I simply shifted the signal and we get different phase gradient while having the same number of sampling points. Also, why is this case tricky? Besides the shift, windowing and sampling is always present when we discretize signal for fft.

[Paul]

What I'm saying is this: If we have a finite duration sequence, x[n], it has a DTFT, call it X(w). The independent variable of the DTFT, w, is a continuous variable, so it's reasonable to define the "gradient" of the phase, except at points of discontinuity where the phase, by convention, wraps around from -pi to +pi (or the other way). The DFT of x[n] exists only at discrete frequency points, and each element of the DFT is a sample of the DTFT at those discrete frequency points (here, I'm only considering that portion of the DFT as computed by fft()). Because the DFT is only defined at discrete points, we can visually "connect the dots" to come up with what looks like a function of a continuous independent variable. However, if the samples aren't close enough to each other in the frequency domain, then we are not gettting the full picture of the underlying DTFT.

Zero-padding "behind the signal" allows us to fill in additional frequency domain samples to get closer and closer to the underlying DTFT. However, zero-padding "in front of the signal" will result in an addittional, linear, discrete phase shift of the DFT relative to the underlying DTFT. That is if we have sequence y, we can zero pad "behind"

Y1 = fft([y zeros(1,300)]);  % frequency domain samples of Y(w)

or in front

Y2 = fft([zeros(1,300) y]);  % not frequency domain samples of Y(w)

in which case Y2 and Y1 will be related by

Y2 = Y1.*exp(-1j*w*300)

where w is the discrete frequency vector of the DFT in units of rad/sample.

I've never seen any reference that talks about zero-padding in front as an analysis technique, but that's just me.

As for the "tricky situation," I used that term because in that example, the continuous-time signal was shifted and the CTFTs of the shifted and unshifted signals are reated by a linear phase factor of exp(-1j*w*T), where here w is the continous-time frequency variable in rad/sec and T is the time shift, which I think think was -4 seconds. However, for both signals the window was same from 0 < t < 10 seconds. As a result, the discrete-time signals are not related by a simple shift. Consequently, their DTFTs aren't related by a linear phase factor at all frequencies. However, if I recall correctly, the DTFTs were essentially related by the linear phase factor in the frequency window 9 < f < 11 where the amplitude of the DTFTs is large, and in that region the effect of that linear phase factor made it look like the "gradient" of the discrete samples of the phase became reversed in that frequency region, very similar to the example I showed in this comment.

[Jiayun Liu]

Yes I totally agree that zero padding infront of the signal is different from the original in the sense that it adds a term to the transform result. Therefore, zero padding 'behind' the signal should give the same result since I am not shifting the signal. Which isn't the case as we can see below.

clear
t = 0:1/40:10;  t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty = fft(y);
dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
figure
subplot(2,2,1)
plot(t,y);  title('Time signal') 
subplot(2,2,3)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

t0 = -7;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);   y = [y,zeros(1,800)];
fty = fft(y);   t = 0:1/40:30;
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
subplot(2,2,2)
plot(t,y);  title('Time signal') 
subplot(2,2,4)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

When the signal is shifted from the right half to the left half, the phase gradient changes. If you are saying this is because of the insufficient sampling points, we can see that if the signal is kept on the left half, the phase gradient doesn't change.

t = 0:1/40:10;  t0 = -3;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty = fft(y);
dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
figure
subplot(2,2,1)
plot(t,y);  title('Time signal') 
subplot(2,2,3)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

t0 = -3;
y = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);   y = [y,zeros(1,800)];
fty = fft(y);   t = 0:1/40:30;
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta = unwrap(angle(fty));
subplot(2,2,2)
plot(t,y);  title('Time signal') 
subplot(2,2,4)
plot(ftx(1:floor(end/2)),theta(1:floor(end/2)));  title('Phase of Fourier components')

All in all, what I observed is the following:

Suppose we have a time signal and we divide it into left half (1:end/2) and right half (end/2+1:end). If the signal is in the left half, we get negative phase gradient. If the signal is in the right half, we get positive phase gradient. This observation is true regardless of sample size (zero padding doesn't affect the gradient). Is there an explanation to this observation in your comment that I missed?

"However, for both signals the window was same from 0 < t < 10 seconds. As a result, the discrete-time signals are not related by a simple shift." You mean a shift in the time signal with the same time window doesn't correspond to a shift in time where the fft of the shifted and unshifted signal can be related through a term? Why is this the case?

[Paul]

As to the first part of your comment, I don't think I can say any more than I already have. The unwrapped phase of the DFT samples of the signal(s) being discussed here depends on where those samples end up on the wrapped phase of the underlying DTFT, as has been demonstrated above. I don't know how to explain the phase of the DFT samples in any fundamental way regarding "right half" or "left half.". So I can't be of any more help in this regard.

As to the second part of your comment about the time shift and windowing ... Let's go back to that case. Taking the code from this comment

t = 0:1/40:10;  t0 = -7;
y1 = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty1 = fft(y1);
dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta1 = unwrap(angle(fty1));
figure
subplot(2,2,1)
plot(t,y1);  title('Time signal') 
subplot(2,2,3)
plot(ftx(1:floor(end/2)),theta1(1:floor(end/2)));  title('Phase of Fourier components')

t0 = -3;
y2 = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty2 = fft(y2);
dt = t(2)-t(1);
ftx = (0:1:length(t)-1)*(1/dt)/length(t);
theta2 = unwrap(angle(fty2));
subplot(2,2,2)
plot(t,y2);  title('Time signal') 
subplot(2,2,4)
plot(ftx(1:floor(end/2)),theta2(1:floor(end/2)));  title('Phase of Fourier components')

My claim is that the sequences y1 and y2 do not have a time-shift relationship that would correspond to the DFTs being related by a linear phase term. I think the easiest way to see this is to apply a linear phase term to the DFT of y1, take the iDFT of the result , and compare the results to fty2 and y2.

First, compute the number of samples in the shift from y1(t) to y2(t)

Nshift = 4/dt % numer of samples the signal was shifted

Nshift = 160

Multiply the DFT fty1 by the linear phase term.

fty2shiftedfromfty1 = fty1.*exp(1j*(ftx*dt*2*pi)*Nshift);

Compare the wrapped angle fty2 and fty2shiftedfromfty1. They do not look the same.

figure
plot(ftx,angle(exp(1j*theta2)),'-o',ftx,(angle(fty2shiftedfromfty1)),'-x')
xlim([0 20])

So, fty2 is not reated to fty1 by a linear phase term.

Compute the ifft of fty2shiftedfromfty1

y2shiftedfromfty1 = ifft(fty2shiftedfromfty1,'symmetric');

Compare it to y2

figure
plot(t,y2shiftedfromfty1 - y2)

So they are close, but not exactly the same, and apparently that small difference in the time domain samples makes a big difference in the figure-of-merit you're looking at.

To enforce a linear phase relationship between fty1 and fty2, we would create y2 as

y2shiftedfromy1 = circshift(y1,-Nshift);
fty2shiftedfromy1 = fft(y2shiftedfromy1);
figure
plot(t,y1,t,y2shiftedfromy1)

y2shiftedfromy1 looks very much like y2 above, but it's not the same:

figure
plot(t,y2shiftedfromy1-y2)

In fact, it is the same as y2shiftedfromfty1

figure
plot(t,y2shiftedfromy1 - y2shiftedfromfty1)

Verify the linear phase relationship

figure
plot(ftx,angle(fty2shiftedfromy1),'-o',ftx,angle(fty1.*exp(1j*(ftx*dt*2*pi)*Nshift)),'-x')

[Jiayun Liu]

Is this result expected or it is strange that 2 seemingly similar signals give different results? If we test this in a controlled way like this, we can explain that the 2 signals are not exactly the same but what happens when we are dealing with measured signals?

[Jiayun Liu]

Hi David,

I wrote a code to investigate where the phase gradient changes sign and it is consistent with exactly what you observed. Using your convention, positive gradient when t0>mid point and negative when t0<mid point.

[Paul]

"Is this result expected or it is strange that 2 seemingly similar signals give different results?"

I'm pretty sure that what's happening here can be explained as follows. As you noted, the signals y2 and y2shiftedfromfty1 are very nearly the same, with their difference plotted above and no larger than about 1e-5

y2 = y2shiftedfromfty1 + dy2

By linearity, we know that the DFT of y2 is the sum of the DFTs of the terms on the rhs. But for frequencies outside the region 9 < f < 11, the amplitude of the DFT of y2shiftedfromfty1 is very small as is the amplitude of the DFT of dy2, which is small everywhere. The phase of the sum of two complex numbers can be quite different from the phase of either if the the magnitudes of each are close enough, even if the magnitudes of each are both small. Note that in the region 9 < f < 11 where their magnitudes are large, the phase of y2 is essentially the same as that of y2shiftedfromfty1.

Can you further clarify the question about "measured signals"?

[Jiayun Liu]

Your explanation attributes the different gradient to the small difference between the two signals. But this doesn't explain why each signal has the phase gradient that they have except the gradients are different because they have a small but finite difference in the time domain.

I am referring to experimentally measured signal e.g. we measure and record the temperature using a thermocouple.

[Paul]

Sure, let's assume that we have a signal that is a sequence of temperature measurements from a thermocouple. What is your question about that signal? Is that signal somehow related to the signal, or variations of it, that we've been discussing up to this point?

[Jiayun Liu]

We saw that in your example 2 seemingly similar signals can have different phase gradient due to some very small difference although we do not know why that is the case or in other words we do not know why they have the phase gradients that they have.

Let's say that in my measurement, the phase gradient of the fft-ed temperature signal is all that I care about. Now we assume that we have 2 signals that look very much the same (just like in your example). I would believe that it is possible that they can have different phase gradient based on your example unless there is some reason your example doesn't apply to measured signal. If that is the cae, I don't see why the result is correct since the phase gradient describes the position of the signal with respect to the imaginary part of the fft kernel i.e. sin function.

To see why this is so, let's consider a simple cos function. Since the phase is atan(y/x) where x and y is the real and imaginary part of the fft signal, we can see that the phase goes to 0 when y=0 and this happens if which happens when mod(N,m)=0 where m is the length of 1 wavelength of the signal. We see that this summation is non zero if the signal is shifted by some . Hence, the phase gradient describes the position of the signal with respect to the first point of your signal based on the documentation in matlab fft help.

Therefore, it would be wrong to have 2 seemingly the same signals to exhibit different phase gradient when the phase gradient is really measuring the relative delay between the different sine waves in the signal (dispersion of the signal). From your example, we see that this is possible.

[Paul]

Yes, the example(s) showed that it is possible that the "phase gradient" of the DFTs are different for two seemingly the same signals for the reasons that I've already explained above: where the frequency samples of the DFTs end up relative to the DTFTs, and that the phase of the sum of complex numbers can be very different than the phase of either if both numbers have similar-ish magnitudes, even if they are very small. Here's another example that illustrates the second point. Since we are now talking about measurements, let's just add some noise to the signal that we've been talking about.

t = 0:1/40:10;  t0 = -7;
y1 = sin(2*pi*10*(t+t0)).*exp(-(t+t0).^2);
fty1 = fft(y1);
theta1 = angle(fty1);

dt = t(2)-t(1);
ftx = (0:1:length(theta1)-1)*(1/dt)/length(theta1);

Add some noise. This noise is quite small relative to the peak of y1, but is larger than the tails of y1. Is y2 "seemingly the same" as y1? Anyway

rng(100)
y2 = y1 + 1e-4*rand(size(y1));  
fty2 = fft(y2);
theta2 = angle(fty2);

figure
subplot(221)
plot(t,y1)
subplot(223)
plot(ftx,unwrap(theta1));xlim([0 20]),ylim([-20 60])
subplot(222)
plot(t,y2);
subplot(224);
plot(ftx,unwrap(theta2));xlim([0 20]),ylim([-20 60])

The unwrapped phases of y2 and y1 aren't the same, but at least kind of have the same trend. But, in the center of the frequency range where the magnitude of the DFTs are large, the phase is the same, because the DFT of the noise is small amplitude at all frequencies.

figure
plot(ftx,theta1,ftx,theta2,'o'),xlim([9 11])

Try running the code above with zero-paded fft's, e.g., to 8192 points, and see how the results change (you'll have to adjust some axis limits) and if that change makes sense.

[Jiayun Liu]

"I've already explained above: where the frequency samples of the DFTs end up relative to the DTFTs" Yeah I guess this is the best explanation. Thanks very much for the discussion :D