plomb() - Signal Processing
example[pxx,f]
= plomb(x,t) returns
the Lomb-Scargle power spectral density (PSD) estimate, pxx,
of a signal, x, that is sampled at the instants
specified in t. t must increase
monotonically but need not be uniformly spaced. All elements of t must
be nonnegative. pxx is evaluated at the frequencies
returned in f.If x is a vector, it is treated
as a single channel.If x is a matrix, then plomb computes
the PSD independently for each column and returns it in the corresponding
column of pxx.x or t can contain NaNs
or NaTs. These values are treated as missing data
and excluded from the spectrum computation.
example[pxx,f]
= plomb(x,fs) treats
the case where the signal is sampled uniformly, at a rate fs,
but has missing samples. Specify the missing data using NaNs.
[pxx,f]
= plomb(___,fmax) estimates
the PSD up to a maximum frequency, fmax, using
any of the input arguments from previous syntaxes. If the signal is
sampled at N non-NaN instants,
and Δt is
the time difference between the first and the last of them, then pxx is
returned at round(fmax / fmin) points,
where fmin = 1/(4 × N × ts) is
the smallest frequency at which pxx is computed
and the average sample time is ts = Δt/(N – 1). fmax defaults
to 1/(2 × ts),
which for uniformly sampled signals corresponds to the Nyquist frequency.
example[pxx,f]
= plomb(___,fmax,ofac) specifies
an integer oversampling factor, ofac. The use
of ofac to interpolate or smooth a spectrum resembles
the zero-padding technique for FFT-based methods. pxx is
again returned at round(fmax/fmin)
frequency points, but the minimum frequency considered in this case
is 1/(ofac × N × ts). ofac defaults
to 4.
example[pxx,fvec]
= plomb(___,fvec) estimates
the PSD of x at the frequencies specified in fvec. fvec must
have at least two elements. The second output argument is the same
as the input fvec.You cannot specify a maximum frequency or an oversampling factor
if you use this syntax.
example[___] = plomb(___,spectrumtype) specifies
the normalization of the periodogram.Set spectrumtype to 'psd',
or leave it unspecified, to obtain pxx as a power
spectral density.Set spectrumtype to 'power' to
get the power spectrum of the input signal.Set spectrumtype to 'normalized' to
get the standard Lomb-Scargle periodogram, which is scaled by two
times the variance of x.
example[___,pth] = plomb(___,'Pd',pdvec) returns
the power-level threshold, pth, such that a peak
with a value larger than pth has a probability pdvec of
being a true signal peak and not the result of random fluctuations. pdvec can
be a vector. Every element of pdvec must be greater
than 0 and smaller than 1. Each row of pth corresponds
to an element of pdvec. pth has
the same number of channels as x. This option
is not available if you specify the output frequencies in fvec.
example[pxx,w]
= plomb(x) returns the PSD estimate
of x evaluated at a set of evenly spaced normalized
frequencies, w, spanning the Nyquist interval.
Use NaNs to specify missing samples. All of the
above options are available for normalized frequencies. To access
them, specify an empty array as the second input.
exampleplomb(___) with no output arguments
plots the Lomb-Scargle periodogram PSD estimate in the current figure
window.
Syntax
<![CDATA[[pxx,f]
= plomb(x,t) example[pxx,f]
= plomb(x,fs) example[pxx,f]
= plomb(___,fmax)[pxx,f]
= plomb(___,fmax,ofac) example[pxx,fvec]
= plomb(___,fvec) example[___] = plomb(___,spectrumtype) example[___,pth] = plomb(___,'Pd',pdvec) example[pxx,w]
= plomb(x) exampleplomb(___) example]]>
Example
<![CDATA[Irregularly Sampled Signal and Signal with Missing SamplesOpen This Example
The Lomb-Scargle method can handle signals that have been sampled unevenly or have missing samples.
Generate a two-channel sinusoidal signal sampled at 1 kHz for about 0.5 s. The sinusoid frequencies are 175 Hz and 400 Hz. Embed the signal in white noise with variance
.Fs = 1000;
f0 = 175;
f1 = 400;
t = 0:1/Fs:0.5;
wgn = randn(length(t),2)/2;
sigOrig = sin(2*pi*[f0;f1]*t)' + wgn;
Compute and plot the periodogram of the signal. Use periodogram with the default settings.periodogram(sigOrig,[],[],Fs)
axisLim = axis;
title('Periodogram')
Use plomb with the default settings to estimate and plot the PSD of the signal. Use the axis limits from the previous plot.plomb(sigOrig,t)
axis(axisLim)
title('Lomb-Scargle')
Suppose the signal is missing 10% of the original samples. Place NaN's in random locations to simulate the missing data points. Use plomb to estimate and plot the PSD of the signal with missing samples.sinMiss = sigOrig;
misfrac = 0.1;
nTime = length(t)*2;
sinMiss(randperm(nTime,round(nTime*misfrac))) = NaN;
plomb(sinMiss,t)
axis(axisLim)
title('Missing Samples')
Sample the original signal, but make the sampling nonuniform by adding jitter (uncertainty) to the time measurements. The first instant continues to be at zero. Use plomb to estimate and plot the PSD of the nonuniformly sampled signal.tirr = t + (1/2-rand(size(t)))/Fs/2;
tirr(1) = 0;
sinIrreg = sin(2*pi*[f0;f1]*tirr)' + wgn;
plomb(sinIrreg,tirr)
axis(axisLim)
title('Nonuniform Sampling')
Periodogram of Data Set with Missing SamplesOpen This Example
Galileo Galilei observed the motion of Jupiter's four largest satellites during the winter of 1610. When the weather allowed, Galileo recorded the satellites' locations. Use his observations to estimate the orbital period of one of the satellites, Callisto.
Callisto's angular position is measured in minutes of arc. Missing data due to cloudy conditions are specified using NaNs. The first observation is dated January 15. Generate a datetime array of observation times.yg = [10.5 NaN 11.5 10.5 NaN NaN NaN -5.5 -10.0 -12.0 -11.5 -12.0 -7.5 ...
NaN NaN NaN NaN 8.5 12.5 12.5 10.5 NaN NaN NaN -6.0 -11.5 -12.5 ...
-12.5 -10.5 -6.5 NaN 2.0 8.5 10.5 NaN 13.5 NaN 10.5 NaN NaN NaN ...
-8.5 -10.5 -10.5 -10.0 -8.0]';
obsv = datetime(1610,1,14+(1:length(yg)));
plot(yg,'o')
ax = gca;
nights = [1 18 32 46];
ax.XTick = nights;
ax.XTickLabel = char(obsv(nights));
Estimate the power spectrum of the data using plomb. Specify an oversampling factor of 10. Express the resulting frequencies in inverse days.[pxx,f] = plomb(yg,obsv,[],10,'power');
f = f*86400;
Use findpeaks to determine the location of the only prominent peak of the spectrum. Plot the power spectrum and show the peak.[pk,f0] = findpeaks(pxx,f,'MinPeakHeight',10);
plot(f,pxx,f0,pk,'o')
xlabel('Frequency (day^{-1})')
Determine Callisto's orbital period (in days) as the inverse of the frequency of maximum energy. The result differs by less than 1% from the value published by NASA.Period = 1/f0
NASA = 16.6890184;
PercentDiscrep = (Period-NASA)/NASA*100
Period =
16.6454
PercentDiscrep =
-0.2613
Periodogram of Data Set with Irregular SamplingOpen This Example
Galileo Galilei discovered Jupiter's four largest satellites in January of 1610 and recorded their locations every clear night until March of that year. Use Galileo's data to find the orbital period of Callisto, the outermost of the four satellites.
Galileo's observations of Callisto's angular position are in minutes of arc. There are several gaps due to cloudy conditions.t = [0 2 3 7 8 9 10 11 12 17 18 19 20 24 25 26 27 28 29 31 32 33 35 37 ...
41 42 43 44 45]';
yg = [10.5 11.5 10.5 -5.5 -10.0 -12.0 -11.5 -12.0 -7.5 8.5 12.5 12.5 ...
10.5 -6.0 -11.5 -12.5 -12.5 -10.5 -6.5 2.0 8.5 10.5 13.5 10.5 -8.5 ...
-10.5 -10.5 -10.0 -8.0]';
plot(t,yg,'o')
xlabel('Day')
Use plomb to compute the periodogram of the data. Estimate the power spectrum up to a frequency of
. Specify an oversampling factor of 10. Choose the standard Lomb-Scargle normalization.[pxx,f] = plomb(yg,t,0.5,10,'normalized');
The periodogram has one clear maximum. Name the peak frequency
. Plot the periodogram and annotate the peak.[pmax,lmax] = max(pxx);
f0 = f(lmax);
plot(f,pxx,f0,pmax,'o')
xlabel('Frequency (day^{-1})')
Use linear least squares to fit to the data a function of the form
The fitting parameters are the amplitudes A, B, and C.ft = 2*pi*f0*t;
ABC = [ones(size(ft)) cos(ft) sin(ft)] \ yg
ABC =
0.4243
10.4444
6.6137
Use the fitting parameters to construct the fitting function on a 200-point interval. Plot the data and overlay the fit.x = linspace(t(1),t(end),200)';
fx = 2*pi*f0*x;
y = [ones(size(fx)) cos(fx) sin(fx)] * ABC;
plot(t,yg,'o',x,y)
xlabel('Day')
Irregular Sampling and AliasingOpen This Example
Sample a 0.8 Hz sinusoid at 1 Hz for 100 s. Embed the sinusoid in white noise with a variance of 1/100. Reset the random number generator for repeatable results.
f0 = 0.8;
rng default
wgn = randn(1,100)/10;
ts = 1:100;
s = sin(2*pi*f0*ts) + wgn;
Compute and plot the power spectral density estimate up to the sample rate. Specify an oversampling factor of 10.plomb(s,ts,1,10)
The aliasing results from the fact that the frequency of the sinusoid is greater than the Nyquist frequency.Repeat the calculation, but now sample the sinusoid at random times. Include frequencies up to 1 Hz. Specify an oversampling factor of 2. Plot the PSD.tn = sort(100*rand(1,100));
n = sin(2*pi*f0*tn) + wgn;
ofac = 2;
plomb(n,tn,1,ofac)
The aliasing disappears. The irregular sampling increases the effective sample rate by shrinking some time intervals.Zoom in on the frequencies around 0.8 Hz. Use a fine grid with a spacing of 0.001 Hz. You cannot specify an oversampling factor or a maximum frequency in this case.df = 0.001;
fvec = 0.7:df:0.9;
hold on
plomb(n,tn,fvec)
legend('ofac = 2','df = 0.001')
Exponential DistributionOpen This ExampleGenerate
samples of white noise with variance
, 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
. 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
. To normalize the histograms, recall that the total number of periodogram samples is
. 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
. Specify the sinusoid amplitude,
, using the relation
. 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
. 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
False-Alarm ProbabilitiesOpen This ExampleGenerate 100 samples of a sinusoidal signal at a sample rate of 1 Hz. Specify an amplitude of 0.75 and a frequency of
. Embed the signal in white noise of variance 0.902. Reset the random number generator for repeatable results.rng default
X0 = 0.75;
f0 = 0.6;
vr = 0.902;
Nsamp = 100;
t = 1:Nsamp;
X = X0*cos(f0*(1:Nsamp))+randn(1,Nsamp)*sqrt(vr);
Discard 10% of the samples at random. Plot the signal.X(randperm(Nsamp,Nsamp/10)) = NaN;
plot(t,X,'o')
Compute and plot the normalized power spectrum. Annotate the levels that correspond to false-alarm probabilities of 50%, 10%, 1%, and 0.01%. If you generate many 90-sample white noise signals with variance 0.902, then half of them have one or more peaks higher than the 50% line, 10% have one or more peaks higher than the 10% line, and so on.Pfa = [50 10 1 0.01]/100;
Pd = 1-Pfa;
[pxx,f,pth] = plomb(X,1:Nsamp,'normalized','Pd',Pd);
plot(f,pxx,f,pth*ones(size(f')))
xlabel('f')
text(0.3*[1 1 1 1],pth-.5,[repmat('P_{fa} = ',[4 1]) num2str(Pfa')])
In this case, the peak is high enough that only about 0.01% of the possible signals can attain it.Use plomb with no output arguments to repeat the calculation. The plot is now logarithmic, and the levels are drawn in terms of detection probabilities.plomb(X,1:Nsamp,'normalized','Pd',Pd)
Lomb-Scargle Periodogram of Noisy SinusoidsOpen This Example
When given a data vector as the only input, plomb estimates the power spectral density using normalized frequencies.
Generate 128 samples of a sinusoid of normalized frequency
rad/sample embedded in white Gaussian noise of variance 1/100.t = (0:127)';
x = sin(2*pi*t/4);
x = x + randn(size(x))/10;
Estimate the PSD using the Lomb-Scargle procedure. Repeat the calculation with periodogram.[p,f] = plomb(x);
[pper,fper] = periodogram(x);
Plot the PSD estimates in decibels. Verify that the results are equivalent.plot(f/pi,pow2db(p))
hold on
plot(fper/pi,pow2db(pper))
axis([0 1 -40 20])
xlabel('\omega / \pi')
ylabel('PSD')
legend('plomb','periodogram')
Estimate the Lomb-Scargle PSD of a three-channel signal composed of sinusoids. Specify the frequencies as
rad/sample,
rad/sample, and
rad/sample. Add white Gaussian noise of variance 1/100. Use plomb with no output arguments to compute and plot the PSD estimate in decibels.x3 = [sin(2*pi*t/3) sin(2*pi*t/4) sin(2*pi*t/5)];
x3 = x3 + randn(size(x3))/10;
figure
plomb(x3)
Compute the PSD estimate again, but now remove 25% of the data at random.x3(randperm(numel(x3),0.25*numel(x3))) = NaN;
plomb(x3)
Power Spectral Density of Signal with Missing SamplesOpen This Example
If you do not have a time vector, use NaN's to specify missing samples in a signal.
Generate 1024 samples of a sinusoid of normalized frequency
rad/sample embedded in white noise of variance 1/100. Estimate the power spectral density using the Lomb-Scargle procedure. Use plomb with no output arguments to plot the estimate.t = (0:1023)';
x = sin(2*pi*t/5);
x = x + randn(size(x))/10;
[pxx,f] = plomb(x);
plomb(x)
Remove every other sample by assigning NaN's. Use plomb to compute and plot the PSD. The periodogram peaks at the same frequency because the time axis is unchanged.xnew = x;
xnew(2:2:end) = NaN;
plomb(xnew)
Remove every other sample by downsampling. The function now estimates the periodicity at twice the original frequency. This is probably not the result you want.xdec = x(1:2:end);
plomb(xdec)
Related ExamplesDetect Periodicity in a Signal with Missing SamplesSpectral Analysis of Nonuniformly Sampled Signals]]>
Output / Return Value
Limitations
Alternatives / See Also
Reference