This gallery provides you with an overview of the timefrequency analysis features available in Signal Processing Toolbox™ and Wavelet Toolbox™. The descriptions and usage examples present various methods that you can use for your signal analysis.
Method  Features  Invertible  Examples 


 
 Yes  
 No  

 
 Yes  
 No 
The shorttime Fourier transform is a linear timefrequency representation useful in the analysis of nonstationary multicomponent signals.
The shorttime Fourier transform is invertible.
The spectrogram is the magnitude squared of the STFT.
You can compute the crossspectrogram of two signals to look for similarities in timefrequency space.
The persistence spectrum of a signal is a timefrequency view that shows the percentage of the time that a given frequency is present in a signal. The persistence spectrum is a histogram in powerfrequency space. The longer a particular frequency persists in a signal as the signal evolves, the higher its time percentage and thus the brighter or "hotter" its color in the display.
The applications of this timefrequency method include, but are not limited to:
Audio signal processing: Fundamental frequency estimation, cross synthesis, spectral envelope extraction, timescale modification, timestretching, and pitch shifting. (See Phase Vocoder with Different Synthesis and Analysis Windows for more details.)
Crack detection: Detect cracks in aluminum plates using dispersion curves of ultrasonic Lamb waves.
Sensor array processing: Sonar exploration, geophysical exploration, and beamforming.
Digital communications: Detection of frequency hopping signal.
stft
computes the shorttime Fourier transform. To invert the shorttime Fourier transform,
use the istft
function.
dlstft
computes the deep learning shorttime Fourier transform. You must have Deep Learning Toolbox™ installed.
pspectrum
or spectrogram
computes the
spectrogram.
xspectrogram
computes the
crossspectrogram of two signals.
You can also use the spectrogram view in Signal Analyzer to view the spectrogram of a signal.
Use the persistence spectrum option in pspectrum
or Signal Analyzer to identify signals hidden in
other signals.
Generate a signal sampled at 5 kHz for 4 seconds. The signal consists of a set of pulses of decreasing duration separated by regions of oscillating amplitude and fluctuating frequency with an increasing trend.
fs = 5000; t = 0:1/fs:41/fs; x = 10*besselj(0,1000*(sin(2*pi*(t+2).^3/60).^5));
Compute and plot the shorttime Fourier transform of the signal. Window the signal with a 200sample Kaiser window with shape factor $\beta =30$.
stft(x,fs,'Window',kaiser(200,30))
Load an audio signal that contains two decreasing chirps and a wideband splatter sound.
load splat
Set the overlap length to 96 samples. Plot the shorttime Fourier transform.
stft(y,Fs,'OverlapLength',96)
Load a file that contains audio data from a Pacific blue whale, sampled at 4 kHz. The file is from the library of animal vocalizations maintained by the Cornell University Bioacoustics Research Program. The time scale in the data is compressed by a factor of 10 to raise the pitch and make the calls more audible.
[w,fs] = audioread('bluewhale.wav');
Compute the spectrogram of the whale song with an overlap percentage equal to eighty percent. Set the minimum threshold for the spectrogram to 50
dB.
pspectrum(w,fs,'spectrogram','Leakage',0.2,'OverlapPercent',80,'MinThreshold',50)
Load an interference narrowband signal embedded within a broadband signal.
load TransientSig
Compute the persistence spectrum of the signal. Both signal components are clearly visible.
pspectrum(x,fs,'persistence', ... 'FrequencyLimits',[100 290],'TimeResolution',1)
The wavelet transform is a linear timefrequency representation that preserves time shifts and time scalings.
The continuous wavelet transform is good at detecting transients in nonstationary signals, and for signals in which instantaneous frequency grows rapidly.
The CWT is invertible.
The CWT tiles the timefrequency plane with variablesized windows. The window automatically widens in time, making it suitable for lowfrequency phenomena, and narrows for high frequency phenomena.
The applications of this timefrequency method include, but are not limited to:
Electrocardiograms (ECG): The most clinically useful information of the ECG signal is found in the time intervals between its consecutive waves and amplitudes defined by its features. The wavelet transform breaks down the ECG signal into scales, making it easier to analyze the ECG signal in different frequency ranges easier to analyze.
Electroencephalogram (EEG): Raw EEG signals suffer from poor spatial resolution, low signaltonoise ratio, and artifacts. Continuous wavelet decomposition of a noisy signal concentrates intrinsic signal information in a few wavelet coefficients having large absolute values without modifying the random distribution of noise. Therefore, denoising can be achieved by thresholding the wavelet coefficients.
Signal demodulation: Demodulate extended binary phase shift keying (EBPSK) using an adaptive wavelet construction method.
Deep learning: The CWT can be used to create timefrequency representations that can be used to train a convolutional neural network. Classify Time Series Using Wavelet Analysis and Deep Learning (Wavelet Toolbox) shows how to classify ECG signals using scalograms and transfer learning.
cwt
(Wavelet Toolbox) computes the continuous wavelet
transform and displays the scalogram. Alternatively, create a CWT filter bank using
cwtfilterbank
(Wavelet Toolbox) and apply the wt
(Wavelet Toolbox)
function. Use this method to run in parallel applications or when computing the
transform for several functions in a loop.
icwt
(Wavelet Toolbox) inverts the continuous wavelet
transform.
Signal Analyzer has a scalogram view to visualize the CWT of a time series.
Load a noisy ECG waveform sampled at 360 Hz.
load ecg
Fs = 360;
Compute the continuous wavelet transform.
cwt(ecg,Fs)
The ECG data is taken from the MITBIH Arrhythmia Database [2].
The WignerVille distribution (WVD) is a quadratic energy density computed by correlating the signal with a time and frequency translated and complexconjugated version of itself.
The WignerVille distribution is always real even if the signal is complex.
Time and frequency marginal densities correspond to instantaneous power and spectral energy density, respectively.
The instantaneous frequency and group delay can be evaluated using local firstorder moments of the Wigner distribution.
The time resolution of the WVD is equal to the number of input samples.
The Wigner distribution can locally assume negative values.
The applications of this timefrequency method include, but are not limited to:
Otoacoustic emissions (OAEs): OAEs are narrowband oscillatory signals emitted by the cochlea (inner ear), and their presence is indicative of normal hearing.
Quantum mechanics: Quantum corrections to classical statistical mechanics, model electron transport, and calculate static and dynamic properties of manybody quantum systems.
wvd
computes the WignerVille distribution.
xwvd
computes the cross WignerVille distribution
of two signals. See Use Cross WignerVille Distribution to Estimate Instantaneous Frequency
for more details.
Load a data file containing otoacoustic emission data sampled at 20 kHz. The emission is produced by a stimulus beginning at 25 milliseconds and ending at 175 milliseconds.
load dpoae
Fs = 20e3;
Compute the smoothedpseudo Wigner Ville distribution of the otoacoustic data. The convenience plot isolates the emission frequency at roughly the expected value 1.2 kHz.
wvd(dpoaets,Fs,'smoothedPseudo',kaiser(511,10),kaiser(511,10),'NumFrequencyPoints',4000,'NumTimePoints',3990)
For more details on otoacoustic emissions, see "Determining Exact Frequency Through the Analytic CWT" in CWTBased TimeFrequency Analysis (Wavelet Toolbox).
Reassignment sharpens the localization of spectral estimates and produces spectrograms that are easier to read and interpret. The technique relocates each spectral estimate to the center of energy of its bin instead of the bin's geometric center. It provides exact localization for chirps and impulses.
The Fourier synchrosqueezed transform starts from the shorttime Fourier transform and "squeezes" its values so that they concentrate around curves of instantaneous frequency in the timefrequency plane.
The wavelet synchrosqueezed transform reassigns the signal energy in frequency.
Both the Fourier synchrosqueezed transform and the wavelet synchrosqueezed transform are invertible.
The reassigned and synchrosqueezing methods are especially suited to track and extract timefrequency ridges.
The applications of this timefrequency method include, but are not limited to:
Audio signal processing: Synchrosqueezing transform (SST) was originally introduced in the context of audio signal analysis.
Seismic data: Analysis of seismic data to find oil and gas traps. Synchrosqueezing can also detect deeplayer weak signals that are usually smeared in seismic data.
Oscillations in power systems: A steam turbine and electric generator can have mechanical subsynchronous oscillation (SSO) modes between the various turbine stages and the generator. The frequency of the SSO is generally between 5 Hz and 45 Hz, and the mode frequencies are often close to each other. The antinoise ability and timefrequency resolution of WSST improves the readability of the timefrequency view.
Deep learning: Synchrosqueezed transforms can be used to
extract timefrequency features and fed into a network that classifies timeseries
data. Waveform Segmentation Using Deep Learning shows how
fsst
outputs can be fed into an LSTM
network that classifies ECG signals.
Use the 'reassigned'
option in spectrogram
, set the 'Reassigned'
argument to
true
in pspectrum
, or check the Reassign box in the
spectrogram view of Signal Analyzer to compute reassigned
spectrograms.
fsst
computes the Fourier synchrosqueezed
transform. Use the ifsst
function to invert the Fourier
synchrosqueezed transform. (See Fourier Synchrosqueezed Transform of Speech Signal for reconstruction of speech signals using
ifsst
.)
wsst
(Wavelet Toolbox) computes the wavelet
synchrosqueezed transform. Use the iwsst
(Wavelet Toolbox) function to invert the wavelet synchrosqueezed transform. (See
Inverse Synchrosqueezed Transform of Chirp (Wavelet Toolbox) for reconstruction of a quadratic chirp using
iwsst
(Wavelet Toolbox).)
Load an echolocation pulse emitted by a big brown bat (Eptesicus Fuscus). The sampling interval is 7 microseconds.
load batsignal
Fs = 1/DT;
Compute the reassigned spectrogram of the signal.
subplot(2,1,1) pspectrum(batsignal,Fs,'spectrogram','TimeResolution',280e6, ... 'OverlapPercent',85,'MinThreshold',45,'Leakage',0.9) subplot(2,1,2) pspectrum(batsignal,Fs,'spectrogram','TimeResolution',280e6, ... 'OverlapPercent',85,'MinThreshold',45,'Leakage',0.9,'Reassign',true)
Thanks to Curtis Condon, Ken White, and Al Feng of the Beckman Center at the University of Illinois for the bat data and permission to use it in this example [3].
Load a file containing the word "strong," spoken by a woman and by a man. The signals are sampled at 8 kHz. Concatenate them into a single signal.
load Strong
x = [her' him'];
Compute the synchrosqueezed Fourier transform of the signal. Window the signal using a Kaiser window with shape factor $\beta =20$.
fsst(x,Fs,kaiser(256,20),'yaxis')
Load the synthetic seismic data sampled at 100 Hz for 1 second.
load SyntheticSeismicData
Compute the wavelet synchrosqueezed transform of the seismic data using the bump wavelet and 30 voices per octave.
wsst(x,Fs,'bump','VoicesPerOctave',30,'ExtendSignal',true)
The seismic signal is generated using the two sinusoids mentioned in "TimeFrequency Analysis of Seismic Data Using Synchrosqueezing Transform" by Ping Wang, Jinghuai Gao, and Zhiguo Wang [4].
Load acceleration measurements recorded on the first floor of a three story test structure under earthquake conditions. The measurements are sampled at 1 kHz.
load quakevib
Fs = 1e3;
Compute the wavelet synchrosqueezed transform of the acceleration measurements. You are analyzing vibration data that exhibit a cyclic behavior. The synchrosqueezed transform allows you to isolate the three frequency components, separated by roughly 11 Hz. The main vibration frequency is at 5.86 Hz, and the equispaced frequency peaks suggest that they are harmonically related. The cyclic behavior of the vibrations is also visible.
wsst(gfloor1OL,Fs,'bump','VoicesPerOctave',48) ylim([0 35])
Load seismograph data recorded during the 1995 Kobe earthquake. The data has a sample rate of 1 Hz.
load kobe
Fs = 1;
Compute the wavelet synchrosqueezed transform that isolates the different frequency components of the seismic data.
wsst(kobe,Fs,'bump','VoicesPerOctave',48) ylim([0 300])
The data are seismograph (vertical acceleration, nm/sq.sec) measurements recorded at Tasmania University, Hobart, Australia on 16 January 1995 beginning at 20:56:51 (GMT) and continuing for 51 minutes at 1 second intervals [5].
Load the subsynchronous oscillation data of a Power System.
load OscillationData
Compute the wavelet synchrosqueezed transform using the bump wavelet and 48 voices per octave. The four mode frequencies are at 15 Hz, 20 Hz, 25 Hz and 32 Hz. Notice that the energies of the modes at 15 Hz and 20 Hz decrease with time, whereas the energy of the modes at 25 Hz and 32 Hz increase gradually over time.
wsst(x,Fs,'bump','VoicesPerOctave',48) ylim([10 50])
This synthetic subsynchronous oscillation data was generated using the equation defined by Zhao et al in "Application of Synchrosqueezed Wavelet Transforms for Extraction of the Oscillatory Parameters of Subsynchronous Oscillation in Power Systems" [6].
The constantQ nonstationary Gabor transform uses windows with different center frequencies and bandwidths such that the ratio of center frequency to bandwidth, the Q factor, remains constant.
The constantQ Gabor transform enables the construction of stable inverses, yielding perfect signal reconstruction.
In frequency space, the windows are centered at logarithmically spaced center frequencies.
The applications of this timefrequency method include, but are not limited to:
Audio signal processing: The fundamental frequencies of the tones in music are geometrically spaced. The frequency resolution of the human auditory system is approximately constantQ, making this technique appropriate for music signal processing.
Load an audio file containing a fragment of Rock music with vocals, drums, and guitar. The signal has a sample rate of 44.1 kHz.
load drums
Set the frequency range over which the CQT has a logarithmic frequency response to be the minimum allowable frequency to 2 kHz. Perform the CQT of the signal using 20 bins per octave.
minFreq = fs/length(audio); maxFreq = 2000; cqt(audio,'SamplingFrequency',fs,'BinsPerOctave',20,'FrequencyLimits',[minFreq maxFreq])
The empirical mode decomposition decomposes the signals into intrinsic mode functions which form a complete and nearly orthogonal basis for the original signal.
The variational mode decomposition decomposes a signal into a small number of narrowband intrinsic mode functions. The method simultaneously calculates all the mode waveforms and their central frequencies by optimizing a constrained variational problem.
The empirical wavelet transform decomposes the signals into multiresolution analysis (MRA) components. The method uses an adaptable wavelet subdivision scheme that automatically determines the empirical wavelet and scaling filters and preserves energy.
The HilbertHuang transform computes the instantaneous frequency of each intrinsic mode function.
The maximal overlap discrete wavelet transform (MODWT) partitions a signal's energy across detail and scaling coefficients. The MODWT is a nondecimated discrete wavelet transform useful for applications that require a shiftinvariant transform. You can obtain multiscale variance and correlation estimates, and invert the transform.
The tunable Qfactor wavelet transform provides a Parseval frame decomposition where energy is partitioned among components, as well as perfect reconstruction of the signal. The tunable Qfactor wavelet transform is a technique that creates an MRA with a userspecified Qfactor. The Qfactor is the ratio of the center frequency to the bandwidth of the filters used in the transform.
These methods combined are useful for analyzing nonlinear and nonstationary signals.
The applications of this timefrequency method include, but are not limited to:
Physiological signal processing: Analyze human EEG response to transcranial magnetic stimulation (TMS) of the brain cortex.
Structural applications: Locate anomalies that appear as cracks, delamination, or stiffness loss in beams and plates.
System identification: Isolate modal damping ratios of structures with closely spaced modal frequencies.
Ocean engineering: Identify transient electromagnetic disturbances caused by humans in underwater electromagnetic environments.
Solar physics: Extract periodic components of sunspot data.
Atmospheric turbulence: Observe stable boundary layer to separate turbulent and nonturbulent motions.
Epidemiology: Assess traveling speed of communicative diseases such as Dengue fever.
emd
computes the empirical mode
decomposition.
vmd
computes the variational mode
decomposition.
ewt
(Wavelet Toolbox)
computes the empirical wavelet transform.
hht
computes the Hilbert Huang spectrum of an
empirical mode decomposition.
modwt
(Wavelet Toolbox) computes the maximal overlap
discrete wavelet transform. To obtain the MRA analysis, use modwtmra
(Wavelet Toolbox).
tqwt
(Wavelet Toolbox)
computes the tunable Qfactor wavelet transform. To obtain the MRA analysis, use
tqwtmra
(Wavelet Toolbox).
Load the vibration signal from a defective bearing generated in the Compute Hilbert Spectrum of Vibration Signal example. The signal is sampled at a rate 10 kHz.
load bearingVibration
Compute the first five intrinsic mode functions (IMFs) of the signal. Plot the Hilbert spectrum of the first and third empirical modes. The first mode reveals increasing wear due to highfrequency impacts on the bearing's outer race. The third mode shows a resonance occurring halfway through the measurement process that caused the defect in the bearing.
imf = emd(y,'MaxNumIMF',5,'Display',0); subplot(2,1,1) hht(imf(:,1),fs) subplot(2,1,2) hht(imf(:,3),fs,'FrequencyLimits',[0 100])
[1] The Pacific blue whale file is obtained from the library of animal vocalizations maintained by the Cornell University Bioacoustics Research Program.
[2] Moody G. B, Mark R. G. The impact of the MITBIH Arrhythmia Database. IEEE Eng in Med and Biol 20(3):4550 (MayJune 2001). (PMID: 11446209)
[3] Thanks to Curtis Condon, Ken White, and Al Feng of the Beckman Center at the University of Illinois for the bat echolocation data.
[4] Wang, Ping, Gao, J., and Wang, Z. TimeFrequency Analysis of Seismic Data Using Synchrosqueezing Transform, IEEE Geoscience and Remote Sensing Letters, Vol 12, Issue 11, Dec. 2014.
[5] Seismograph (vertical acceleration, nm/sq.sec) of the Kobe earthquake, recorded at Tasmania University, Hobart, Australia on 16 January 1995 beginning at 20:56:51 (GMTRUE) and continuing for 51 minutes at 1 second intervals.
[6] Zhao et al. Application of Synchrosqueezed Wavelet Transforms for Extraction of the Oscillatory Parameters of Subsynchronous Oscillation in Power Systems MDPI Energies; Published 12 June 2018.
[7] Boashash, Boualem. TimeFrequency Signal Analysis and Processing: A Comprehensive Reference Elsevier, 2016.
cqt
(Wavelet Toolbox)  cwt
(Wavelet Toolbox)  cwtfilterbank
(Wavelet Toolbox)  dlstft
 emd
 ewt
(Wavelet Toolbox)  fsst
 hht
 icqt
(Wavelet Toolbox)  icwt
(Wavelet Toolbox)  ifsst
 istft
 iwsst
(Wavelet Toolbox)  kurtogram
 modwt
(Wavelet Toolbox)  modwtmra
(Wavelet Toolbox)  pkurtosis
 pspectrum
 spectrogram
 stft
 tqwt
(Wavelet Toolbox)  tqwtmra
(Wavelet Toolbox)  vmd
 wsst
(Wavelet Toolbox)  wt
(Wavelet Toolbox)  xspectrogram
 wvd
 xwvd