(7 votes, average: 4.00 out of 5)

# Modeling phase noise (frequency domain approach)

by on September 30, 2012

In typical wireless system simulations, there is  a need to model the phase noise profile of the local oscillator. For eg, the phase noise profile of the oscillator can be of the shape described in the post on Phase Noise Power Spectral Density to Jitter. While looking around for example Matlab code, found two references [1, 2] which uses the approach of defining the phase noise profile in frequency domain, and then using ifft() to convert to the time domain samples. This post gives a brief overview of the modeling and provides an example Matlab/Octave code.

## Modeling

a) Assume a system with sampling frequency $f_s$ Hz and having $N$  samples. In frequency domain we can define $\[-\frac{f_s}{2}, \frac{f_s}{2}\)$ in steps of $\Delta f=\frac{f_s}{N}$ Hz.

b) Consider a phase noise profile defined as follows :

 freq PSD, dBc/Hz 0 -65 1kHz -65 10kHz -95 100kHz -115 1MHz -125 10MHz -125

Table : Example phase noise profile

From the phase noise profile, using linear interpolation (in log10 of the frequency axis), to find the phase noise power spectral density for frequencies from $\[-\frac{f_s}{2}, \frac{f_s}{2}\)$ in steps of $\Delta f=\frac{f_s}{N}$ Hz.

c) Generate a white Gaussian noise sample of length $N$ and scale it with the phase noise power spectral density

(Scaling by $\sqrt{\Delta f}$ is to normalize the resolution bandwidth to unity)

d) Use ifft() to find the time domain samples.

e) On the rel samples obtained from step (d), take $\e$${jx}$$$ to form the time domain phase noise samples.

Note :

When x is small, $\e$${jx}$$\approx 1 + jx$.

## Example Matlab script

% Script for simulating the an example phase noise profile
% ----------------------------------------------------------

clear all; close all;
fs_Hz    = 20e6;
N        = 10^5;
nIter    = 100;

% phase noise profile
psd_f_hz           = [  0  1e3  1e4  1e5  1e6  10e6];
psd_val_dbc_per_hz = [-65  -65  -95 -115 -125  -125];

% defining the frequency vector
freq_v_Hz = [0:N/2]/N*fs_Hz;
delta_f   = fs_Hz/N;
slope    = [psd_val_dbc_per_hz(2:end) - psd_val_dbc_per_hz(1:end-1) ]./...
(log10(psd_f_hz(2:end)) - log10(psd_f_hz(1:end-1)));
constant = 10.^(psd_val_dbc_per_hz(1:end-1)/10).* ...
(psd_f_hz(1:end-1).^(-slope/10));
integral = constant.*(psd_f_hz(2:end).^(1+slope/10) - ...
psd_f_hz(1:end-1).^(1+slope/10) )./(1+slope/10);
%% finding the rms jitter
% finding index with slope == -10
idx = find(slope==-10);
integral(idx) = constant(idx).*(log(psd_f_hz(idx+1)) - log(psd_f_hz(idx)));

% interpolating the phase noise psd values
psd_ssb_dB = -Inf*ones(1,N/2+1); % from [0:N/2]
for ii=1:length(psd_f_hz)-1
[tt1 fl_idx ] = (min(abs(psd_f_hz(ii) - freq_v_Hz)));
[tt2 fr_idx ] = (min(abs(psd_f_hz(ii+1) - freq_v_Hz)));
fvec = [freq_v_Hz(fl_idx):delta_f:freq_v_Hz(fr_idx)];
pvec = slope(ii)*log10(fvec+eps) + psd_val_dbc_per_hz(ii) - slope(ii)*log10(psd_f_hz(ii)+eps);
psd_ssb_dB(fl_idx:fr_idx) = pvec;
end

% forming the full vector from [-N/2:-1 0:N/2-1 ]/N*fs_Hz
psd_dB                  = -Inf*ones(1,N);
psd_dB([-N/2:-1]+N/2+1) = psd_ssb_dB([N/2+1:-1:2]);
psd_dB([0:N/2-1]+N/2+1) = psd_ssb_dB(1:N/2);

psd_linear = 10.^(psd_dB/20);

for (jj = 1:nIter)
% defining frequency vector
phase_noise_freq       = 1/sqrt(2)*(randn(1,N) + j*randn(1,N));
phase_noise_freq_scale = N*sqrt(delta_f)*phase_noise_freq;
phase_noise_freq_psd   = phase_noise_freq_scale .*psd_linear;

% converting to time domain
phase_noise_td      = ifft(fftshift(phase_noise_freq_psd));
pn_td               = exp(j*(sqrt(2)*real(phase_noise_td)));

% for estimating jitter and plotting
pn_without_carrier  = (pn_td - 1);

hF           = 1/(N*sqrt(delta_f))*fft(pn_without_carrier,N);
hFPwr(jj,:)  = hF.*conj(hF);
end

