Why the calculated phase differs so much from the real phase of the signal?

Hi guys,

I have a small question regarding my code below. In my code I create a complex signal E1 solely consisting of the phase noise phi. My question is how to obtain the correct values for phi? My idea was to create the FFT of my signal E1 and then simply calculating the phase of the FFT by using the angle() function.

But plotting phi and phi_ang (phase obtain with angle function) shows that both phase vectors differ enormously.

I hope you can help me with my question.

Kind regards

Janko

The code:

clear
clc
N = 10000000;
wn =  randn(N,1);
phi = cumsum(wn);
E1 =  exp(1i*phi);
E1_fft = fftshift(fft(E1));
phi_ang = angle(E1_fft);
plot(phi)
hold on
plot(unwrap(phi_ang))

Answer 1

[William Rose]

@Janko,

The DFT expresses a signal x(k) as a sum of sinusoids, each of which is sampled at regular intervals. The phase angles of the elements of the DFT are the phases of the individual sinusoidal components, at t=0.

Your signal can be thought of as a regularly sampled signal whose phase advances by a random amount (forward or backward) with each timestep. Thus phi is a random walk with Gaussian steps. The phase angles (phi) of the random walk are not equal to the phase angles of the sinusoids which, when added up, will equal x(k).

[William Rose]

@Janko,

your signal

where is a random walk

The DFT of x(n) is

Then

[Janko]

Hi @William Rose,

thank you very much for the detailed answer!

My main aim is to reconstruct the phase of the signal as good as possible. From your answer I conclude that increasing the fft size will give me better results?

Besides, are there better methods to obtain the exact phase information of a random input signal?

[William Rose]

@Janko,

Increasing the size of the Fourier transform will not help. The phases of the sinusoids that make up a random signal are themselves random. May I ask why you wan these phase angles? It is hard for me to imagine how they will be useful to you.

N = 100;
wn =  randn(N,1);
phi = cumsum(wn);
E1 =  exp(1i*phi);
E1_fft = fftshift(fft(E1));
phi_ang = angle(E1_fft);
w=2*pi*(-N/2:N/2-1)/N;
figure;
subplot(211), plot(0:N-1,phi,'-r.');
grid on; xlabel('n=Time'); ylabel('\phi(n)'); title('x(n)=e^{i*\phi(n)}');
subplot(212), plot(w,unwrap(phi_ang),'-b.'); 
grid on; xlabel('\omega =Frequency'); ylabel('Phase(X(k))'); 
xlim([-pi,pi]); title('X(k)');

OK

[Janko]

@William Rose

thanks again for the detailed answer!

Regarding your question, my aim is to obtain information of the signal phase noise. Therefore I have to somehow obtain (approximate) the exact phase vector of the signal. The PSDs of phi and phi_ang differ unfortunetaly and therefore I cant use phi_ang for further calculations.

To obtain the linewidth of a signal one method would be to use the one-sided frequency noise spectrum of the signal. The contribution of white noise, which determines the intrinsic linewidth of the signal, should be constant. The code below can be used to calculate and plot the phase & frequency noise spectra. I have also added 1/f (or pink) frequency noise to the signal in order to see the effects (slopes) of the different contributions in the output spectra.

N = 100000;
wn = randn(N,1);
pn =  pinknoise(N);
phi = cumsum(wn+ pn);
E1 = exp(1i*phi);
E1_fft = fftshift(fft(E1));
phi_ang = unwrap(angle(E1_fft));
omega =2*pi*(0:N/2-1)/N;
S_phi = 2*abs(fftshift(fft(phi))).^2;
S_phi_ang = 2*abs(fftshift(fft(phi_ang))).^2;
S_phi(1:N/2) = [];
S_phi_ang(1:N/2) = [];
figure
loglog(omega,S_phi)
hold on
loglog(omega, S_phi_ang)
freq = 1/(2*pi) * diff(phi) *(2*omega(end)/2/pi);
freq_ang = 1/(2*pi) * diff(phi_ang) *(2*omega(end)/2/pi);
freq = [freq; 1/(2*pi)*(phi(end)-phi(1))*(2*omega(end)/2/pi)];
freq_ang = [freq_ang; 1/(2*pi) *(phi_ang(end)-phi_ang(1))*(2*omega(end)/2/pi)];
S_freq = 2* abs(fftshift(fft(freq))).^2;
S_freq_ang = 2* abs(fftshift(fft(freq_ang))).^2;
S_freq(1:N/2) = [];
S_freq_ang(1:N/2) = [];
figure
loglog(omega,S_freq)
hold on
loglog(omega, S_freq_ang)

