pkurtosis

(To be removed) Spectral kurtosis from signal or spectrogram

pkurtosis will be removed in a future release. Use spectralKurtosis instead. For more information, see Version History.

Syntax

sk = pkurtosis(x)

sk = pkurtosis(x,sampx)

sk = pkurtosis(xt)

sk = pkurtosis(___,window)

sk = pkurtosis(s,sampx,f,window)

[sk,fout] = pkurtosis(___)

[___,thresh] = pkurtosis(___,'ConfidenceLevel',p)

pkurtosis(___)

Description

sk = pkurtosis(x) returns the spectral kurtosis of vector x as the vector sk. pkurtosis uses normalized frequency (evenly spaced frequency vector spanning [0 π]) to compute the time values. pkurtosis computes the spectrogram of x using pspectrum with default window size (time resolution in samples), and 80% window overlap.

sk = pkurtosis(x,sampx) returns the spectral kurtosis of vector x sampled at rate or time interval sampx.

sk = pkurtosis(xt) returns the spectral kurtosis of single-variable timetable xt in the vector sk. xt must contain increasing finite time samples.

sk = pkurtosis(___,window) returns the spectral kurtosis using the time resolution specified in window for the pspectrum spectrogram computation. You can use window with any of the input arguments in previous syntaxes.

sk = pkurtosis(s,sampx,f,window) returns the spectral kurtosis using the spectrogram or power spectrogram s, along with:

• Sample rate or time, sampx, of the original time-series signal that was transformed to produce s

• Spectrogram frequency vector f

• Spectrogram time resolution window

Use this syntax when you want to customize the options for pspectrum, rather than accept the default pspectrum options that pkurtosis applies. You can specify sampx as empty to default to normalized frequency. Although window is optional for previous syntaxes, you must supply a value for window when using this syntax.

[sk,fout] = pkurtosis(___) returns the spectral kurtosis sk along with the frequency vector fout. You can use these output arguments with any of the input arguments in previous syntaxes.

[___,thresh] = pkurtosis(___,'ConfidenceLevel',p) returns the spectral kurtosis threshold thresh using the confidence level p. thresh represents the range within which the spectral kurtosis indicates a Gaussian stationary signal, at the optional confidence level p that you either specify or accept as default. Specifying p allows you to tune the sensitivity of the spectral kurtosis thresh results to behavior that is non-Gaussian or nonstationary. You can use the thresh output argument with any of the input arguments in previous syntaxes. You can also set the confidence level in previous syntaxes, but it has no effect unless you are returning or plotting thresh.

pkurtosis(___) plots the spectral kurtosis, along with confidence level and thresholds, without returning any data. You can use this syntax with any of the input arguments in previous syntaxes.

Examples

Plot Spectral Kurtosis of Nonstationary Signal Using Different Confidence Levels

Plot the spectral kurtosis of a chirp signal in white noise, and see how the nonstationary non-Gaussian regime can be detected. Explore the effects of changing the confidence level, and of invoking normalized frequency.

Create a chirp signal, add white Gaussian noise, and plot.

fs = 1000;
t = 0:1/fs:10;
f1 = 300;
f2 = 400;

xc = chirp(t,f1,10,f2);
x = xc + randn(1,length(t));

plot(t,x)
title('Chirp Signal with White Gaussian Noise')

Plot the spectral kurtosis of the signal.

pkurtosis(x,fs)
title('Spectral Kurtosis of Chirp Signal with White Gaussian Noise')

The plot shows a clear extended excursion from 300–400 Hz. This excursion corresponds to the signal component which represents the nonstationary chirp. The area between the two horizontal red-dashed lines represents the zone of probable stationary and Gaussian behavior, as defined by the 0.95 confidence interval. Any kurtosis points falling within this zone are likely to be stationary and Gaussian. Outside of the zone, kurtosis points are flagged as nonstationary or non-Gaussian. Below 300 Hz, there are a few additional excursions slightly above the above the zone threshold. These excursions represent false positives, where the signal is stationary and Gaussian, but because of the noise, has exceeded the threshold.

Investigate the impact of the confidence level by changing it from the default 0.95 to 0.85.

pkurtosis(x,fs,'ConfidenceLevel',0.85)
title('Spectral Kurtosis of Chirp Signal with Noise at Confidence Level of 0.85')

