%% Exponential Distribution
% Generate N=1024 samples of white noise with variance \sigma=1, given a sample rate of 1 Hz. Compute the power spectrum of the white noise. Choose the Lomb-Scargle normalization and specify an oversampling factor {\tt ofac}=15. Reset the random number generator for repeatable results.

rng default

N = 1024;
t = (1:N)';
wgn = randn(N,1);

ofac = 15;
[pwgn,f] = plomb(wgn,t,[],ofac,'normalized');

%
% Verify that the Lomb-Scargle power spectrum estimate of white noise has an exponential distribution with unit mean. Plot a histogram of the values of pwgn and a histogram of a set of exponentially distributed random numbers generated using the distribution function f(z|1)=e^{-z}. To normalize the histograms, recall that the total number of periodogram samples is N\times{\tt ofac}/2. Specify a bin width of 0.25. Overlay a plot of the theoretical distribution function.

dx = 0.25;
br = 0:dx:5;

Nf = N*ofac/2;

hpwgn = histcounts(pwgn,br)';

hRand = histcounts(-log(rand(Nf,1)),br)';

bend = br(1:end-1);

bar(bend,[hpwgn hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('wgn','Empirical pdf','Theoretical pdf')
hold off

%
% Embed in the noise a sinusoidal signal of frequency 0.1 Hz. Use a signal-to-noise ratio of \xi=0.01. Specify the sinusoid amplitude, x_0, using the relation x_0=\sigma\sqrt{2\xi}. Compute the power spectrum of the signal and plot its histogram alongside the empirical and theoretical distribution functions.

SNR = 0.01;
x0 = sqrt(2*SNR);
sigsmall = wgn + x0*sin(2*pi/10*t);

[psigsmall,f] = plomb(sigsmall,t,[],ofac,'normalized');

hpsigsmall = histcounts(psigsmall,br)';

bar(bend,[hpsigsmall hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('sigsmall','Empirical pdf','Theoretical pdf')
hold off

%
% Repeat the calculation using \xi=1. The distribution now differs noticeably.

SNR = 1;
x0 = sqrt(2*SNR);
siglarge = wgn + x0*sin(2*pi/10*t);

[psiglarge,f] = plomb(siglarge,t,[],ofac,'normalized');

hpsiglarge = histcounts(psiglarge,br)';

bar(bend,[hpsiglarge hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('siglarge','Empirical pdf','Theoretical pdf')
hold off

% Copyright 2015 The MathWorks, Inc.