The code for the pink noise contribution:

% +------------------------------------------------------+
% |   Pink Noise Generation with MATLAB Implementation   | 
% |                                                      |
% | Author: Ph.D. Eng. Hristo Zhivomirov        07/31/13 | 
% +------------------------------------------------------+
%
% function: x = pinknoise(N)
%
% Input:
% N - number of samples to be returned in the noise column vector
% alpha - PSD spectral slope
% 
% Output:
% x - column vector of pink noise samples with unity  
%     standard deviation and zero mean value 
%
% The function generates a column vector of pink (flicker) noise 
% samples. In terms of power at a constant bandwidth, the pink  
% noise falls off by -3 dB/oct i.e., -10 dB/dec.
function x = pinknoise2(N)
% input validation
validateattributes(N, {'double'}, ...
                      {'scalar', 'integer', 'nonnan', 'finite'}, ...
                      '', 'N', 1)
% set the PSD slope
alpha = -1; 
% convert from PSD (power specral density) slope 
% to ASD (amplitude spectral density) slope
alpha = alpha/2;
% generate AWGN signal
%x =  lpfilter(randn(1,N)*sqrt(variance), N/iOverSamp, iOverSamp);
x =  randn(1,N);
% calculate the number of unique fft points
NumUniquePts = ceil((N+1)/2);
% take fft of x
X = fft(x);
% fft is symmetric, throw away the second half
X = X(1:NumUniquePts);
% prepare a vector with frequency indexes 
n = 1:NumUniquePts;
% manipulate the left half of the spectrum so the spectral 
% amplitudes are proportional to the frequency by factor f^alpha
X = X.*(n.^alpha);
% perform ifft
if rem(N, 2)	% odd N excludes Nyquist point 
    % reconstruct the whole spectrum
    X = [X conj(X(end:-1:2))];
    
    % take ifft of X
    x = real(ifft(X));   
else            % even N includes Nyquist point  
    % reconstruct the whole spectrum
    X = [X conj(X(end-1:-1:2))];
    
    % take ifft of X
    x = real(ifft(X));  
end
% ensure zero mean value and unity standard deviation 
x = x - mean(x);
x = x/std(x, 1);
x = x(:);
end

[Janko]

After further investigation I realized that calculating the FFT is completely unnecessary for me. I should simply calculate the angle of my signal by calculating unwrap(angle(E1)). But then again the calculated an original phase vectors differ very much and I don't understand why.

Updated code:

N = 10000000;
wn = randn(N,1);
pn =  pinknoise(N);
phi = cumsum(wn+ pn);
E1 = exp(1i*phi);
E1_fft = fftshift(fft(E1));
phi_ang = unwrap(angle(E1));

figure
plot(phi)
hold on
plot(phi_ang)

omega =2*pi*(0:N/2-1)/N;
S_phi = 2*abs(fftshift(fft(phi))).^2;
S_phi_ang = 2*abs(fftshift(fft(phi_ang))).^2;
S_phi(1:N/2) = [];
S_phi_ang(1:N/2) = [];
figure
loglog(omega,S_phi)
hold on
loglog(omega, S_phi_ang)
freq = 1/(2*pi) * diff(phi) *(2*omega(end)/2/pi);
freq_ang = 1/(2*pi) * diff(phi_ang) *(2*omega(end)/2/pi);
freq = [freq; 1/(2*pi)*(phi(end)-phi(1))*(2*omega(end)/2/pi)];
freq_ang = [freq_ang; 1/(2*pi) *(phi_ang(end)-phi_ang(1))*(2*omega(end)/2/pi)];
S_freq = 2* abs(fftshift(fft(freq))).^2;
S_freq_ang = 2* abs(fftshift(fft(freq_ang))).^2;
S_freq(1:N/2) = [];
S_freq_ang(1:N/2) = [];
figure
loglog(omega,S_freq)
hold on
loglog(omega, S_freq_ang)

[William Rose]

@Janko,

