I have a problem with the angle function

I'm trying to plot the phase of a fft result using the functions angle and unwrap, but I noticed the result was weird, so I tried reducing the problem to the minimum. Now I'm just defining a phase P and plotting that vector P and the unwrapped angle of exp(i*P), so both plots should be the same, but the second one does something weird at the end having a maximum and then decreasing. What am I doing wrong?

w0 = 2.3607e15;
v=3e8;
lmax=835e-9;
lmin=765e-9;
num=255;
lstep=(lmax-lmin)/num;
l=lmin:lstep:lmax;

w=flip(2*pi*v./l);

TO = 0.8*10^-39;
P = TO/6*(w-w0).^3;
Ew = exp(1i*P);


plot(w,unwrap(angle(Ew)))
hold on
plot(w,P)

[the cyclist]

I don't have anything to add to the two answers, but note that it is happening at the very beginning of your curves as well:

w0 = 2.3607e15;
v=3e8;
lmax=835e-9;
lmin=765e-9;
num=255;
lstep=(lmax-lmin)/num;
l=lmin:lstep:lmax;

w=flip(2*pi*v./l);

TO = 0.8*10^-39;
P = TO/6*(w-w0).^3;
Ew = exp(1i*P);

% Plot just the first few points
plot(w(1:4),unwrap(angle(Ew(1:4))))

Answer 1

[Daniel]

Toward the end of the plot, your angle starts increasing by more than π between successive samples. unwrap just compresses the differences to lie within the range, so when the desired angle increase is greater than π, unwrap will roll that to a smaller (negative) value rather than the positive value you want.

The unwrapped angle starts to drift negative when the difference in angles exceeds π in the graphs below.

w0 = 2.3607e15;
v=3e8;
lmax=835e-9;
lmin=765e-9;
num=255;
lstep=(lmax-lmin)/num;
l=lmin:lstep:lmax;
w=flip(2*pi*v./l);
TO = 0.8*10^-39;
P = TO/6*(w-w0).^3;

subplot(2,1,1)
plot(diff(P))
hold on
plot([0 255],[pi pi])
legend('diff(P)','pi')
title('Angle difference in P')
subplot(2,1,2)
plot(P)
hold on
plot(unwrap(P))
legend('P','unwrap(P)')
title('Unwrapped version of P')

[FranM]

So when I don't have P, how could I plot the unwrapped phase of a vector Ew without having this problem?

[Daniel]

Unfortunately, "unwrapped phase" becomes ambiguous when the angle varies this quickly. It's an ambiguity that's intrinsic to how angle is measured. You can kind of get around it by collecting data that's more closely spaced, but that may or may not work as you wish depending on the application.

When I expect phase to follow a particular curve, sometimes I'll do the unwrapping on my own. For instance, I might try to optimize for a smooth second-order difference, rather than holding the first-order difference to the range like unwrap does. But it would probably help at this point to see what your original issue was, in the fft domain. That will give us a better idea what you're looking for and what sort of solutions (or explanations) will help.

[FranM]

The original issue is that I have a gaussian function (Ew) whosephase I can define by selecting the second and third order components of a Taylor expansion. I want to plot the amplitude and phase of its IFT function (Et).

w0 = 2.3607e15;
v = 3e8;
lmax = 835e-9;
lmin = 765e-9;
num = 255;
lstep = (lmax-lmin)/num;
l = lmin:lstep:lmax;

w = flip(2*pi*v./l);

cw = 1.2508e13;
Aw = exp(-(w-w0).^2/(2*cw^2));

% I can change these numbers to see how it affects the amplitude and phase of the IFT
SO = 23000*10^-30;
TO = 0.8*10^-39;

P=SO/2*(w-w0).^2+TO/6*(w-w0).^3;

Ew=Aw.*exp(1i*P);

Et=ifftshift(ifft(Ew));

figure
yyaxis left
plot(abs(Et).^2)
hold on
yyaxis right
plot(unwrap(angle(Et)))
hold off

[Daniel]

Thanks, the context helps.

Normally I would recommend taking more interior points. In this case, increasing num results in a faster-varying IFFT, due to the time-frequency duality involved. But you can increase the number of points in the IFFT directly by calling ifft(Ew,4096) (or whatever) rather than ifft(Ew). This gives you a finer time resolution, more interpolated points, and to your issue, decreases the phase difference so it's always within the range.

There's still a discontinuity around x=3000. I think this is because the phase shifts by π at that location, due to the abs function changing its behavior as the function crosses the origin. unwrap wraps to multiples of , so it will always show points where a function crosses through the origin as discontinuities. EDIT: Looking more closely at the behavior of the graph, this assumption is definitely not quite correct, but if you look at the spacing of the zero-crossings between real(Et) and imag(Et), something interesting definitely happens around 2750 (using num=255 and a 4096-point IFFT).

I'm not sure why the unwrapped phase here doesn't match P, but I would guess that it has to do with relative scaling as you go through the IFFT, and that the domain of P might be wider than the domain of angle(Et). But I might be wrong. I'm not familiar with the underlying mathematics here...

Hopefully that helps work through the mechanics of it. It looks like there are further problems to sort out after making sure phase was unwrapped as correctly as possible, but I don't think those can be attributed to phase unwrapping.

Answer 2

[Walter Roberson]

Your angle(Ew) oscillates around 0 so quickly that unwrap() ends up just adding more and more full cycles to try to keep up. Eventually you happen upon a change that can be explained in-phase as a descrease instead of needing a wrap, and at the point unwrap() switches to descreasing phases.