PSD – Power Spectral Density

Block SymbolLicensing group: ADVANCED

Function Description
The function block PSD computes a power spectral density (PSD, periodogram) using the PocketFFT package [17]. The PocketFFT C99 implementation is based on FFTPack (in Fortran) [18], which is based on the chapter in [19].

PSD is most often calculated from segments of length N samples of the input signal. The data of these segments is multiplied by a window function (e.g. Hamming window), the FFT is calculated from them, and the PSD of each segment is calculated from its result. Finally, the navg of such results is averaged (the navg parameter). The input signal can be divided into segments either sequentially or with a half-overlap (Welch method [20], [21]). The resulting PSD can be calculated in decibels (dB) if the parameter DB = on.

Input data (a vector or a matrix) of the block are referenced by the uc input. If the input uc refers to a column vector (number of columns m = $1$) with the number of elements n, the PSD will be calculated from a single signal. If uc refers to a matrix with n rows and m columns, the PSD will be computed m times (for each column). If the input n > 0 then the number of samples per segment N = n else N is set to the number of rows of the array referenced by uc.

The input uf references the internal vector of dimension $\ge \mathtt{\text{N}}$ for FFT computing. The resulting PSD with npsd rows and m columns is referenced by upsd. Number of PSD frequencies $\mathtt{\text{npsd}}=\mathtt{\text{N}}∕2+1$ for even N or $\mathtt{\text{npsd}}=\left(\mathtt{\text{N}}+1\right)∕2$ for odd N. Working array referenced by uwrk has the same dimension as the array referenced by upsd. If the input data contains more elements than can be processed in one execution of the PSD block, the remaining data is copied into an array referenced by urem that has N rows and m columns.

Input data is processed according to the mode and mtrig parameters (see below). The mode parameter specifies the way the PDS is calculated, whether the data in the input vector referenced by uc is real or complex, and whether segment overlapped by half will be used. The input data sampled with the frequency connected to the input fs may be acquired in a different task (e.g. with a significantly shorter sampling period) than the one in which the PSD is calculated. The method of synchronizing the start of the calculation is provided by the mtrig trigger mode. The calculation can be triggered every PSD block execution period or when a sufficient number of samples are available (at least equal to N) or by the rising edge of the START input.

A vector of dimension N is connected to the uwnd input to store the window function selected by the iwin parameter. Window functions have (optionally) an integer parameter lwnd and/or a real parameter rwnd. Some window functions are used in two variants:

Symmetric (for lwnd=0). The window of length $L=\mathtt{\text{N}}$ samples is computed, the first and last elements are equal. This variant is suitable for FIR filter design.

Periodic (for lwnd=1). The window of length $L=\mathtt{\text{N}}+1$ is computed but only the first N samples is stored. The periodic version is the preferred method for spectral analysis because the discrete Fourier transform assumes periodic extension of the input vector.

Window functions are defined by the following expressions:

iwin=1

Bartlett-Hahn Window

$$w\left(n\right)=0.62-0.48\left|\left(\frac{n}{N-1}-0.5\right)\right|+0.38cos\left[2\pi \left(\frac{n}{N-1}-0.5\right)\right],0\le n\le N-1.$$

iwin=2

Bartlett Window