You say "the calculated and original phase vectors differ very much and I don't understand why". Is that a question? Do you mean "Why doesn't phase(X(1))=phi(1), and phase(X(2))=phi(2), etc?" I tried to answer that before, but I did not do a great job, and I could do better, if that is still your question.

You also wrote

"my aim is to obtain information of the signal phase noise. Therefore I have to somehow obtain (approximate) the exact phase vector of the signal."

I assume you mean the phase of X(k), where X(k)=DFT of x(n). I still don't understand what you expect to learn from the vector of phase angles. Are you studying errors that occur in phase shift keying? It does not seem so.

You write: "The PSDs of phi and phi_ang differ unfortunetaly and therefore I cant use phi_ang for further calculations."

Those PSDs are unusual quantities to compute. Even if they were equal, what would you do with them?

Your code includes

freq = 1/(2*pi) * diff(phi) *(2*omega(end)/2/pi);

This appears to be an attempt to compute

or something proportional. That equation is useful for understanding a steady sine or cosine, or a wave packet, or an amplitude-modulated carrier wave. But it is probably not appropriate for understanding the sequence , where is a random sequence. And the power spectrum of that vector of frequencies, which you also compute, has no usefulness that I can think of.

I'm no expert on white and pink noise, but I have some familiarity with both. I have generated pink noise and white noise waveforms and have used both in experiments which were published in good journals. But I do not understand what you are doing, or trying to do, in the code you have posted. I encourage you to consult with somone knowledgeable in signal processing, to assure your effort is well-spent.

[Paul]

Hi Janko,

"I should simply calculate the angle of my signal by calculating unwrap(angle(E1)). But then again the calculated an original phase vectors differ very much and I don't understand why."

They won't match because phi can change by much more than +-pi from one sample to the next.

rng(100)
N = 10000000;
wn = randn(N,1);
pn =  pinknoise(N);
figure
plot(wn + pn)

So the unwrapped angle won't match the original angle after the first instance where such a large jump occurs. Don't need to go through the exp function to see this effect.

phi = cumsum(wn+ pn);
figure
plot(phi),hold on, plot(unwrap(phi))

Zoom in at the start

figure
plot(phi),hold on, plot(unwrap(phi))
xlim([0 1000])

figure
plot(abs(diff(phi))),xlim([0 1000]),yline(pi)

Taking a round trip through exp doesn't change anything

E1 = exp(1i*phi);
%E1_fft = fftshift(fft(E1));
phi_ang = unwrap(angle(E1));
figure
plot(phi)
hold on
plot(phi_ang)

[Janko]

Hi @William Rose

thank you very much again!

My main aim is to obtain is to approximate the phase noise within one signal as good as possible in order to calculate the corresponding frequency noise by simply differentiating it. Afterwards I want to use the obtained frequency noise to compute its PSD, so I get the corresponding Frequency Noise Spectrum. Within this spectrum there are various parts with different slopes, which are related to the different noise contributions within the signal. White Noise contributions lead to constant part within the Frequency Noise Spectrum and thus is proportional to the intrinsic linewidth of the signal. So all in all I simply want to calculate the linewidth of a signal by using time-domain data measured with an oscilloscope.

Basically, I want to recreate of what's described within this paper:

https://opg.optica.org/directpdfaccess/8cf27d94-951d-479b-bc9dd3aba7009582_227860/oe-20-5-5291.pdf?da=1&id=227860&seq=0&mobile=no

[Janko]

Hi @Paul

thank you very much aswell, now I understand why both phases differ so much.

Can I wrap my phi to 2 Pi and unwrap it and do something like this:

N = 1000000;
wn = randn(N,1);
pn =  pinknoise(N);
phi = unwrap(wrapTo2Pi(cumsum(wn+ pn)));
E1 = exp(1i*phi);
E1_fft = fftshift(fft(E1));
phi_ang = unwrap(angle(E1));

This way both phases are now the same. Is this legitimate?

[Paul]

After reading this thread I'm sill not sure what the goal is, so I can't say what's legitimate or not relative to the intent. I was only trying to explain why phi and phi_ang don't match.

[William Rose]

@Janko, I share @Paul's sentiment.

