**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

14 views (last 30 days)

Show older comments

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))

### Accepted Answer

William Rose
on 5 Feb 2024

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).

##### 19 Comments

Janko
on 5 Feb 2024

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
on 5 Feb 2024

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
on 6 Feb 2024

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
on 6 Feb 2024

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
on 6 Feb 2024

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
on 7 Feb 2024

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
on 7 Feb 2024

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:

Janko
on 7 Feb 2024

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
on 8 Feb 2024

William Rose
on 8 Feb 2024

William Rose
on 8 Feb 2024

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
on 8 Feb 2024

Edited: Janko
on 8 Feb 2024

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
on 9 Feb 2024

Paul
on 9 Feb 2024

Edited: Paul
on 9 Feb 2024

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
on 9 Feb 2024

Edited: Paul
on 9 Feb 2024

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
on 11 Feb 2024

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
on 13 Feb 2024

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

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)