1
$\begingroup$

I have a system that I am simulating, and I want to add phase noise to the process. I didn't know how to do it, so I have been conducting an literature review.

I found a good reference with an "easy to understnad" example. This example is on page 145 and 146, where there are a pair of equations:

$$osc(t[n])=\sin(\omega_0 t[n] + \Delta \phi(t[n])$$ $$\Delta \phi(t[n])=\phi_{synch}(t[n])+\sum^n_{k=1}\phi_{acc}(t[k])$$

In these equations, $osc(t[n])$ is the oscillator output, $\omega_0$ is the frequency of oscillation, $\Delta \phi(t[n])$ is the phase noise, and $\phi_{synch}(t[n])$ and $\phi_{acc}(t[k])$ represent the synchronous and accumulating components of the phase deviations.

The authors used these equations in conjunction with some matlab code to publish a nice graph on page 146:

enter image description here

This is exactly what I need! However, I am confused by the matlab code, and by the exact nature of $\phi_{synch}(t[n])$ and $\phi_{acc}(t[k])$. Specifically, the authors say:

$\phi_{synch}(t[n])$ and $\phi_{acc}(t[k])$ are random variables with Gaussian distributions.

These are my questions:

How does one choose the variance and mean for the phase noise Gaussian distributions?

How do these values relate to the system itself, to Johnson-Nyquist thermal noise, and to the oscillator frequency $\omega_0$?

And what purpose do the filters in the code (shown below) serve?

Bonus quote from this fantastic book:

In order to derive an expression for $\Delta \phi(t[n])$, we note that an oscillator is an autonomous system in which phase deviations accumulate over time. However, the inherent accumulation of phase deviations does not affect the noise contributions of external sources that are added directly to the output of the oscillator , as shown in [4]. Hence ,the phase of the output signal of an oscillator can include a noise component which is not accumulated.This noise component can be considered synchronous with the output of the oscillator. Accumulating and synchronous phase deviations are related to accumulating and synchronous jitter...


Code from appendix A.4 page 173:

clc
close all
clear all

tstep = 1;
nsamples = 2^23
time = 0:tstep:nsamples -1;
period = 2^9;
ThreeSigma = period/20;
Sigma = ThreeSigma /3;
varjitter = Sigma^2;
Nperiods = length(time)/period +10;

jitterWhite = Sigma*randn(1,nsamples);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% FILTER DESIGN: Single POLE %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms z z1 outz inz p1;
fs = tstep/period;
s = 2 * fs * (z-1)/(z+1);
Hz = p1/(s+p1);
[N,D] = numden(Hz);
Cnum = coeffs(N, z);
Cden = coeffs(D, z);
Cnum1 = Cnum(1);
Cnum2 = Cnum(2);
Cden1 = Cden(1);
Cden2 = Cden(2);
fsval = 2*pi*fs; % sampling frequency in the discrete time domain
p1val = fsval /100; % frequency of the pole
Cnum1valextra = double(subs(Cnum1 , {p1, fs},{p1val , fsval}));
Cnum2valextra = double(subs(Cnum2 , {p1, fs},{p1val , fsval}));
Cden1valextra = double(subs(Cden1 , {p1, fs},{p1val , fsval}));
Cden2valextra = double(subs(Cden2 , {p1, fs},{p1val , fsval}));
Aextra = (Cnum2valextra/ Cden2valextra);
Bextra = (Cnum1valextra/ Cden2valextra);
Cextra = - (Cden1valextra/ Cden2valextra);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% FILTER: Implementation %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jitterFiltered = zeros(1,length(jitterWhite));
for n = 2:1: length(jitterWhite)
   jitterFiltered(n) = (Cnum2valextra/ Cden2valextra) * jitterWhite(n) + (Cnum1valextra/ Cden2valextra) * jitterWhite(n-1) - (Cden1valextra/ Cden2valextra) * jitterFiltered(n-1) ;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% SECOND FILTER: Design and Implementation %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jitterWhiteNew = Sigma*randn(1,nsamples);
syms z z1 outz inz p1;
fs = tstep/period;
s = 2 * fs * (z-1)/(z+1);
Hz = p1/(s+p1);
[N,D] = numden(Hz);
Cnum = coeffs(N, z);
Cden = coeffs(D, z);
Cnum1 = Cnum(1);
Cnum2 = Cnum(2);
Cden1 = Cden(1);
Cden2 = Cden(2);
fsval = 2*pi*fs; % sampling frequency in the discrete time domain
p1val = fsval /10000; % frequency of the pole (Hz)
Cnum1valextra = double(subs(Cnum1 , {p1, fs},{p1val , fsval}));
Cnum2valextra = double(subs(Cnum2 , {p1, fs},{p1val , fsval}));
Cden1valextra = double(subs(Cden1 , {p1, fs},{p1val , fsval}));
Cden2valextra = double(subs(Cden2 , {p1, fs},{p1val , fsval}));
Aextra = (Cnum2valextra/ Cden2valextra);
Bextra = (Cnum1valextra/ Cden2valextra);
Cextra = - (Cden1valextra/ Cden2valextra);
jitterFilteredNew = zeros(1,length(jitterWhite));
for n = 2:1: length(jitterWhite)
   jitterFilteredNew(n) = (Cnum2valextra/ Cden2valextra) * jitterWhiteNew(n) + (Cnum1valextra/Cden2valextra) * jitterWhiteNew(n-1) - (Cden1valextra/ Cden2valextra) * jitterFilteredNew(n-1) ;
end
jitter = 0.01*( 0.01* cumsum(jitterFilteredNew) + 0.01* jitterFiltered );
sinsignalNoisy = sin(2*pi*(1./ period).*time+ jitter);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.