Strange disagreement between analytical phase and unwrap(angle())

I consider a simple pulse sequence composed of four discete pulses that start at sampel and include zeros interleaved between the pulses:

Matlab uses the following convention for the discrete Fourier Tranform in function fft.m:

I compute this by hand:

using the Dirichlet kernel to simplify the expression for the amplitude and phase; here , , , and . I'd like to numerically compute the amplitude and phase of my analytic solution with Matlab's output. I need to understand the discprepancy between Matlab's output phase and the analytical version ; I emphasize that the "analytical" version is the discrete Fourier transform version, using the Matlab convention. The following code shows that the fourier domain signal appears to agree with my solution (above), but that the phase output by Matlab, separately, does not:

%Time series parameters
L       = 4;%number of pulses
dn      = 20;%zeros between pulses
N	= 400;%total length of time series
n0	= 2^4;%first non-zero pulse
Ampl    = 15; %pulse sequence amplitude

%Define signal and it's spectral form
k           = 1:N;
k           = k(:);
signal      = zeros(N,1);
ind         = n0 + dn*(0:(L-1)); %non-zero pulse indices
signal(ind) = Ampl; %asign at non-zero pulse indices
Xk          = fft(signal); %fft of signal using convention I wrote above

%compare MATLAB output to analytical solutions

%MATLAB ampltiude and phase
absXk       = abs(Xk); %magnitude of signal
phiXk       = unwrap(angle(Xk)); %phase of signal

%analytical amplitude and phase
Xkabs       = abs( Ampl*sin((pi*L/N)*(k-1)*dn)./sin((pi/N)*(k-1)*dn) ); %magnitude above
Xkabs(1)    = Ampl*L; %divide by zero correction
Xkphi       = -(2*pi/N)*(k-1)*(n0 - 1 + 1/2*(L-1)*dn); %analytical phase term from equation above

%Compare amplitude (they match)
figure;
plot(k, absXk,'-k','linewidth',6); hold on;
plot(k, Xkabs,'--r','linewidth',2);
title('Compare Amplitude Spectrum'); %Checks!
axis square;

%Compare phase (they do not match)
figure;
plot(k, phiXk, '-k','linewidth',6); hold on;
plot(k, Xkphi,'--r','linewidth',4);
title('Compare Phase'); %Disagreement
axis square;

%Compare total signals (they match)
figure;
plot(Xk,'-k','linewidth',6); hold on;
plot(Xkabs.*exp(1j*Xkphi),'-r','linewidth',2); hold on;
title('Compare Complex Signals'); %Agreement
axis square;

The three figures from each plot illustrate the issue with unwrap( angle ( ) ), which I show below. I'm unclear what the problem is. Is it circularity in the phase? I require that the phase output by Matlab can be modified to agree with a model.

Thank you.

Amplitude, phase, and total signal shown in order:

Signal:

Answer 1

[David Goodmanson]

Hello Joshua,

There a a couple of issues with the analytic calculation.

Xkabs   = abs( Ampl*sin((pi*L/N)*(k-1)*dn)./sin((pi/N)*(k-1)*dn) ); %magnitude above
Xkphi       = -(2*pi/N)*(k-1)*(n0 - 1 + 1/2*(L-1)*dn);

[1] Either or both of the sine terms can be negative. When one is positive and one is negative, there is a phase factor of +-pi that is not taken into account by the Xkphi expression. It probably can be done to figure out when the angle is +pi and when it is -pi, but it's much easier to use an analytic expression that stops one step short of deriving a ratio of sine terms.

[2] Since dn = 20 evenly divides N = 400, there are values of k such that the argument of one or both of the sines are exact multiples of pi, so that the sines are zero. Whether numerically the sines are exactly zero is a whole different question. But when one or both of the sines differ from zero by even a tiny amount, the derived angle can vary by a lot in an unintended way. Rather than deal with that issue, it's easier to sidestep the whole thing it by changing N from 400 to 399, so that dn does not divide N.

The resulting calculation shows agreement between the fft and analytic results. Since I was not keen on trying to remember which was which between a variable called Xkabs and another one called absXk, I used Z for the fft result and A for the analytic result.

L   = 4;       % number of pulses
dn  = 20;      % zeros between pulses
N	= 399;     % total length of time series
n0	= 2^4;     % first non-zero pulse
Ampl    = 15;  % pulse sequence amplitude
% Define signal and its spectral form
k           = 1:N;
k           = k(:);
signal      = zeros(N,1);
ind         = n0 + dn*(0:(L-1)); % non-zero pulse indices
signal(ind) = Ampl;              % assign at non-zero pulse indices

Z = fft(signal);         % fft of signal using convention written above
% compare MATLAB output to analytical solutions
% analytical amplitude and phase, but not converted to sines
A  = Ampl*exp(-2*pi*i/N*(k-1)*(n0-1))  ...
    .*(1-exp(-2*pi*i*L/N*(k-1)*dn))./(1-exp(-2*pi*i/N*(k-1)*dn))
figure(1)
plot(k,abs(Z),k,abs(A)+2)      % overlay, add offset for visual purposes
grid on
figure(2)
plot(k,unwrap(angle(Z)),k,unwrap(angle(A))+2) % overlay, add offset for visual purposes
grid on  

[Joshua Carmichael]

Well, color me a bit embarrassed. I was almost certain I had the analytical solution considered carefully, and you demonstrated I had not. I appreciate it!

I did a cut and paste and see that the plot of the phase that you generated does indeed agree with that of the analytical solution. If I could buy you a beer (or coffee), I would! Thank you. -Josh

[David Goodmanson]

Beer is better! Like many things these days it will have to be virtual, but thanks very much for the thought.