The lower confidence level implies more sensitive detection of nonstationary or non-Gaussian frequency components. Reducing the confidence level shrinks the thresh-delimited zone. Now the low-level excursions — false alarms — have increased in both number and amount. Setting the confidence level is a balancing act between achieving effective detection and limiting the number of false positives.

You can accurately determine and compare the zone width for the two cases by using the pkurtosis form that returns it.

[sk1,~,thresh95] = pkurtosis(x);
[sk2,~,thresh85] = pkurtosis(x,'ConfidenceLevel',0.85);
thresh = [thresh95 thresh85]

thresh = 1×2

    0.3578    0.2628

Plot the spectral kurtosis again, but this time, omit the sample time information so that pkurtosis plots normalized frequency.

pkurtosis(x,'ConfidenceLevel',0.85)
title('Spectral Kurtosis using Normalized Frequency')

The frequency axis has changed from Hz to a scale from 0 to π rad/sample.

Plot Spectral Kurtosis Using a Customized Window Size

The pkurtosis function uses the default pspectrum window size (time resolution). You can specify the window size to use instead. In this example, use the function kurtogram to return an optimal window size and use that result for pkurtosis.

Create a chirp signal with white Gaussian noise.

fs = 1000;
t = 0:1/fs:10;
f1 = 300;
f2 = 400;
x = chirp(t,f1,10,f2)+randn(1,length(t));

Plot the spectral kurtosis with the default window size.

pkurtosis(x,fs)
title('Spectral Kurtosis with Default Window Size')

Now compute the optimal window size using kurtogram.

kurtogram(x,fs)

The kurtogram plot also illustrates the chirp between 300 and 400 Hz, and shows that the optimum window size is 256. Feed w0 into pkurtosis.

w0 = 256;
pkurtosis(x,fs,w0)
title('Spectral Kurtosis with Optimum Window Size of 256')

The main excursion has higher kurtosis values. The higher values improve the differentiation between stationary and nonstationary components, and enhance your ability to extract the nonstationary component as a feature.

Plot Spectral Kurtosis Using a Customized Spectrogram

When using signal input data, pkurtosis generates a spectrogram by using pspectrum with default options. You can also create the spectrogram yourself if you want to customize the options.

Create a chirp signal with white Gaussian noise.

fs = 1000;
t = 0:1/fs:10;
f1 = 300;
f2 = 400;
x = chirp(t,f1,10,f2)+randn(1,length(t));

Generate a spectrogram that uses your specification for window, overlap, and number of FFT points. Then use that spectrogram in pkurtosis.

window = 256;
overlap = round(window*0.8);
nfft = 2*window;
[s,f,t] = spectrogram(x,window,overlap,nfft,fs);
figure
pkurtosis(s,fs,f,window)

The magnitude of the excursion is higher, and therefore better differentiated, than with default inputs in previous examples. However, the excursion magnitude here is not as high as it is in the kurtogram-optimized window example.

Input Arguments

x — Time-series signal

vector

Time-series signal from which pkurtosis returns the spectral kurtosis, specified as a vector.

sampx — Sample rate or sample time of signal

normalized frequency (default) | positive numeric scalar | duration scalar | numeric vector in seconds | duration array | datetime array

Sample rate or sample time, specified as one of the following::

• Positive numeric scalar — frequency in hertz

• duration scalar — time interval between consecutive samples of X

• Vector, duration array, or datetime array — time instant or duration corresponding to each element of x

When sampx represents a time vector, time samples can be nonuniform, with the pspectrum constraint that the median time interval and the mean time interval must obey:

$$\frac{1}{100}<\frac{\text{Median time interval}}{\text{Mean time interval}}<100.$$

If you specify sampx as empty, then pkurtosis uses normalized frequency. In other words, it assumes an evenly spaced frequency vector spanning [0 π].

xt — Signal timetable

timetable

Signal timetable from which pkurtosis returns the spectral kurtosis, specified as a timetable that contains a single variable with a single column. xt must contain increasing, finite row times. If the timetable has missing or duplicate time points, you can fix it using the tips in Clean Timetable with Missing, Duplicate, or Nonuniform Times. xt can be nonuniformly sampled, with the pspectrum constraint that the median time interval and the mean time interval must obey:

