# How to accurately compute the Winger-Ville Distribution of an exponential

I am using MATLAB for this question so hopefully you can help me out.

I am trying to compute the Wigner-Ville distribution (WVD) of a sinusoidal signal defined as $$$$x(t) = e^{-i\omega_0 t} e^{i\phi(t)},$$$$ where $$\omega_0$$ is the central frequency (rad/s) and $$\phi(t)$$ is a Brownian phase noise function, defined as a random walk $$\phi(t) = \sqrt{\Delta\Omega}W(t)$$ and $$W(t)$$ is a Wiener process $$W(t)\sim \sqrt{|t|} N(0,1)$$. The parameter $$\Delta\Omega$$ represents the linewidth of my laser. The idea is that I can compute the auto-correlation function of $$x(t)$$ $$$$r_x(t+\tau/2,t-\tau/2) = E\left[x(t+\tau/2)x^*(t-\tau/2)\right] = e^{-i\omega_0 \tau} e^{-\Delta\Omega |\tau|/2 },$$$$ where $$E[\cdot]$$ is the expectation operator or ensemble average. Once I compute this, I can then compute the average Wigner-Ville Distribution using $$$$E[W(t,\omega)] = \int_{-\infty}^{\infty} r_x(t+\tau/2,t-\tau/2)e^{i\omega \tau}\ d\tau,$$$$ Or essentially as the "Fourier transform" of the function $$r_x$$. This gives me the analytic formula $$$$E[W(t,\omega)] = 2\frac{\Delta\Omega/2}{(\Delta\Omega/2)^2 + (\omega-\omega_0)^2}.$$$$

I can then compute the power spectral density (PSD) of $$x(t)$$ by integrating $$E[W(t,\omega)]$$ over $$t$$.

The problem now comes in when I try to replicate this result numerically, that is simulating $$N$$ random walks for $$\phi(t)$$, and creating a matrix $$X$$ which contains the signal on a time domain $$t\in [0,T]$$ for each of those random walks. I haven't been able to find any MATLAB functions that can compute the autocorrelation function $$r_x$$ for me, so instead I am trying to use the built-in function wvd() (Wigner-Ville distribution) directly on each of my $$x$$ simulations contained in the matrix $$X$$. The problem is that the result I get from the wvd() function has a different line-width in frequency space than the analytic (in fact the line-width of my function should be exactly $$\Delta\Omega$$, or simply $$\Delta\nu = \Delta\Omega/2\pi$$ when plotted on $$f$$ (Hz)), but that is not the case here (see plot below):

Does anyone know why this happens? Here is my MATLAB code:


%% --- INPUTS --- %%
% User-defined parameters
T = 20e-9;    % Time duraction of signal (s)
f0 = 200e12;  % Laser center frequency (Hz)
A = 1;        % Laser amplitude
FWHM = T/5;   % Full-Width at Half Maximum of Gaussian pulse
Nt = 1001;    % Sample points in time
Nf = Nt;      % Sample points in frequency
N_MC = 20;  % Monte Carlo Simulations
dNu = 30e6;   % Laser linewidth (Hz)

% Calculated parameters
w0 = 2*pi*f0;
dOm = 2*pi*dNu;
t = linspace(0,T,Nt);
dt = t(2)-t(1);
a = 4*log(2)/FWHM^2;
t0 = T/3;
df = 1/dt/Nf;
% f = f0 + df.*linspace(-Nf/2,Nf/2,Nf); % shift array to center f
% w = 2*pi*f;

% Signal
phi = @(t) sqrt(dOm)*cumsum(sqrt(dt)*randn(size(t))); % phase noise
x = @(t) exp(-1i*w0*t).*exp(1i*phi(t));

%% --- Analytic --- %%
% Signal properties
r_x = @(tau) A*exp(-1i*w0*tau).*exp(-dOm/2*abs(tau));
PSD_analytic = @(f) A*dOm./((dOm/2)^2 + 4*pi^2*(f-f0).^2);

% Phase properties
phi_var_analytic = @(t,dOm) dOm*abs(t);
phi_mgf_analytic = @(t,dOm) exp(-1/2*dOm*abs(t));

