Baseband-equivalent phase noise model

Markus Nentwig December 18, 20114 comments Coded in Matlab

Models a noisy local oscillator as found in wireless communications

Example output

On a logarithmic frequency axis:

On a linear frequency axis:


In wireless communications, signals are translated to and from radio frequency by multiplication with a local oscillator signal.

Ideally, the local oscillator is a pure sine wave. In reality, it includes a phase noise spectrum that will cause degradation of the wanted signal.
Unfortunately, modeling phase noise is not as straightforward as the usual ”additive white Gaussian noise”:
First, it requires information about the phase noise spectrum, which strongly depends on the application. For example, it is quite different for a base station, a mobile station or a measurement instrument due to different implementation constraints (cost, size, power consumption etc).
Then, given the spectrum, the generation of the phase noise signal as a ”streaming” signal source is not exactly straightforward, because it is difficult to design a reasonably-sized filter that approximates the phase noise spectrum accurately enough

This code snippet

The code below implements a phase noise source. Based on a phase noise spectrum (and throwing in some discrete spurious tones for added realism), a cyclic phase noise signal is generated and written to file. Full-length FFT filtering avoids the need to design a filter in order to approximate the PN-spectrum.

Cyclic phase noise

The generated phase noise can be used in a simulation or the like. If the data file is too short, it can be wrapped around as many times as needed. A cycle length of a few seconds should be sufficient in typical applications.

Very-low-frequency phase noise is usually not of interest, as it doesn't shift unwanted signals in frequency, and a receiver will remove it from the wanted signal via carrier recovery, frequency-offset correction etc. If necessary, a random-walk process for phase or the like could be added on the simulation side.

Using the phase noise model

Multiply a complex-valued sample stream with the generated phase noise signal.

If the argument includeCarrier is true, the product includes signal and phase noise.

Otherwise (includeCarrier == false), the product will be the isolated phase noise component, for example to visualize the spectrum of the introduced error.
Use with caution, as there may be correlation between PN product and the original signal if the instantaneous phase "wanders" too far from 0 degrees.


Be aware that phase noise products which fall beyond the Nyquist bandwidth will fold back and cause aliasing. Depending on what is being simulated (blockers etc), the simulation may require some ”frequency planning”.

The parameter fMax_Hz will set an upper limit on the generated phase noise (typically ~30 dB reduction, the remaining power beyond fMax_Hz is caused by higher-order sidebands from the nonlinear modulation).


The model can generate discrete spurious tones, if desired. Since a continuous-wave signal is infinitely narrow in frequency, the absolute power is specified in dBc (dB relative to carrier) instead of "dBc in 1 Hz" for the phase noise power spectral density.

Finding model parameters

One possibility is to study data sheets and publications on oscillator or synthesizer design. The included example data is by no means intended as "the" phase noise model. As the main purpose is to demonstrate phase noise, it may turn out to be too noisy for some applications.

A division by 2 ideally improves phase noise by 6 dB. For example, a phase noise spectrum published at 1800 MHz could be lowered by 6 dB as a first approximation, if using it at 900 MHz.

As a simple "sanity check", the total phase noise acting on the wanted signal needs to be low enough to reach the maximum SNR required by the radio standard (as a first guess, look up the SNR per symbol here for the modulation in question and an error rate of 1/1000).

Integrated RMS phase error

If needed, this code snippet can be used to calculate the integrated RMS phase error.

Running the code

Copy to "pn_generator.m" and run pn_generator.

The spectrum analyzer (description) can be optionally downloaded here. Copy it into the same directory as pn_generator.m.


% ****************************************************************
% baseband equivalent source of local oscillator with phase noise
% Markus Nentwig, 18.12.2011
% ****************************************************************
function pn_generator()
    close all;
    % ****************************************************************
    % PN generator configuration
    % ****************************************************************
    srcPar = struct();
    srcPar.n = 2 ^ 18; % generated number of output samples
    srcPar.rate_Hz = 7.68e6; % sampling rate
    srcPar.f_Hz = [0, 10e3, 1e6, 9e9]; % phase noise spectrum, frequencies
    srcPar.g_dBc1Hz = [-80, -80, -140, -140]; % phase noise spectrum, magnitude
    srcPar.spursF_Hz = [300e3, 400e3, 700e3]; % discrete spurs (set [] if not needed)
    srcPar.spursG_dBc = [-50, -55, -60]; % discrete spurs, power relative to carrier
    % ****************************************************************
    % run PN generator
    % ****************************************************************
    s = PN_src(srcPar);
    if false
        % ****************************************************************
        % export phase noise baseband-equivalent signal for use in simulator etc
        % ****************************************************************
        tmp = [real(s); imag(s)] .';
        save('phaseNoiseSample.dat', 'tmp', '-ascii');
    if exist('spectrumAnalyzer', 'file')
        % ****************************************************************
        % spectrum analyzer configuration
        % ****************************************************************
        SAparams = struct();
        SAparams.rate_Hz = srcPar.rate_Hz; % sampling rate of the input signal
        SAparams.pRef_W = 1e-3; % unity signal represents 0 dBm (1/1000 W)
        SAparams.pNom_dBm = 0; % show 0 dBm as 0 dB;
        SAparams.filter = 'brickwall';
        SAparams.RBW_window_Hz = 1000; % convolve power spectrum with a 1k filter
        SAparams.RBW_power_Hz = 1; % show power density as dBc in 1 Hz
        SAparams.noisefloor_dB = -250; % don't add artificial noise
        SAparams.logscale = true; % use logarithmic frequency axis
        % plot nominal spectrum
        figure(1); grid on; 
        h = semilogx(max(srcPar.f_Hz, 100) / 1e6, srcPar.g_dBc1Hz, 'k+-');
        set(h, 'lineWidth', 3);
        hold on;
        spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'fig', 1);
        ylabel('dBc in 1 Hz');
        legend('nominal PSD', 'output spectrum');
        title('check match with nominal PSD');
        spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'RBW_power_Hz', 'sine', 'fig', 2);
        title('check carrier level (0 dBc); check spurs level(-50/-55/-60 dBc)');
        ylabel('dBc for continuous-wave signal');
        spectrumAnalyzer('signal', s, SAparams, 'fig', 3, 'logscale', false);
        ylabel('dBc in 1 Hz');    