$$\frac{1}{100}<\frac{\text{Median time interval}}{\text{Mean time interval}}<100.$$

window — Window time resolution

positive scalar

Window time resolution to use for the internal pspectrum spectrogram computation, specified as a positive scalar in samples. window is required for syntaxes that use an existing spectrogram as input, and optional for the rest. You can use the function kurtogram to determine the optimal window size to use. pspectrum uses 80% overlap by default.

s — Spectrogram or power spectrogram of signal

complex matrix | real nonnegative matrix

Power spectrogram or spectrum of a signal, specified as a matrix (spectrogram) or a column vector (spectrum).

• If s is complex, then pkurtosis treats s as a short-time Fourier transform (STFT) of the original signal (spectrogram).

• If s is real, then pkurtosis treats s as the square of the absolute values of the STFT of the original signal (power spectrogram). Thus, every element of s must be nonnegative.

If you specify s, pkurtosis uses s rather than generate its own spectrogram or power spectrogram. For an example, see Plot Spectral Kurtosis Using a Customized Spectrogram.

Data Types: single | double
Complex Number Support: Yes

f — Frequencies for s

vector

Frequencies for spectrogram or power spectrogram s when s is supplied explicitly to pkurtosis, specified as a vector in hertz. The length of f must be equal to the number of rows in s.

'ConfidenceLevel', p — Confidence level

0.95 (default) | [0 to 1]

Confidence level used to determine whether signal is likely to be Gaussian and stationary, specified as a numeric scalar value from 0 to 1. p influences the thresh range where the spectral kurtosis value indicates a Gaussian and stationary signal. The confidence level therefore provides a detection-sensitivity tuning parameter. Kurtosis values outside of this range indicate, with a probability of (1-p), non-Gaussian or nonstationary behavior. For an example, see Plot Spectral Kurtosis of Nonstationary Signal Using Different Confidence Levels.

Output Arguments

sk — Spectral kurtosis

double vector

Spectral Kurtosis, returned as a double vector. The spectral kurtosis is a statistical quantity that contains low values where data is stationary and Gaussian, and high values where transients occur. One use of the spectral kurtosis is to detect and locate nonstationary or non-Gaussian behavior that could result from faults or degradation. The high-valued kurtosis data reveals such signal components.

fout — frequencies for sk

double vector

Frequencies associated with sk values, returned as a vector in hertz.

thresh — Spectral kurtosis band size for stationary Gaussian behavior

scalar

Spectral kurtosis band size for stationary Gaussian behavior, returned as a numeric scalar representing the thickness of the band centered at the sk = 0 line, given confidence level p. Excursions outside the thresh-delimited band indicate possible nonstationary or non-Gaussian behavior. Confidence level p directly influences the thickness of the band and the sensitivity of the results. For an example, see Plot Spectral Kurtosis of Nonstationary Signal Using Different Confidence Levels.

More About

Spectral Kurtosis

Spectral kurtosis (SK) is a statistical tool that can indicate and pinpoint nonstationary or non-Gaussian behavior in the frequency domain, by taking:

• Small values at frequencies where only stationary Gaussian noise is present

• High positive values at frequencies where transients occur

This capability makes SK a powerful tool for detecting and extracting signals associated with faults in rotating mechanical systems. On its own, SK can identify features or conditional indicators for fault detection and classification. As preprocessing for other tools such as envelope analysis, SK can supply key inputs such as optimal band [1], [2],.

The spectral kurtosis, or K(f), of a signal x(t) can be computed based on the short-time Fourier transform (STFT) of the signal, S(t,f):

$$S\left(t,f\right)={\displaystyle {\int }_{-\infty }^{\infty }x\left(\tau \right)w\left(\tau -t\right)e^{-j2\pi f\tau }d\tau },$$

where w(t) is the window function used in the STFT. K(f) is calculated as

$$K(f)=\frac{{\langle {\left|S(t,f)\right|}^4\rangle }_t}{{\langle {\left|S(t,f)\right|}^2\rangle }_t^2}-2,f\ne 0,$$

where 〈·〉_t is the time-average operator.

If the signal x(t) contains only stationary Gaussian noise, then K(f) at each frequency f has an asymptotic normal distribution with 0 mean and variance 4/M, where M is the number of elements along the time axis in S(t,f). Hence, a statistical threshold s_α given a confidence level α is