$$w\left(n\right)=\left\{\begin{array}{cc}\frac{2n}{N},\hfill & 0\le n\le \frac{N}{2},\hfill \\ 2-\frac{2n}{N},\hfill & \frac{N}{2}\le n\le N.\hfill \end{array}\right.$$

iwin=3

Blackman Window

$$w\left(n\right)=0.42-0.5cos\left(\frac{2\pi n}{L-1}\right)+0.08cos\left(\frac{4\pi n}{L-1}\right),0\le n\le M-1,$$

where $M$ is $N∕2$ when $N$ is even and $\left(N+1\right)∕2$ when $N$ is odd.

iwin=4

Blackman-Harris Window

$$w\left(n\right)=a_0-a_{1}cos\left(\frac{2\pi n}{L-1}\right)+a_{2}cos\left(\frac{4\pi n}{L-1}\right)-a_{3}cos\left(\frac{6\pi n}{L-1}\right),0\le n\le N-1$$

where $a_0=0.35875$, $a_1=0.48829$, $a_2=0.14128$ and $a_3=0.01168$.

iwin=5

Bohman Window

$$w\left(x\right)=\left(1-\left|x\right|\right)cos\left(\pi \left|x\right|\right)+\frac{1}{\pi }sin\left(\pi \left|x\right|\right),-1\le x\le 1$$

iwin=6

Chebyshev Window

Using Chebyshev polynomials of $n$-th order

$$T_n\left(x\right)\triangleq \left\{\begin{array}{cc}cos\left[n{cos}^{-1}\left(x\right)\right],\hfill & \left|x\right|\le 1,\hfill \\ cosh\left[n{cosh}^{-1}\left(x\right)\right],\hfill & \left|x\right|>1,\hfill \end{array}\right.$$

the Fourier transform of the Chebyshev window is

$$W\left(k\right)=\frac{T_{N-1}\left[x_{0}cos\left(\pi k∕N\right)\right]}{T_{N-1}\left(x_0\right)}$$

where

$$x_0=cosh\left[\frac{{cosh}^{-1}\left(10^{r∕20}\right)}{N-1}\right].$$

iwin=7

Flat Top Window

$$w\left(n\right)=a_0-a_{1}cos\left(\frac{2\pi n}{L-1}\right)+a_{2}cos\left(\frac{4\pi n}{L-1}\right)-a_{3}cos\left(\frac{6\pi n}{L-1}\right)+a_{4}cos\left(\frac{8\pi n}{L-1}\right),$$

where $0\le n\le N-1$ and $a_0=0.21557895$, $a_1=0.41663158$, $a_2=0.277263158$, $a_3=0.083578947$ and $a_4=0.006947368$.

iwin=8

Gauss Window

$$w\left(n\right)=e^{-\frac{1}{2}{\left[\alpha \frac{n}{\left(N-1\right)∕2}\right]}^2}=e^{-n^2∕\left(2{\sigma }^2\right)}$$

where $-\left(N-1\right)∕2\le n\le \left(N-1\right)∕2$ and $\sigma =\left(N-1\right)∕\left(2\alpha \right)$ is the standard deviation of a Gaussian probability density function.

iwin=9

Hamming Window

$$w\left(n\right)=0.54-0.46cos\left(2\pi \frac{n}{L-1}\right),0\le n\le N-1,$$

iwin=10

Hann Window

$$w\left(n\right)=0.5\left[1-cos\left(2\pi \frac{n}{L-1}\right)\right],0\le n\le N-1,$$

iwin=12

Nuttall’s Blackman-Harris Window

$$w\left(n\right)=a_0-a_{1}cos\left(\frac{2\pi n}{L-1}\right)+a_{2}cos\left(\frac{4\pi n}{L-1}\right)-a_{3}cos\left(\frac{6\pi n}{L-1}\right),0\le n\le N-1$$

where $a_0=0.3635819$, $a_1=0.4891775$, $a_2=0.1365995$ and $a_3=0.0106411$.

iwin=13

Parzen Window

$$w\left(n\right)=\left\{\begin{array}{cc}1-6{\left(\frac{\left|n\right|}{N∕2}\right)}^2+6{\left(\frac{\left|n\right|}{N∕2}\right)}^3,\hfill & 0\le \left|n\right|\le \left(N-1\right)∕4,\hfill \\ 2{\left(1-\frac{\left|n\right|}{N∕2}\right)}^3,\hfill & \left(N-1\right)∕4<\left|n\right|\le \left(N-1\right)∕2.\hfill \end{array}\right.$$

iwin=14

Rectangular Window

$$w\left(n\right)=1,0\le n\le N-1.$$

Multiplication of all input data elements by 1 is not performed in this block, therefore the uwnd input may not be connected in this case.

iwin=16

Triangular Window

For $N$ odd:

$$w\left(n\right)=\left\{\begin{array}{cc}\frac{2n+2}{N+1},\hfill & 0\le n\le \left(N-1\right)∕2,\hfill \\ 2-\frac{2n+2}{N+1},\hfill & \left(N-1\right)∕2<n\le \left(N-1\right).\hfill \end{array}\right.$$

For $N$ even:

$$w\left(n\right)=\left\{\begin{array}{cc}\frac{2n+1}{N},\hfill & 0\le n\le N∕2-1,\hfill \\ 2-\frac{2n+1}{N},\hfill & N∕2\le n\le \left(N-1\right).\hfill \end{array}\right.$$

iwin=17

Tukey Window

$$w\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}\left\{1+cos\left(\frac{2\pi }{r}\left[x-r∕2\right]\right)\right\},\hfill & 0\le x<\frac{r}{2},\hfill \\ 1,\hfill & \frac{r}{2}\le x<1-\frac{r}{2},\hfill \\ \frac{1}{2}\left\{1+cos\left(\frac{2\pi }{r}\left[x-1+r∕2\right]\right)\right\},\hfill & 1-\frac{r}{2}\le x<1,\hfill \end{array}\right.$$

where $r$ is so called cosine fraction (set in the parameter rwnd, default value $0.5$).

iwin=18

Exponential Window

$$w\left(x\right)=e^{-\frac{\left|x\right|}{\tau }},-1\le x\le 1,$$

where $\tau$ is a window parameter (set in the parameter rwnd).

iwin=19

Welch Window

$$w\left(x\right)=1-x^2,-1\le x\le 1.$$

iwin=20

Externally Defined Window

If the user is not satisfied with any of the predefined windows, he/she can set his own window consisting of N samples in the vector referenced by uwnd before executing the PSD block.

The HLD input allows the user to temporarily stop the PSD calculation.

This block does not propagate the signal quality. More information can be found in the 1.4 section.

2024 © REX Controls s.r.o., www.rexygen.com