If you want understand the noise spectrum, just compute the power spectrum, with fft(), or with pwelch(), or maybe with an ARMA model. Don't make things more complicated than necessary. Don't try to estimate the phases and differentiate those phases to get frequencies. Maybe I'll think differently if I read the paper you read - but the link you gave yields only "PDF no longer available on the URL requested". An example, which I mention only to illustrate that I have some practical experience: I analyzed muscle force from normal subjects and subjects with tremor. Tremor manifests itself as a sharp peak in the force spectrum. I estimated and removed white and pink noise components in the power spectrum, in order to get the tremor spectrum "by itself". I did not do any of the phase angle stuff which you are doing.

[William Rose]

You said in your most recent comment "White Noise ... is proportional to the intrinsic linewidth of the signal", and you said something similar earlier. I am not convinced that that is right. Consider this counter-example: two signal, identical except the white noise power is 4 time larger in one. And yet the linewidth is the same, as the plot shows. That seems to contradict your statement.

N=2^14;                 % number of points
t=0:N-1;                % time vector
fs=(t(end)-t(1))/(N-1); % sampling rate
wn1=randn(1,N);         % white noise, variance=1
wn2=2*randn(1,N);       % white noise, variance=4
f0=0.2;                 % oscillation frequency (Hz)
x1=cos(2*pi*f0*t)+wn1;  % signal x1
x2=cos(2*pi*f0*t)+wn2;  % signal x2
% COmpute power spectral densities
[pxx1,f1] = pwelch(x1,hann(256,'periodic'),128,256,fs);
[pxx2,f2] = pwelch(x2,hann(256,'periodic'),128,256,fs);
% Plot results
figure; plot(f1,pxx1,'-r',f2,pxx2,'-b');
grid on; xlabel('Frequency (Hz)'); ylabel('P.S.D.');
legend('PSD(x1)','PSD(x2)');

The linewidths are the same, even though the white noise power is 4 times greater in signal x2 (blue).

[Janko]

Hi @William Rose,

I meant "white frequency noise" and not "additive white gaussian noise". The white noise you have added in the code above introduces basically an offset in the amplitude of the signal, thus changing its power. The white frequency noise I am speaking about is cumulated over time and creates phase noise. Therefore you have to take the cumulative sum of the white noise contributions (creates a random walk process) and place those within the arguments of the cos functions.

It should look like this:

N=2^14;                 % number of points
t=0:N-1;                % time vector
fs=10*(t(end)-t(1))/(N-1); % sampling rate
wn1=randn(1,N);         % white noise, variance=1
wn2=4*randn(1,N);       % white noise, variance=4
f0=0.2;                 % oscillation frequency (Hz)
x1=cos(2*pi*f0*t+ cumsum(wn1)*1/fs);  % signal x1
x2=cos(2*pi*f0*t+ cumsum(wn2)*1/fs);  % signal x2
% COmpute power spectral densities
[pxx1,f1] = pwelch(x1,hann(256,'periodic'),128,256,fs);
[pxx2,f2] = pwelch(x2,hann(256,'periodic'),128,256,fs);
% Plot results
figure; plot(f1,pxx1,'-r',f2,pxx2,'-b');
grid on; xlabel('Frequency (Hz)'); ylabel('P.S.D.');
legend('PSD(x1)','PSD(x2)');

See how the linewidth increases when increasing the underlying white frequency noise?

The introduced phase noise does not change the power of the signal (contrary to AWGN), instead the power is distributed over a broader frequency range.

[William Rose]

@Janko,Thank you for introducing me to frequency noise. I have not heard of it before Here is a short summary which I found helpful. I will re-look at your comments above. Please send a working link to the paper whose results you wish to recreate.

[Paul]

This example threw me off when I saw that f0 = 0.2, but the peak in the plots is at 2Hz (no 0.2 Hz). Then I realized why I was confused

N=2^14;                 % number of points

Is there a reason the time vector and fs are defined like this? I've not seen this approach, and it seems like a peculiar way to get fs = (10*(N-1)/(N-1)

t=0:N-1;                % time vector
fs=10*(t(end)-t(1))/(N-1) % sampling rate

fs = 10

Why not

fs = 10;        % sampling rate
t = (0:N-1)/fs; % time vector
wn1=randn(1,N);         % white noise, variance=1
wn2=4*randn(1,N);       % white noise, variance=4

Change f0

%f0=0.2;                 % oscillation frequency (Hz)
f0=2;