$$s_{\alpha }={\Phi }^{-1}\left(\alpha \right)\frac{2}{\sqrt{M}},$$

where $${\Phi }^{-1}(\alpha )=\sqrt{2}{\mathrm{erf}}^{-1}(2\alpha -1)$$ is the quantile function of the standard normal distribution.

It is important to note that the STFT window length N_w directly drives frequency resolution, which is f_s/N_w, where f_s is the sample rate. The window size must be shorter than the spacing between transient impulses but longer than the individual transient impulses.

References

[1] Antoni, J. "The Spectral Kurtosis: A Useful Tool for Characterising Non-Stationary Signals." Mechanical Systems and Signal Processing. Vol. 20, Issue 2, 2006, pp. 282–307.

[2] Antoni, J., and R. B. Randall. "The Spectral Kurtosis: Application to the Vibratory Surveillance and Diagnostics of Rotating Machines." Mechanical Systems and Signal Processing. Vol. 20, Issue 2, 2006, pp. 308–331.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The pkurtosis function supports GPU array input with these usage notes and limitations:

• If the input is a pure tone signal, the result may show small numerical differences with respect to in-memory computation of spectral kurtosis.

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2018a

R2024b: pkurtosis will be removed

The pkurtosis function will be removed in a future release. Use spectralKurtosis instead. You must update your code to use spectralKurtosis.

Consider a signal x sampled at a rate fs:

fs = 1000;
Ts = seconds(1/fs);
t = (0:1/fs:10)';

x = chirp(t,300,10,400);
xT = timetable(x,SampleRate=fs);

function y = winlen(x)
    wdiv = 2.^[1 3:7];
    y = ceil(x/wdiv(find(x < 2.^[6 11:14 Inf],1)));
end

wnd = 500;

[Sstft,Fstft] = stft(x,fs, ...
    FrequencyRange="onesided", ...
    Window=triang(600),FFTLength=1512,OverlapLength=421);

To get the spectral kurtosis with spectralKurtosis, make these updates to your code..

Original Code in R2024a or Earlier	Updated Code in R2024b
sk = pkurtosis(xT);	[S,F] = pspectrum(xT,"spectrogram", ... FrequencyResolution= ... xT.Properties.SampleRate/winlen(height(xT)), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x);	[S,F] = pspectrum(x,1,"spectrogram", ... FrequencyResolution= ... 1/winlen(length(x)), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,fs);	[S,F] = pspectrum(x,fs,"spectrogram", ... FrequencyResolution= ... fs/winlen(length(x)), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,Ts);	[S,F] = pspectrum(x,Ts,"spectrogram", ... FrequencyResolution= ... 1/(seconds(Ts)*winlen(length(x))), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,t);	[S,F] = pspectrum(x,t,"spectrogram", ... FrequencyResolution= ... 1/(mean(diff(t))*winlen(length(x))), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(xT,wnd);	[S,F] = pspectrum(xT,"spectrogram", ... FrequencyResolution= ... xT.Properties.SampleRate/wnd, ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,[],wnd);	[S,F] = pspectrum(x,1,"spectrogram", ... FrequencyResolution=1/wnd, ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,fs,wnd);	[S,F] = pspectrum(x,fs,"spectrogram", ... FrequencyResolution=fs/wnd, ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,Ts,wnd);	[S,F] = pspectrum(x,Ts,"spectrogram", ... FrequencyResolution= ... 1/(seconds(Ts)*wnd), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(x,t,wnd);	[S,F] = pspectrum(x,t,"spectrogram", ... FrequencyResolution= ... 1/(mean(diff(t))*wnd), ... OverlapPercent=80); sK = spectralKurtosis(S,F,Scaled=false);
sk = pkurtosis(Sstft,fs,Fstft,wnd);	sK = spectralKurtosis(abs(Sstft).^2,Fstft, ... Scaled=false);
sk = pkurtosis(Sstft,fs,Fstft,wnd);	sK = spectralKurtosis(x,fs, ... Window=triang(600),FFTLength=1512, ... OverlapLength=421,Scaled=false);

• You can use the ConfidenceLevel name-value argument in spectralKurtosis just like in pkurtosis.

R2023b: Use timetables for code generation

The pkurtosis function supports MATLAB^® timetables for code generation.