PSD – Power Spectral Density block

Block SymbolLicensing group: ADVANCED
PIC

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

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 [17], [18]). 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 N for FFT computing. The resulting PSD with npsd rows and m columns is referenced by upsd. Number of PSD frequencies npsd = N2 + 1 for even N or npsd = (N + 1)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 = 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 = 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(n) = 0.62 0.48 n N 1 0.5 + 0.38cos 2π n N 1 0.5, 0 n N 1.

iwin=2
Bartlett Window
w(n) = 2n N , 0 n N 2 , 2 2n N , N 2 n N.

iwin=3
Blackman Window
w(n) = 0.42 0.5cos 2πn L 1 + 0.08cos 4πn L 1, 0 n M 1,

where M is N2 when N is even and (N + 1)2 when N is odd.

iwin=4
Blackman-Harris Window
w(n) = a0 a1 cos 2πn L 1 + a2 cos 4πn L 1 a3 cos 6πn L 1, 0 n N 1

where a0 = 0.35875, a1 = 0.48829, a2 = 0.14128 and a3 = 0.01168.

iwin=5
Bohman Window
w(x) = 1 xcos π x + 1 πsin π x,  1 x 1

iwin=6
Chebyshev Window

Using Chebyshev polynomials of n-th order

Tn(x) cos ncos1(x), |x| 1, cosh ncosh1(x), |x| > 1,

the Fourier transform of the Chebyshev window is

W(k) = TN1 x0 cos(πkN) TN1(x0)

where

x0 = cosh cosh1 10r20 N 1 .

iwin=7
Flat Top Window
w(n) = a0 a1 cos 2πn L 1 + a2 cos 4πn L 1 a3 cos 6πn L 1 + a4 cos 8πn L 1,

where 0 n N 1 and a0 = 0.21557895, a1 = 0.41663158, a2 = 0.277263158, a3 = 0.083578947 and a4 = 0.006947368.

iwin=8
Gauss Window
w(n) = e1 2 α n (N1)22 = en2(2σ2)

where (N 1)2 n (N 1)2 and σ = (N 1)(2α) is the standard deviation of a Gaussian probability density function.

iwin=9
Hamming Window
w(n) = 0.54 0.46cos 2π n L 1, 0 n N 1,

iwin=10
Hann Window
w(n) = 0.5 1 cos 2π n L 1, 0 n N 1,

iwin=12
Nuttall’s Blackman-Harris Window
w(n) = a0 a1 cos 2πn L 1 + a2 cos 4πn L 1 a3 cos 6πn L 1, 0 n N 1

where a0 = 0.3635819, a1 = 0.4891775, a2 = 0.1365995 and a3 = 0.0106411.

iwin=13
Parzen Window
w(n) = 1 6 |n| N22 + 6 |n| N23, 0 |n| (N 1)4, 2 1 |n| N23, (N 1)4 < |n| (N 1)2.

iwin=14
Rectangular Window
w(n) = 1, 0 n 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(n) = 2n + 2 N + 1 , 0 n (N 1)2, 2 2n + 2 N + 1 , (N 1)2 < n (N 1).

For N even:

w(n) = 2n + 1 N , 0 n N2 1, 2 2n + 1 N , N2 n (N 1).

iwin=17
Tukey Window
w(x) = 1 2 1 + cos 2π r x r2, 0 x < r 2, 1, r 2 x < 1 r 2, 1 2 1 + cos 2π r x 1 + r2, 1 r 2 x < 1,

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

iwin=18
Exponential Window
w(x) = e|x| τ ,  1 x 1,

where τ is a window parameter (set in the parameter rwnd).

iwin=19
Welch Window
w(x) = 1 x2,  1 x 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.

Inputs

uc

Input reference to input data

Reference

uf

Input reference to internal FFT vector

Reference

upsd

Input reference to output PSD vector/matrix

Reference

urem

Input reference to remaining (not processed) data

Reference

uwnd

Input reference to window function vector

Reference

uwrk

Input reference to working vector/matrix

Reference

n

Number of signal samples for FFT segment

Long (I32)

fs

Sampling frequency in [Hz]

Double (F64)

START

Starting signal (rising edge)

Bool

HLD

The HLD=on freezes the computation

Bool

Parameters

mode

PSD computation mode   1  4 1

Long (I32)

1 ....

Real forward

2 ....

Real backward

3 ....

Complex forward

4 ....

Complex backward

mtrig

Computation trigger mode   1  3 1

Long (I32)

1 ....

Each period

2 ....

Enough samples

3 ....

START rising edge

4 ....

Complex overlap

iwnd

Window function   1  9 1

Long (I32)

1 ....

Bartlett-Hann

2 ....

Bartlett

3 ....

Blackman

4 ....

Blackman-Harris

5 ....

Bohman

6 ....

Chebyshev

7 ....

Flat top

8 ....

Gaussian

9 ....

Hamming

10 ...

Hann

11 ...

12 ...

Nuttall

13 ...

Parzen

14 ...

Rectangular

15 ...

16 ...

Triangular

17 ...

Tukey

18 ...

Exponential

19 ...

Welch

20 ...

External

lwnd

Parameter lwnd of some window   1  9 1

Long (I32)

rwnd

Parameter rwnd of some window   1  9 1

Long (I32)

navg

Number of FFT segments for averaging   1 10

Long (I32)

DB

Computed PSD is converted do decibels [dB]

Bool

Output

yc

Output reference to input data

Reference

yf

Output reference to internal FFT vector

Reference

ypsd

Output reference to output PSD vector/matrix

Reference

yrem

Output reference to remaining (not processed) data

Reference

ywnd

Output reference to window function vector

Reference

ywrk

Output reference to working vector/matrix

Reference

E

Error indicator

Bool

DONE

Averaging Completion Flag

Bool

iavg

Current index of segment for averaging

Long (I32)

lcomp

Number of successfully computed PSDs

Large (I64)

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