%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All rights reserved by Krishna Pillai, http://www.dsplog.com
% The file may not be re-distributed without explicit authorization
% from Krishna Pillai.
% Checked for proper operation with Octave Version 3.0.0
% Author        : Krishna Pillai
% Email         : krishna@dsplog.com
% Version       : 1.0
% Date          : 22nd  June 2012
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Script for computing the integrated phase jitter 
% from the power spectral density (dBc/Hz) of phase noise
% 

clear all; close all;
fc_hz    = 10e6;
psd_f_hz           = [1e3  1e4  1e5  1e6  10e6];
psd_val_dbc_per_hz = [-65  -95 -115 -125  -125];

% computations
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 index with slope == -10
idx = find(slope==-10);
integral(idx) = constant(idx).*(log(psd_f_hz(idx+1)) - log(psd_f_hz(idx)));

% results
rms_jitter_radians  = sqrt(2*integral);
rms_jitter_seconds  = rms_jitter_radians/(2*pi*fc_hz)

integrated_jitter_radians = sqrt(2*sum(integral))
integrated_jitter_seconds = integrated_jitter_radians/(2*pi*fc_hz)
	
% plotting
semilogx(psd_f_hz,psd_val_dbc_per_hz,'bs-','Linewidth',2);
axis([1e3 1e7 -130 -60]); grid on
xlabel('freq, Hz'); ylabel('psd, dbc/hz');
title_str = sprintf('example phase noise profile (fc = %d MHz, rms jiiter = %4.6f ps)', ...
		fc_hz/1e6, integrated_jitter_seconds*1e12);
title(title_str);







