pyyeti.psd.area

pyyeti.psd.area(spec)

Compute the area under a PSD random specification.

Parameters:

spec (2d ndarray or 2-element tuple/list) – If ndarray, its columns are [Freq, PSD1, PSD2, ... PSDn]. Otherwise, it must be a 2-element tuple or list, eg: (Freq, PSD) where PSD is: [PSD1, PSD2, ... PSDn]. In the second usage, PSD can be 1d; in the first usage, PSD is always considered 2d. The frequency vector must be monotonically increasing.

Returns:

area (1d array) – Area(s) under the PSD specification.

Notes

This routine assumes a constant db/octave slope for all segments. With this assumption, it computes the area for each segment with the following formula (derivation below). The segment goes from (f1, p1) to (f2, p2):

s = log(p2/p1)/log(f2/f1)
if s == -1:
    area_segment = p1*f1*log(f2/f1)
else:
    area_segment = (f2*p2 - f1*p1)/(s + 1)

Warning

This routine is only for specifications with all constant db/octave segments. Do not use for general freq vs. psd curves, such as output from analysis; use something like numpy.trapezoid() (or trapz() for NumPy versions < 2.0).

The following derives the equations for computing the area under the curve. Each segment is assumed to have a constant db/octave slope. We’ll start by computing the slope (\(m\)) of each segment. To get the slope, we’ll need the number of dbs \(d\) and the number of octaves \(n\). First, the number of dbs:

\[d = 10 \log_{10} (p_2 / p_1)\]

The number of octaves (\(n\)) is defined by \(f_2 = 2^n f_1\). Solving for \(n\) gives:

\[n = \frac{\log_{10} (f_2 / f_1)}{\log_{10} (2)}\]

Therefore, the slope \(m = d / n\) for the segment is:

\[m = 10 \log_{10} (2) \frac{\log_{10} (p_2 / p_1)} {\log_{10} (f_2 / f_1)}\]

To simplify the equations, we’ll define the variable \(s\) as:

\[s = \frac{m}{10 \log_{10} (2)} = \frac{\log_{10} (p_2 / p_1)} {\log_{10} (f_2 / f_1)}\]

Solving for \(p_2\) gives:

\[p_2 = p_1 ( f_2 / f_1 )^s\]

Since the segment has a constant db/octave slope, that equation can be generalized for any frequency value \(x\) from \(f_1\) to \(f_2\):

\[p(x) = p_1 ( x / f_1 )^s\]

The area under the segment is simply the integral of that equation:

\[area_{segment} = \int_{f_1}^{f_2} p_1 ( x / f_1 )^s dx\]

If \(s = -1\):

\[area_{segment} = p_1 f_1 \ln ( f_2 / f_1 )\]

Otherwise, if \(s \neq -1\):

\[\begin{split}\begin{aligned} area_{segment} &= \frac{p_1 ( f_2^{s+1} - f_1^{s+1} )} {f_1^s (s+1)} \\ &= ( p_2 f_2 - p_1 f_1 ) / (s + 1) \end{aligned}\end{split}\]

Examples

Compute a 3-sigma peak value given a test spec:

>>> import numpy as np
>>> from pyyeti import psd
>>> spec = np.array([[20, .0053],
...                  [150, .04],
...                  [600, .04],
...                  [2000, .0036]])
>>> 3*np.sqrt(psd.area(spec))
array([ 18.43...])

For comparison, use a brute-force technique:

>>> f = np.arange(20, 2000.1, 0.1)
>>> p = psd.interp(spec, f, linear=False)
>>> 3*np.sqrt(np.trapezoid(p, f, axis=0))
array([ 18.43...])