x1=cos(2*pi*f0*t+ cumsum(wn1)*1/fs);  % signal x1
x2=cos(2*pi*f0*t+ cumsum(wn2)*1/fs);  % signal x2
% COmpute power spectral densities
[pxx1,f1] = pwelch(x1,hann(256,'periodic'),128,256,fs);
[pxx2,f2] = pwelch(x2,hann(256,'periodic'),128,256,fs);
% Plot results
figure; plot(f1,pxx1,'-r',f2,pxx2,'-b');
grid on; xlabel('Frequency (Hz)'); ylabel('P.S.D.');
legend('PSD(x1)','PSD(x2)');

[Paul]

If the intial goal is to estimate the additive phase noise to a cos signal, maybe hilbert is what you're looking for.

N = 2^14;                 % number of points

fs = 10;        % sampling rate
t = (0:N-1)/fs; % time vector

rng(100)
wn1 = randn(1,N);         % white noise, variance=1

Use the same sequence for apples-to-apples comparison

%wn2 = 4*randn(1,N);       % white noise, variance=4
wn2 = 4*wn1;

f0 = 2;

phi1 = cumsum(wn1)*1/fs;
x1   = cos(2*pi*f0*t + phi1);  % signal x1

phi2 = cumsum(wn2)*1/fs;
x2   = cos(2*pi*f0*t+ phi2);  % signal x2

Reconstruct the phase noise and overlay with the true phase noise.

X1 = hilbert(x1);
figure
plot(t,unwrap(atan2(imag(X1),real(X1))) - 2*pi*f0*t)
hold on
plot(t,phi1)
xlabel('Time');ylabel('rad')

Plot the difference

figure
plot(t,unwrap(atan2(imag(X1),real(X1)))-2*pi*f0*t - phi1)
xlabel('Time');ylabel('rad')

Do the same for x2

X2 = hilbert(x2);
figure
plot(t,unwrap(atan2(imag(X2),real(X2)))-2*pi*f0*t)
hold on
plot(t,phi2)
xlabel('Time');ylabel('rad')

figure
plot(t,unwrap(atan2(imag(X2),real(X2))) - 2*pi*f0*t - phi2)
xlabel('Time');ylabel('rad')

Based on the very extremely little I know about the Hilbert transform, I'm kind of surprised it worked as well as it did for x2. But the reconstruction with x2 is not as good as it is for x1, presumably because the factor of 4 on the phase noise is effectively increasing dphi/dt, in general, and pushing it closer to 2*pi*f0.

Don't know how well this approach will work in general, just thought it might be an interesting avenue to pursue. Your mileage may vary.

[William Rose]

@Janko,

I am very impressed with @Paul's ideas and knowledge.

Here is a script that implements frequency noise, using the random walk in phase that was your original idea. The (detrended, so it has zero mean value and zero trend) random walk is added to a phase angle that increases linearly with time. Your original example did not include the linearly increasing component of phase.

The script makes four figures, below.

Above: Phase (top) and detrended phase (bottom) versus time. The three phase plots in the top panel are indistinguishable, because the phase noise is very small copmpared to the linear trend in phase. Note difference in y-axis scale for the upper and lower panels.

Above: Power spectra of phase (top) and power spectra of instantaneous frequency (bottom). The phase spectra, for the phases with noise, are pink (power proportional to 1/f^2). The instantaneous frequency spectra, for the noisy frequencies, are white.

Above: Instantaneous frequency versus time. The mean frequency is 250 kHz for all three signals.

Above: Power spectrum of three signals, two of which have frequency noise added. Upper panel has al linear vertical axis; lower panel has a logarithmic vertical axis. Frequency noise causes broadening of the spectral peak.

[Janko]

Hi @Paul and @William Rose

thank you very much for the help and the detailed answers!

I will for sure look into the hilbert transformation!

And thank you for providing me the frequencyNoise script!

Here is the link to the paper I've mentioned above: https://doi.org/10.1364/OE.20.005291

(Kazuro Kikuchi, "Characterization of semiconductor-laser phase noise and estimation of bit-error rate performance with low-speed offline digital coherent receivers," Opt. Express 20, 5291-5302 (2012))

Thank you again for your time!

Kind regards

Janko

[William Rose]

@Janko, you're welcome. Good luck with your work.