Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language

Function Repository Resource:

 WelchSpectralEstimate

Source Notebook

Fast numerical estimation of the power spectral density or the cross spectral density

Contributed by: Julien Kluge
Quantum Optical Metrology; Joint Lab Integrated Quantum Sensors
Department of Physics
Humboldt-Universität zu Berlin
julien@physik.hu-berlin.de

ResourceFunction["WelchSpectralEstimate"][data,rate]

estimates the data power spectral density at a given rate.

ResourceFunction["WelchSpectralEstimate"][data1,data2,rate]

estimates the data cross spectral density at a given rate.

Details and Options

ResourceFunction["WelchSpectralEstimate"] accepts numerical data, TemporalData objects or TimeSeries objects.
ResourceFunction["WelchSpectralEstimate"] always assumes evenly spaced data points. For TimeSeries or TemporalData sets, the rate is estimated as the number of data points divided by the duration between the first and last time stamp.
Supported options are:
"SegmentSize"Automaticthe number of datapoints to use in a Fourier transform
"OverlapOffset"Automaticthe offset between segments to average or a fraction of the possible offset (defaults to half of the "SegmentSize")
FourierParameters{0,1}Fourier convention parameters
"Window"Automaticthe spectral window to apply to the data (if appropriate, the HannWindow will be chosen)
"OneSided"TrueTrue to evaluate only the first half of the Fourier transform, otherwise False for evaluating the whole range
"Reduction"Meanthe function to average the segments together
"DensityScaling"TrueTrue for scaling as a density, False for spectrogram scaling
"Detrend""Constant"the level of detrending to apply to the data segments, before transformation

Examples

Basic Examples (1) 

Calculate the spectral density of a random data sample:

In[1]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogLogPlot[ResourceFunction["WelchSpectralEstimate"][data, 100.0], Joined -> True]
Out[2]=

Scope (2) 

Evaluate the Cross Spectral Density (CSD) of two sinusoidal, noisy datasets and display amplitude and phase:

In[3]:=
data1 = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; data2 = 0.3 Sin[15*2 Pi*Range[0, 10, 0.01] + 0.5] + Sin[10*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; csd = ResourceFunction["WelchSpectralEstimate"][data1, data2, 100.0]; {ListLogLogPlot[Abs[csd], Joined -> True], ListLogLinearPlot[{#1, Arg@#2} & @@@ csd, Joined -> True]}
Out[6]=

Compare the influence of different spectral windows on the PSD:

In[7]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.001], 1001]; windows = {DirichletWindow, HannWindow, HammingWindow, BlackmanHarrisWindow, KaiserWindow}; ListLogLogPlot[ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Window" -> #] & /@ windows, Joined -> True, PlotLegends -> windows]
Out[9]=

Options (12) 

SegmentSize (1) 

"SegmentSize" sets number of data points from which to calculate the Fourier transform; smaller numbers means more averaging but lower resolution:

In[10]:=
data = Sin[15*2 Pi*Range[0, 10, 0.001]] + RandomVariate[NormalDistribution[0, 0.1], 10001]; ListLogLogPlot[ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "SegmentSize" -> #], Joined -> True, PlotLabel -> #] & /@ {200, 500, 1000, 5000, 10000}
Out[11]=

OverlapOffset (2) 

"OverlapOffset" allows to set the distance of the starting point between two partitioned segments; a smaller integer number means closer overlaps and thus more averaging:

In[12]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogLogPlot[ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "OverlapOffset" -> #], Joined -> True, PlotLabel -> #] & /@ {1, 511}
Out[13]=

"OverlapOffset" also allows the specification of an overlap fraction of the number of data points minus the segment size, as a real number between zero and one:

In[14]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogLogPlot[ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "OverlapOffset" -> #], Joined -> True, PlotLabel -> #] & /@ {0.01, 0.99}
Out[15]=

FourierParameters (1) 

Set the FourierParameters option for the respective Fourier transforms:

In[16]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogLogPlot[ ResourceFunction["WelchSpectralEstimate"][data, 100.0, FourierParameters -> {1, -1}], Joined -> True]
Out[17]=

Window (3) 

Use the "Window" option to specify windowing functions:

In[18]:=
data = Sin[15*2 Pi*Range[0, 10, 0.001]] + RandomVariate[NormalDistribution[0, 0.1], 10001]; ListLogLogPlot[{ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Window" -> DirichletWindow], ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Window" -> HannWindow] }, Joined -> True, PlotLegends -> {"Dirichlet", "Hann"}]
Out[19]=

A list can be applied as the window function, provided the length is equal to its segment size:

In[20]:=
data = Sin[15*2 Pi*Range[0, 10, 0.001]] + RandomVariate[NormalDistribution[0, 0.1], 10001]; segmentSize = 1024; windowPoints = PDF[NormalDistribution[0, 0.1], Subdivide[-1, 1, segmentSize - 1]]; ListLogLogPlot[{ResourceFunction["WelchSpectralEstimate"][data, 100.0, "SegmentSize" -> segmentSize, "Window" -> windowPoints]}, Joined -> True]
Out[23]=

Set a custom function as the window:

In[24]:=
data = Sin[15*2 Pi*Range[0, 10, 0.001]] + RandomVariate[NormalDistribution[0, 0.1], 10001]; ListLogLogPlot[{ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Window" -> (Exp[-50 #^2] &)], ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Window" -> (HannWindow[#, 0.4] &)] }, Joined -> True]
Out[25]=

OneSided (1) 

Set "OneSided" to False to get a two sided evaluation:

In[26]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogPlot[{ResourceFunction["WelchSpectralEstimate"][data, 100.0, "OneSided" -> False]}, Joined -> True, PlotRange -> All]
Out[27]=

Reduction (1) 

Specify Median as a reduction function:

In[28]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogLogPlot[{ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Reduction" -> Median]}, Joined -> True]
Out[29]=

DensityScaling (1) 

Set "DensityScaling" to False to get a spectrum scaling:

In[30]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.1], 1001]; ListLogLogPlot[{ResourceFunction["WelchSpectralEstimate"][data, 100.0, "DensityScaling" -> False]}, Joined -> True, PlotRange -> All]
Out[31]=

Detrend (2) 

The "Detrend" option can be used to detrend the initial data segments before transformations occur. The option can be specified with string arguments for constant or linear detrends:

In[32]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/08505304-8fb8-4f9a-8650-67d7c32ff920"]
Out[33]=

A list of integers can be specified instead which polynomial orders should be detrended:

In[34]:=
data = Sin[15*2 Pi*Range[0, 10, 0.001]] + 10*Range[0, 10, 0.001] + 10*Range[0, 10, 0.001]^2 + 10*Range[0, 10, 0.001]^3 + 100 + RandomVariate[NormalDistribution[0, 0.1], 10001]; ListLogLogPlot[{ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Detrend" -> None], ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Detrend" -> {0, 1, 2, 3}] }, Joined -> True, PlotLegends -> {"None", "Detrend: \!\(\*SuperscriptBox[\(x\), \(3\)]\), \!\(\*SuperscriptBox[\(x\), \(2\)]\), x, \!\(\*SuperscriptBox[\(x\), \(0\)]\)"}]
Out[35]=

Applications (2) 

Investigate the frequency response on white noise data for different kind of filters:

In[36]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/07e1b8a6-d841-4a80-99fc-593947cfa61c"]
Out[40]=

Estimate the system transfer function of a set state and its measured response:

In[41]:=
time = Range[0, 10, 0.05]; data = Table[Round[Mod[x, 1]] - 0.5, {x, time}]; dataResponse = NDSolveValue[{u'[ t] == (Interpolation[Transpose[{time, data}], InterpolationOrder -> 0][t] - u[t]), u[0] == 0}, u, {t, Min[time], Max[time]}] /@ time; csd = ResourceFunction["WelchSpectralEstimate"][data, dataResponse, 0.05, "Window" -> KaiserWindow]; psd = ResourceFunction["WelchSpectralEstimate"][data, 0.05, "Window" -> KaiserWindow]; transferFunction = Transpose[{csd[[All, 1]], csd[[All, 2]]/psd[[All, 2]]}]; phase = Arg[transferFunction[[All, 2]]]; { ListLogLogPlot[{#1, Abs@#2} & @@@ transferFunction, Joined -> True, PlotLabel -> "Abs of frequency response", ImageSize -> Medium], ListLogLinearPlot[{#1, 180/Pi*Arg@#2} & @@@ transferFunction, Joined -> True, PlotRange -> All, PlotLabel -> "Phase of frequency response", ImageSize -> Medium] }
Out[48]=

Possible Issues (1) 

Even for complex input, the "OneSided" option will still remain True by default and thus return only the front of the Fourier transform:

In[49]:=
ResourceFunction["WelchSpectralEstimate"][ RandomVariate[NormalDistribution[0, 0.001], 10] + I*RandomVariate[NormalDistribution[0, 0.001], 10], 1.0]
Out[49]=

Neat Examples (1) 

Estimate the uncertainty of the Fourier deviations by just choosing a different reduction function:

In[50]:=
data = Sin[15*2 Pi*Range[0, 10, 0.01]] + RandomVariate[NormalDistribution[0, 0.001], 1001]; ListLogLogPlot[ ResourceFunction["WelchSpectralEstimate"][data, 100.0, "Reduction" -> MeanAround]]
Out[51]=

Publisher

Julien Kluge

Version History

  • 1.0.0 – 10 February 2021

Source Metadata

  • Citation:
    • Welch, P. D., "The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms." IEEE Transactions on Audio and Electroacoustics, AU-15 (2), 70–73, 1967. DOI:10.1109/TAU.1967.1161901

License Information