Wolfram Language & System Documentation Center

PowerSpectralDensity

PowerSpectralDensity[data,ω]

estimates the power spectral density for data.

PowerSpectralDensity[data,ω,sspec]

estimates the power spectral density for data with smoothing specification sspec.

PowerSpectralDensity[tproc,ω]

represents the power spectral density of a time series process tproc.

Details and Options

PowerSpectralDensity is also known as the energy spectral density.
PowerSpectralDensity[tproc,ω] is defined for weakly stationary time series processes as , where denotes CovarianceFunction[proc,h].
The following smoothing specifications sspec can be given:
c use c as a cutoff

w use a window function w

{c,w} use both a cutoff and a window function
For a window function w and positive integer c, PowerSpectralDensity[data,ω,{c,w}] is computed as , where is defined as CovarianceFunction[data,h].
By default, the cutoff c is chosen to be , where is the length of data, and the window function is DirichletWindow.
A window function is an even function such that , , for , including standard windows such as HammingWindow, ParzenWindow, etc.
A window function can be given as a list of values {w₀,…}, where , and it will be applied symmetrically in the vector case.
PowerSpectralDensity takes the FourierParameters option. Common settings for FourierParameters include:
{1,1} default setting

{-1,1} often used for time series

{a,b} general setting

Examples

open all close all

Basic Examples (3)

Estimate the power spectral density for some data:

Wolfram Language code: PowerSpectralDensity[Range[10], ω]

Calculate the power spectral density for a univariate time series:

Wolfram Language code: PowerSpectralDensity[ARProcess[{a}, σ^2], ω]

The sample power spectral density for a random sample from autoregressive time series:

Wolfram Language code: data = RandomFunction[ARProcess[{.2}, .1], {1, 1000}];

Calculate power spectral density with cutoff:

Wolfram Language code:

cutoffs = {5, 10, 17};
spec = Table[PowerSpectralDensity[data, ω, i], {i, cutoffs}];

Wolfram Language code: Plot[Evaluate@spec, {ω, -π, π}, PlotLegends -> cutoffs]

Scope (14)

Empirical Estimates (4)

Estimate the power spectral density for a univariate time series:

Wolfram Language code:

data = TemporalData[TimeSeries, {{{-0.0614065187043199, 2.6306604608936857, 3.9465133828074777, 
    3.745705410210661, 3.049768958525431, 1.4003183889805526, -0.755044374280679, 
    -2.0834316975526366, -1.3583280202011117, 0.14932145182161072, -0.8588219015885377}}, 
  {{0, 10, 1}}, 1, {"Continuous", 1}, {"Discrete", 1}, 1, 
  {ValueDimensions -> 1, ResamplingMethod -> None}}, False, 10.1];

Wolfram Language code: PowerSpectralDensity[data, ω]

Wolfram Language code: Plot[%, {ω, -π, π}, Filling -> Axis, PlotRange -> All]

Power spectral density for a vector time series:

Wolfram Language code:

data = TemporalData[TimeSeries, {{{{-0.8805412843977076, 0.5330873177260956}, 
    {-1.091904325250007, 2.8824771870568764}, {0.4712280848587277, -2.563788099997642}, 
    {-1.0458566604266817, 1.2863067750997783}, {0.15937341588902115, 0.0582751907559111 ... 0441, -1.4624411116452787}, 
    {-1.2038000803608786, 0.6266434204415561}, {-1.5532848855025663, -1.4784686973103551}}}, 
  {{0, 10, 1}}, 1, {"Continuous", 1}, {"Discrete", 1}, 2, 
  {ValueDimensions -> 2, ResamplingMethod -> None}}, False, 10.1];

Wolfram Language code: psd = PowerSpectralDensity[data, ω]//Simplify

Power spectral density for each component:

Wolfram Language code: Table[Plot[psd[[i, i]], {ω, -π, π}], {i, 1, 2}]

Cross power spectral density between components:

Wolfram Language code:

Table[ParametricPlot[{Re[#], Im[#]}&@psd[[i, Mod[i, 2] + 1]], {ω, -π, π}, ColorFunction -> Function[{x, y, ω}, (ColorData["Rainbow"][ω])]], {i, 1, 2}]

Estimate the power spectral density for an ensemble of paths:

Wolfram Language code: data = RandomFunction[ARProcess[{.8}, 1], {0, 100}, 5];

Wolfram Language code: psd = PowerSpectralDensity[data, ω];

Wolfram Language code: Plot[psd, {ω, -π, π}, PlotRange -> All, Filling -> Axis]

Compare empirical and theoretical power spectral densities functions:

Wolfram Language code: proc = MAProcess[{.4, .3, .5, .6, .3}, 1.];

Wolfram Language code: data = RandomFunction[proc, {0, 10 ^ 2}];

Wolfram Language code:

Plot[#, {ω, -π, π}, Filling -> Axis, PlotRange -> All]& /@ {PowerSpectralDensity[proc, ω], PowerSpectralDensity[data, ω]}

Smoothing (5)

Obtain a smoothed estimate using a cutoff at 5:

Wolfram Language code: data = Range[10];

Wolfram Language code: sspec = PowerSpectralDensity[data, ω, 5]

Compare the smoothed spectrum to the original:

Wolfram Language code: Plot[{Evaluate@PowerSpectralDensity[data, ω], sspec}, {ω, -π, π}, Filling -> Axis, PlotRange -> All]

Compute the power spectral density using a NuttallWindow:

Wolfram Language code: data = Range[10];

Wolfram Language code: sspec = PowerSpectralDensity[data, ω, NuttallWindow]

Compare the smoothed spectrum to the original:

Wolfram Language code: Plot[{Evaluate@PowerSpectralDensity[data, ω], sspec}, {ω, -π, π}, Filling -> Axis, PlotRange -> All]

Define a window using a pure function:

Wolfram Language code: data = Range[10];

Wolfram Language code: sspec = PowerSpectralDensity[data, ω, HannPoissonWindow[#, (1/2)]&]

Compare the smoothed spectrum to the original:

Wolfram Language code: Plot[{Evaluate@PowerSpectralDensity[data, ω], sspec}, {ω, -π, π}, Filling -> Axis, PlotRange -> All]

Estimate the power spectral density using specified window function values:

Wolfram Language code:

n = 10;
data = Range[n];
win = TukeyWindow[(#/2 * (n - 1)), (1/3)]& /@ Range[0, n - 1];

Compare to power spectral density with explicit TukeyWindow:

Wolfram Language code: PowerSpectralDensity[data, ω, win] === PowerSpectralDensity[data, ω, TukeyWindow]

Compare the smoothed spectrum to the original:

Wolfram Language code:

Plot[{Evaluate@PowerSpectralDensity[data, ω], Evaluate@PowerSpectralDensity[data, ω, TukeyWindow]}, {ω, -π, π}, Filling -> Axis, PlotRange -> All]

Compute the power spectral density, given a cutoff and a window function:

Wolfram Language code: data = Range[10];

Wolfram Language code: sspec = PowerSpectralDensity[data, ω, {5, BartlettWindow}]

Compare the smoothed spectrum to the original:

Wolfram Language code: Plot[{Evaluate@PowerSpectralDensity[data, ω], sspec}, {ω, -π, π}, Filling -> Axis, PlotRange -> All]

Random Processes (5)

Power spectral density for an ARProcess:

Wolfram Language code: proc[a_] = ARProcess[{a}, 1];

Wolfram Language code: Plot[PowerSpectralDensity[#, w], {w, -π, π}, Filling -> Axis, PlotRange -> All]& /@ {proc[0.8], proc[-0.8]}

Wolfram Language code: PowerSpectralDensity[proc[a], w]

Vector ARProcess:

Wolfram Language code:

a = {{1 / 7, 0}, {1 / 3, 0}};
Σ = {{1, -1 / 3}, {-1 / 3, 1}};
proc = ARProcess[{a}, Σ];

Wolfram Language code:

psd = PowerSpectralDensity[proc, ω];
psd//MatrixForm

Cross spectral density:

Wolfram Language code:

ParametricPlot[{Re[#], Im[#]}&@psd[[1, 2]], {ω, -π, π}, ColorFunction -> Function[{x, y, ω}, (ColorData["Rainbow"][ω])]]

Power spectral density for an MAProcess:

Wolfram Language code: proc[b_] = MAProcess[{b}, 1];

Wolfram Language code: Plot[PowerSpectralDensity[#, w], {w, -π, π}, Filling -> Axis]& /@ {proc[0.9], proc[-0.9]}

Wolfram Language code: PowerSpectralDensity[proc[b], w]

Vector MAProcess:

Wolfram Language code:

a = {{1 / 10, 0}, {1 / 3, 0}};
Σ = {{1, 1 / 3}, {1 / 3, 1}};
proc = MAProcess[{a}, Σ];

Wolfram Language code:

psd = PowerSpectralDensity[proc, ω];
psd//MatrixForm

Cross spectral density:

Wolfram Language code:

ParametricPlot[{Re[#], Im[#]}&@psd[[1, 2]], {ω, -π, π}, ColorFunction -> Function[{x, y, ω}, (ColorData["Rainbow"][ω])], AspectRatio -> 1]

Power spectral density for an ARMAProcess:

Wolfram Language code: proc[a_, b_] = ARMAProcess[{a}, {b}, 1];

Wolfram Language code: Plot[PowerSpectralDensity[#, w], {w, -π, π}, Filling -> Axis, PlotRange -> All]& /@ {proc[0.3, 0.2], proc[-0.3, -0.2]}

Wolfram Language code: PowerSpectralDensity[proc[a, b], w]

Vector ARMAProcess:

Wolfram Language code:

a = {{1 / 9, 0}, {1 / 3, 1 / 2}};
b = {{1, 3 / 4}, {1 / 2, -1 / 4}};
Σ = {{1, 1 / 3}, {1 / 3, 1}};
proc = ARMAProcess[{a}, {b}, Σ];

Wolfram Language code:

psd = PowerSpectralDensity[proc, ω];
psd//FullSimplify//MatrixForm

Cross spectral density:

Wolfram Language code:

ParametricPlot[{Re[#], Im[#]}&@psd[[1, 2]], {ω, -π, π}, ColorFunction -> Function[{x, y, ω}, (ColorData["Rainbow"][ω])], AspectRatio -> 1]

Power spectral density for a fractionally integrated time series:

Wolfram Language code: Plot[PowerSpectralDensity[FARIMAProcess[{.3, -.2}, -.4, {1}, 1], w], {w, -π, π}, Filling -> Axis, PlotRange -> All]

Wolfram Language code: PowerSpectralDensity[FARIMAProcess[{.3, .2}, .4, {1}, 1], w]

Vector FARIMAProcess:

Wolfram Language code:

a = {{1 / 9, 0}, {1 / 3, 1 / 2}};
b = {{1, 1 / 4}, {1 / 2, 0}};
Σ = {{1, 0}, {0, 1}};
proc = FARIMAProcess[{a}, {1 / 3, -1 / 4}, {b}, Σ];

Wolfram Language code:

psd = PowerSpectralDensity[proc, ω];
psd//FullSimplify//MatrixForm

Cross spectral density:

Wolfram Language code:

ParametricPlot[{Re[#], Im[#]}&@psd[[1, 2]], {ω, -π, π}, ColorFunction -> Function[{x, y, ω}, (ColorData["Rainbow"][ω])], AspectRatio -> 1]

Power spectral density for a seasonal time series:

Wolfram Language code: Plot[PowerSpectralDensity[SARMAProcess[{.2}, {.4}, {12, {.3}, {1}}, 1], w], {w, -π, π}, Filling -> Axis]

Wolfram Language code: PowerSpectralDensity[SARMAProcess[{a}, {b}, {12, {g}, {h}}, σ^2], ω]

Vector SARMAProcess:

Wolfram Language code:

a = {{1 / 9, 0}, {1 / 3, 1 / 2}};
b = {{1 / 3, 1 / 4}, {1 / 2, 0}};
Σ = {{1, 0}, {0, 1}};
proc = SARMAProcess[{}, {b}, {4, {b}, {}}, Σ];

Wolfram Language code:

psd = PowerSpectralDensity[proc, ω];
psd//FullSimplify//MatrixForm

Cross spectral density:

Wolfram Language code:

ParametricPlot[{Re[#], Im[#]}&@psd[[1, 2]], {ω, -π, π}, ColorFunction -> Function[{x, y, ω}, (ColorData["Rainbow"][ω])], AspectRatio -> 1]

Options (2)

The default value of FourierParameters:

Wolfram Language code: PowerSpectralDensity[Range[10], ω, FourierParameters -> {1, 1}]

Wolfram Language code: PowerSpectralDensity[Range[10], ω]

Wolfram Language code: % - %%

Change FourierParameters:

Wolfram Language code: PowerSpectralDensity[Range[10], ω, FourierParameters -> {-1, 1}]

It is the default value scaled:

Wolfram Language code: PowerSpectralDensity[Range[10], ω]

Wolfram Language code: % / %%//Simplify

Applications (1)

Use power spectral density for estimating time series processes:

Wolfram Language code: data = RandomFunction[ARProcess[{.3, .2}, 1], {10 ^ 3}];

Wolfram Language code: EstimatedProcess[data, ARProcess[2], ProcessEstimator -> "SpectralEstimator"]

Use a smoothing window:

Wolfram Language code: EstimatedProcess[data, ARProcess[2], ProcessEstimator -> {"SpectralEstimator", "Window" -> 5}]

Properties & Relations (11)

Power spectral density of a time series is a transform of the CovarianceFunction:

Wolfram Language code:

proc = ARProcess[c, {a}, σ^2];
cf = CovarianceFunction[proc, h]

Use FourierSequenceTransform:

Wolfram Language code: ft = FourierSequenceTransform[cf, h, z]

Compare to the power spectrum:

Wolfram Language code: PowerSpectralDensity[proc, z]

Wolfram Language code: FullSimplify[% - ft]

For a vector time series:

Wolfram Language code: proc = ARProcess[{{{1 / 3, 11 / 10}, {1 / 5, 1 / 10}}}, {{1, 0}, {0, 1}}];

Wolfram Language code: cov = CovarianceFunction[proc, h]//PiecewiseExpand;

Wolfram Language code: fst = FourierSequenceTransform[cov, h, w]

Wolfram Language code: fst - PowerSpectralDensity[proc, w]//Simplify

Power spectral density of data is a transform of the sample CovarianceFunction:

Wolfram Language code:

sample = Range[20];
n = Length[sample];
cov = CovarianceFunction[sample, {-n + 1, n - 1}];

Apply ListFourierSequenceTransform:

Wolfram Language code: ListFourierSequenceTransform[cov, w, -n + 1]//Simplify

Compare to SamplePowerSpectralDensity:

Wolfram Language code: PowerSpectralDensity[sample, w]//Simplify

Wolfram Language code: % - %%//FullSimplify

For a vector values time series:

Wolfram Language code:

data = RandomFunction[ARProcess[{{{.3, .1}, {.2, .1}}}, {{1, 0}, {0, 1}}], {10}];
n = data["PathLength"];

Wolfram Language code: cov = CovarianceFunction[data, {-n + 1, n - 1}]["Values"];

Wolfram Language code: lfst = Table[ListFourierSequenceTransform[cov[[All, i, j]], w, -n + 1], {i, 1, 2}, {j, 1, 2}]

Wolfram Language code: PowerSpectralDensity[data, w] - lfst//Simplify//Chop

Power spectrum of white noise:

Wolfram Language code: FourierSequenceTransform[Piecewise[{{1, x == 0}}, 0], x, w]

Compare to special case of an MAProcess:

Wolfram Language code: AbsoluteCorrelationFunction[MAProcess[{}, 1], h]

Wolfram Language code: PowerSpectralDensity[MAProcess[{}, 1], w]

Integrate to find the variance:

Wolfram Language code:

Integrate[PowerSpectralDensity[ARProcess[{a}, σ^2], w, FourierParameters -> {-1, 1}], {w, -π, π}, Assumptions -> σ > 0 && -1 ≤ a ≤ 1]

Compare to the variance of the time series:

Wolfram Language code: CovarianceFunction[ARProcess[{a}, σ^2], 0]

Wolfram Language code: % - %%//Simplify

Integrate to find the sample second moment:

Wolfram Language code:

data = TemporalData[TimeSeries, {{{0.30013118796154015, 0.4256944865636253, -1.2546039329698115, 
    0.09639437311894994, -1.9173223832884778, 0.16259596973060586, -0.001109688881539693, 
    0.6619148052033135, 2.232751727745077, 0.06794816837802453, 1. ... 5330770230747, 0.22924102933021667, 0.1509993303482705, 
    -0.4613641689767734}}, {{1, 100, 1}}, 1, {"Continuous", 1}, {"Discrete", 1}, 1, 
  {ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1}, ValueDimensions -> 1}}, False, 
 10.1];

Wolfram Language code: NIntegrate[PowerSpectralDensity[data, w, FourierParameters -> {-1, 1}], {w, -π, π}]

Compare to the sample second moment:

Wolfram Language code: CovarianceFunction[data, 0]

Wolfram Language code: % - %%//Chop

Power spectral density for harmonic frequencies is related to PeriodogramArray:

Wolfram Language code:

data = Range[20];
n = Length[data];
spec = Table[PowerSpectralDensity[data, 2π j / n], {j, 0, n - 1}]//N

Compare with PeriodogramArray:

Wolfram Language code: periodogram = PeriodogramArray[data]

For zero frequency:

Wolfram Language code: First[periodogram] - n * Mean[data] ^ 2

For nonzero frequencies:

Wolfram Language code: Rest[spec] - Rest[periodogram]//Chop

Diagonal elements of the power spectral density for vector data:

Wolfram Language code:

data = TemporalData[TimeSeries, {{{{0.22321490696028234, -0.09955096694103477}, 
    {-1.2036012283524826, 0.21356213785864084}, {-2.170758245872314, -2.172602353385546}, 
    {-2.0417769517615256, -0.9854153357651101}, {-1.8679940757580273, 0.11399938753 ... 3335557631977, 1.4574116828470514}, {-0.6330195762966424, -0.6361997810358189}}}, 
  {{0, 10, 1}}, 1, {"Continuous", 1}, {"Discrete", 1}, 2, 
  {ResamplingMethod -> {"Interpolation", InterpolationOrder -> 1}, ValueDimensions -> 2}}, False, 
 10.1];

Wolfram Language code: data["ValueDimensions"]

Wolfram Language code: psd = PowerSpectralDensity[data, ω];

Compare to univariate power spectral density for each data component:

Wolfram Language code: dens = Table[PowerSpectralDensity[data["Values"][[All, i]], ω], {i, 2}];

Wolfram Language code:

Table[Plot[{psd[[i, i]], dens[[i]]}, {ω, -π, π}, PlotStyle -> {Automatic, Dotted}, PlotLegends -> {"Vector Component", "Univariate"}], {i, 2}]

Power spectral density of a vector process is conjugate symmetric about zero:

Wolfram Language code: proc = ARProcess[{{{1 / 3, 1 / 10}, {1 / 5, 1 / 10}}}, {{1, 0}, {0, 1}}];

Wolfram Language code: Conjugate /@ PowerSpectralDensity[proc, -2]

Wolfram Language code: PowerSpectralDensity[proc, 2]

Wolfram Language code: % - %%

Power spectral density of a univariate process is symmetric about zero:

Wolfram Language code: proc = ARMAProcess[{a}, {b}, σ^2];

Wolfram Language code: PowerSpectralDensity[proc, w] == PowerSpectralDensity[proc, -w]

Power spectral density of a vector process is Hermitian:

Wolfram Language code: proc = ARMAProcess[{{{1 / 3, 1 / 10}, {1 / 5, 1 / 10}}}, {{{1 / 8, -3 / 10}, {3 / 5, 1 / 6}}}, {{1, -.3}, {-.3, 1}}];

Wolfram Language code: HermitianMatrixQ[PowerSpectralDensity[proc, RandomReal[{-π, π}]]]

Also non-negative definite:

Wolfram Language code: PositiveSemidefiniteMatrixQ[PowerSpectralDensity[proc, RandomReal[{-π, π}]]]

The magnitude of the sample cross spectral density is given by each component:

Wolfram Language code:

data = TemporalData[EventSeries, 
 {{{{0.4294925809742698255912415200470240013516218314553046886407`30., 
     -1.7225955472058514396692544902832573500729728703324221048`30.}, 
    {0.0548388212316276096199767991813010875210668838810634829247`30., 
     0 ... 255531611908244000251642161167860416322414`30., 
     0.0404300321804073415502129760700924797739015008346027682413`30.}}}, {{0, 20, 1}}, 1, 
  {"Continuous", 1}, {"Discrete", 1}, 2, {ResamplingMethod -> None, ValueDimensions -> 2}}, False, 
 10.1];

Wolfram Language code: data["ValueDimensions"]

Wolfram Language code: psd = PowerSpectralDensity[data, 3]//Chop

Wolfram Language code: Abs[psd[[1, 2]]] ^ 2 - Times@@Diagonal[psd]

The determinant of the sample power spectral density is constant equal to zero:

Wolfram Language code: FindMaximum[Abs@Det[PowerSpectralDensity[data, w]], {w, -3, 3}]

Use TransferFunctionModel to calculate PowerSpectralDensity of a time series:

Wolfram Language code:

proc = ARMAProcess[{a}, {b}, σ^2];
psd[ω_] = PowerSpectralDensity[proc, ω]

Define transfer function:

Wolfram Language code: g[z_] = TransferFunctionModel[proc][z][[1, 1]]

Calculate spectral density:

Wolfram Language code: σ^2g[Exp[I ω]]g[Exp[-I ω]]//ComplexExpand//TrigExpand//Simplify

Check:

Wolfram Language code: % - psd[ω]//Simplify

Neat Examples (1)

Plot a product of two power spectral densities in 3D:

Wolfram Language code:

proc = SARMAProcess[{.2}, {.4}, {24, {.3}, {.1}}, 1];
Plot3D[Evaluate[PowerSpectralDensity[proc, w] * PowerSpectralDensity[proc, z]], {w, -π, π}, {z, -π, π}, PlotRange -> All, ColorFunction -> "IslandColors", ViewPoint -> {-3, 3, .5}]

Top

More Learning

Tech Support

Wolfram Solutions

Wolfram Solutions For Education

Get Started

Grow Your Skills

Work with Us

Educational Programs for Adults

Educational Programs for Youth

Read

PowerSpectralDensity

Details and Options

Examples

Basic Examples (3)

Scope (14)

Empirical Estimates (4)

Smoothing (5)

Random Processes (5)

Options (2)

Applications (1)

Properties & Relations (11)

Neat Examples (1)

Text

CMS

APA

BibTeX

BibLaTeX

	{1,1}		default setting
	{-1,1}		often used for time series
	{a,b}		general setting

PowerSpectralDensity

Details and Options

Examples

Basic Examples (3)

Scope (14)

Empirical Estimates (4)

Smoothing (5)

Random Processes (5)

Options (2)

Applications (1)

Properties & Relations (11)

Neat Examples (1)

See Also

Related Guides

Related Links

History

Text

CMS

APA

BibTeX

BibLaTeX