%% --- Numeric --- %%
% Phase properties
Tmatrix = columnVector(t).*ones(Nt,N_MC);
phi0 = phi(Tmatrix);
phi_var = mean( phi0.^2, 2) - mean(phi0,2).^2;
phi_mgf = mean( exp(1i*phi0), 2);

% Signal properties
x_samples = x(Tmatrix);
[f,Spectrum,PSD,comp_time] = WVS(t,x_samples,f0);
disp(['Computation time = ' num2str(comp_time/60) ' min']);

%% --- Plots --- %%
close all;
FS = '\fontname{Palatino} ';
fontS = 18;
colors = {'red','black'};
LW = 2;
factor0 = 50;
f00 = f0/1e12;
fmin0 = min(f)/1e12;
fmax0 = max(f)/1e12;
dfmax = (fmax0-fmin0)/2;
fmin = f00 - dfmax/factor0;
fmax = f00 + dfmax/factor0;

% ---> Phase Plots
figure('units','normalized','outerposition',[0 0 1 1]);
set(gcf,'color','w');

subplot(221);
plot(t*1e9,phi_var_analytic(t,dOm),'color',colors{1},'LineWidth',LW);
hold on;
plot(t*1e9,phi_var,'marker','.','color',colors{2},...
'LineStyle','none','LineWidth',LW+1);
hold off;
legend boxoff; legend([FS 'Analytic'], [FS 'Numeric']);
axis tight;
xlabel([FS 't (ns)']);
title([FS 'Variance in the phase']);
set(gca,'FontSize',fontS);

subplot(223);
plot(t*1e9,phi_mgf_analytic(t,dOm),'color',colors{1},'LineWidth',LW);
hold on;
plot(t*1e9,phi_mgf,'marker','.','color',colors{2},...
'LineStyle','none','LineWidth',LW+1); hold off;
legend boxoff; legend([FS 'Analytic'], [FS 'Numeric']);
axis tight;
xlabel([FS 't (ns)']);
ylabel([FS '\langle' 'e^{i\phi(t)}\rangle']);
title([FS 'MGF of the phase']);
set(gca,'FontSize',fontS);

% ---> Spectral Plots (Normalized)
subplot(222);
plot(f/1e12,PSD_analytic(f)./max(PSD_analytic(f)),'color',colors{1},...
'LineWidth',LW); hold on;
plot(f/1e12,PSD./max(PSD),'marker','.','color',colors{2},...
'LineStyle','none','LineWidth',LW+1); hold off;
legend boxoff; legend([FS 'Analytic'], [FS 'Numeric']);
axis tight;
xlabel([FS 'f (THz)']);
ylabel([FS '\langle' '|X(f)|^2' '\rangle (normalized units)']);
title([FS 'PSD (linear scale)']);
set(gca,'FontSize',fontS);

subplot(224);
plot(f/1e12,10*log10(PSD_analytic(f)./max(PSD_analytic(f))),'color',...
colors{1},'LineWidth',LW); hold on;
plot(f/1e12,10*log10(PSD./max(PSD)),'marker','.','color',colors{2},...
'LineStyle','none','LineWidth',LW+1); hold off;
legend boxoff; legend([FS 'Analytic'], [FS 'Numeric']);
axis tight;
xlabel([FS 'f (THz)']);
ylabel([FS '\langle' '|X(f)|^2' '\rangle (normalized units)']);
title([FS 'PSD (dB scale)']);
set(gca,'FontSize',fontS);

%% --- Functions --- %%
function [freq,Spectrum,PSD,comp_time] = WVS(t_array,x_array,f0)
% Ensembled Averaged Wigner-Ville Distribution
% x_array: each column is a single ensemble of the random signal x(t)
N = length(t_array);
dt = abs(t_array(2)-t_array(1));
fs = 1/dt;
[~,N_MC] = size(x_array);
W = 0;
tic;
for n = 1:N_MC % Do ensemble averaging over multiple Monte Carlo samples
[W0,f,~] = wvd(ifftshift(x_array(:,n)),fs);
Nf = length(f);
df = 1/dt/Nf;
W = W + dt*W0;
if rem(n,10)==0
disp(n);
end
end
freq = f0 + f - max(f)/2;
Spectrum = abs(fftshift(W/N_MC));
PSD = dt*trapz(Spectrum,2);
comp_time = toc;
end