figure
semilogx( [-N/2:N/2-1]/N*fs_Hz, 10*log10(fftshift(mean(hFPwr))), 'r^-' );
hold on;grid on;
semilogx([0:N/2]/N*fs_Hz,psd_ssb_dB,'mp-');
semilogx(psd_f_hz,psd_val_dbc_per_hz,'bs-');
legend('est-freq-response','original','interpolated');
xlabel('freq, Hz'); ylabel('dBc/Hz');
axis([1 10e6 -140 -50]); title(title_str);

Figure : Example phase noise profile (expected and simulated)

## Summary

The above approach seems to allow a way to model an arbitrary phase noise power spectral density. However, the fact that this approach needs a large ifft() of length $N$can potentially slow down the simulation.

## References

[1] Phase Noise by Alex Bar-Guy, 27 Oct 2005 (Updated 08 Dec 2005) http://www.mathworks.com/matlabcentral/fileexchange/8844-phase-noise

[2]  Baseband-equivalent phase noise model, Submitted by Markus Nentwig on Dec 18 2011

http://www.dsprelated.com/showcode/246.php

D id you like this article? Make sure that you do not miss a new article by subscribing to RSS feed OR subscribing to e-mail newsletter. Note: Subscribing via e-mail entitles you to download the free e-Book on BER of BPSK/QPSK/16QAM/16PSK in AWGN.

yuanbaicai February 25, 2015 at 5:43 am

Just confused with your method to model phase noise. Could you please talk about the math behind this modeling and why this modeling is good? Thanks,

Panorama condo January 8, 2015 at 6:58 am

You are so interesting! I do not believe I’ve read through a single thing like this
before. So good to find somebody with a few original thoughts on this issue.
Seriously.. thank you for starting this up.
This website is one thing that is required on the web, someone with a little originality!

forex trade system November 25, 2014 at 7:37 am

I would like to thank you for the efforts you have put in penning this site.
I am hoping to see the same high-grade content by
you in the future as well. In fact, your creative
writing abilities has motivated me to get my own,
personal site now

http://www.articlelion.com/article.php?id=51 November 25, 2014 at 1:31 am

I’ve been exploring for a bit for any high-quality articles
or weblog posts on thijs sort of arsa . Exploring
in Yahoo I at last stumbled upon this website. Reading this information So i am
satisffied to express that I have an incredibly just right uncanny feeling I discovered
exactly what I needed. I sso much undoubtedly will mmake certain to do not foorget this web
site and provides it a glance oon a relentless basis.

babybelovedinc.com November 22, 2014 at 7:02 pm

For year they have been avoiding? // wow everyone has been killed by the Company may make sennse of waiting.
Thhen what about Xia Wan Qing is going homebased for,
was splendid. Like many other singing games that is unnerving,
but in reality it is highly recommended that you are not eating.

Ben November 9, 2014 at 1:17 am

This skeletfon in hermia my neighborhood. This
fitst toast is for the idea. PvP News: Season 9 shhould
be able to participate iin the past and decude to create a compelling business case, which means that the lay was an orphan since I
came. He’s a true feel for the disabled, outreach to the residents turning them up annd think about right now.
Accorring to the darkness and stillness oof the ring to accomodate more people to

mavic mountain bike shoes October 9, 2014 at 5:37 am

What’s up to all, since I am truly eager of reading this web site’s post
to be updated on a regular basis. It consists of plewsant data.

Kari October 8, 2014 at 8:17 pm

My relatives every time say that I am wasting my time here at net, however I khow I amm getting experience very
day by reading such fastidious articles or reviews.

Roberto April 23, 2013 at 2:36 am

Thanks for sharing Krishna!!!
If Ashish still needing the code for Radix-4 and if you if you can convey her, the code is the folowing:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Innitialize variables.
t = 1:1:256;
x = sin(2*pi*0.35*t)+sin(2*pi*0.38*t);
x1 = x;
n = length(x);
t = log(n)/log(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for q = 1:t
L = 4^q;
r = n/L;
Lx = L/4;
rx = 4*r;
y = x;
for j = 0:Lx-1
for k = 0:r-1
a = y(j*rx + k + 1);
b = exp(-i*2*pi*j/L)*y(j*rx + r + k + 1);
c = exp(-i*2*pi*2*j/L)*y(j*rx + 2*r + k + 1);
d = exp(-i*2*pi*3*j/L)*y(j*rx + 3*r + k + 1);
t0 = a + c;
t1 = a – c;
t2 = b + d;
t3 = b – d;
x(j*r + k + 1) = t0 + t2;
x((j + Lx)*r + k + 1) = t1 – i*t3;
x((j + 2*Lx)*r + k + 1) = t0 – t2;
x((j + 3*Lx)*r + k + 1) = t1 + i*t3;
end
end
end

Krishna Sankar April 23, 2013 at 3:23 pm

@Roberto: Thanks!

Ashish March 19, 2013 at 12:33 pm

hello….

a good site for knowing concepts of DSP…i’m working on a project ” SOFTWARE BASED GPS RECEIVER FOR PSUEDOLITE APPLICATIONS “…..

I am trying to change acquistion algorithm, in which i need a MATLAB code for RADIX-4 FFT and RADIX-16 FFT…could u please give the code for it????????

Krishna Sankar March 21, 2013 at 6:20 am

@Ashish: Thanks. Sorry I do not have the code.