function pn_td = PN_src(varargin)
    def = {'includeCarrier', true, ...
           'spursF_Hz', [], ...
           'spursG_dBc', [], ...
           'fMax_Hz', []};
    p = vararginToStruct(def, varargin);
    % length of signal in the time domain (after ifft)
    len_s = p.n / p.rate_Hz

    % FFT bin frequency spacing
    deltaF_Hz = 1 / len_s
    % construct AWGN signal in the frequency domain     
    % a frequency domain bin value of n gives a time domain power of 1
    % for example ifft([4 0 0 0]) => 1 1 1 1 
    % each bin covers a frequency interval of deltaF_Hz
    mag = p.n;
    % scale "unity power in one bin" => "unity power per Hz":
    % multiply with sqrt(deltaF_Hz): 
    mag = mag * sqrt(deltaF_Hz);
    % Create noise according to mag in BOTH real- and imaginary value
    mag = mag * sqrt(2);

    % both real- and imaginary part contribute unity power => divide by sqrt(2)
    pn_fd = mag / sqrt(2) * (randn(1, p.n) + 1i * randn(1, p.n));
    % frequency vector corresponding to the FFT bins (0, positive, negative)
    fb_Hz = FFT_freqbase(p.n, deltaF_Hz);
    % interpolate phase noise spectrum on frequency vector
    % note: interpolate dB on logarithmic frequency axis
    H_dB = interp1(log(p.f_Hz+eps), p.g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');

    % dB => magnitude
    H = 10 .^ (H_dB / 20);
    %    H = 1e-6; % sanity check: enforce flat -120 dBc in 1 Hz
    % apply filter to noise spectrum
    pn_fd = pn_fd .* H;

    % set spurs
    for ix = 1:numel(p.spursF_Hz)
        fs = p.spursF_Hz(ix);
        u = abs(fb_Hz - fs);
        ix2 = find(u == min(u), 1);

        % random phase
        rot = exp(2i*pi*rand());
        % bin value of n: unity (carrier) power (0 dBc)
        % scale with sqrt(2) because imaginary part will be discarded
        % scale with sqrt(2) because the tone appears at both positive and negative frequencies
        smag = 2 * p.n * 10 ^ (p.spursG_dBc(ix) / 20);
        pn_fd(ix2) = smag * rot;
    % limit bandwidth (tool to avoid aliasing in an application
    % using the generated phase noise signal)
    if ~isempty(p.fMax_Hz)
        pn_fd(find(abs(fb_Hz) > p.fMax_Hz)) = 0;
    % convert to time domain
    pn_td = ifft(pn_fd);

    % discard imaginary part 
    pn_td = real(pn_td);

    % Now pn_td is a real-valued random signal with a power spectral density 
    % as specified in f_Hz / g_dBc1Hz.
    % phase-modulate to carrier
    % note: d/dx exp(x) = 1 near |x| = 1
    % in other words, the phase modulation has unity gain for small phase error
    pn_td = exp(i*pn_td);
    if ~p.includeCarrier
        % remove carrier
        % returns isolated phase noise component
        pn_td = pn_td - 1;

% returns a vector of frequencies corresponding to n FFT bins, when the
% frequency spacing between two adjacent bins is deltaF_Hz
function fb_Hz = FFT_freqbase(n, deltaF_Hz)
    fb_Hz = 0:(n - 1);
    fb_Hz = fb_Hz + floor(n / 2);
    fb_Hz = mod(fb_Hz, n);
    fb_Hz = fb_Hz - floor(n / 2);
    fb_Hz = fb_Hz * deltaF_Hz;

% *************************************************************
% helper function: Parse varargin argument list
% allows calling myFunc(A, A, A, ...)
% where A is
% - key (string), value (arbitrary) => result.key = value
% - a struct => fields of A are copied to result
% - a cell array => recursive handling using above rules
% *************************************************************
function r = vararginToStruct(varargin)
% note: use of varargin implicitly packs the caller's arguments into a cell array
% that is, calling vararginToStruct('hello') results in
%   varargin = {'hello'}
    r = flattenCellArray(varargin, struct());

function r = flattenCellArray(arr, r)
    ixMax = numel(arr);
    while ix <= ixMax
        e = arr{ix};
        if iscell(e)
            % cell array at 'key' position gets recursively flattened
            % becomes struct
            r = flattenCellArray(e, r);
        elseif ischar(e)
            % string => key.
            % The following entry is a value
            ix = ix + 1;
            v = arr{ix};
            % store key-value pair
            r.(e) = v;
        elseif isstruct(e)
            names = fieldnames(e);
            for ix2 = 1:numel(names)
                k = names{ix2};
                r.(k) = e.(k);
    end % while