function output = columnVector(input_vector)
[sx,sy] = size(input_vector);
if sx < sy
output = input_vector.';
else
output = input_vector;
end
end


PS: I know that my function $$x(t)$$ in this case is stationary (since $$r_x$$ depends only on $$\tau=t-t'$$ and not on two separate times $$(t,t')$$; and so I could calculate the PSD using the Wiener-Khinchin theorem, but I want to use the WVD instead so that my code can be generalized for non-stationary functions $$x(t)$$ (e.g. such as Gaussian pulses)

EDIT: I have changed my definition of phi(t) in the code so that it now is computed as a random walk with $$$$\phi(t+\Delta t) = \phi(t) + \sqrt{\Delta\Omega}\sqrt{\Delta t}N(0,1),$$$$ but this still does not solve the problem of the numerical PSD being a lot wider than the analytic PSD, and I can't figure out why

Your formula for $$\displaystyle W(t)$$ does not represent a Wiener process.

By definition, the Wiener process $${\displaystyle W}$$ has Gaussian increments: $${\displaystyle W_{t+u}-W_{t}}$$ is normally distributed with mean $${\displaystyle 0}$$ and variance $${\displaystyle u}$$, $${\displaystyle W_{t+u}-W_{t}\sim {\mathcal {N}}(0,u)}$$. Just because $$\displaystyle W_t^2 − t$$ is a martingale, does not result in your formula $$\displaystyle W(t)\sim \sqrt{|t|} N(0,1)$$.

You can easily see it, comparing two plots:

• fill in an array $$GN[1:N]$$ with a realization of Gaussian process.
• draw a plot for $$sqrt(k)·GN[k]$$
• construct an array of partial sums, $$W[k] = \sum_{i=1}^{i=k}GN[k]$$, which is a realization of the Wiener process by definition
• draw a plot for $$W[k]$$.

The jumps in the sequence $$sqrt(k)·GN[k]$$ increase with $$k$$ ($$k$$ represents discrete time); the jumps in the partial sums $$W[k]$$ vary around a Gaussian process variance and do not increase with $$k$$.

Notice also that there exists a representation of a Wiener process via i.i.d. Gaussian variable $${\displaystyle ξ_{n}}$$ in terms of a random Fourier series: $${\displaystyle W_{t}={\sqrt {2}}\sum _{n=1}^{\infty }ξ_{n}{\frac {\sin \left(\left(n-{\frac {1}{2}}\right)\pi t\right)}{\left(n-{\frac {1}{2}}\right)\pi }}}$$ on $${\displaystyle [0,1]}$$. Comparing it to your formula for $$W(t)$$ is not a trivial task, but it indicates that you cannot just multiply the Gaussian noise realization with the time square root and have the Wiener process realization.

• Yes this is correct, I had defined phi(t) wrong in my code. I have adjusted it so that it now is computed as a random walk via the equation: phi(t+dt) = phi(t) + sqrt(DeltaOmega)*sqrt(dt)*randn. However, that sill doesn't solve my problem. The main issue is that the frequency spectrum that results from using the Wigner-Ville Distribution function in matlab (wvd()) gives me a spectrum which is a lot wider than it should be, compared to the analytic result and I can't figure out if my implementation of wvd() is wrong – Oscar Andres Nieves Feb 13 at 2:32
• Not "a lot wider than it should be" (as you write): in fact, your corrected numerical plot coincides with your analytic plot in the region where the normalized PSD is greater than 0.2. As for the shoulders of the distribution, have you ever heard of an ensemble averaging? You fail to approach the analytic plot with your simulation because the naive power spectral density obtained from the signal's Fourier transform is a biased estimate of the true spectral content, when done over a single realization of the random process, as is the case in your simulation. – V.V.T Feb 13 at 6:26
• To overcome this trouble, you can average over an ensemble of simulations, or learn and apply the multitaper method (en.wikipedia.org/wiki/Multitaper). – V.V.T Feb 13 at 6:27