TimeFrequency Gallery
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  


ShortTime Fourier Transform (Spectrogram)
Description
The shorttime Fourier transform is a linear timefrequency representation useful in the analysis of nonstationary multicomponent signals.
The spectrogram is the magnitude squared of the STFT. For more information about computing the spectrogram, see Spectrogram Computation with Signal Processing Toolbox (Signal Processing Toolbox).
The shorttime Fourier transform is invertible.
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.
Potential Applications
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 (Signal Processing Toolbox) 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.
How to Use
stft
(Signal Processing Toolbox) computes the shorttime Fourier transform. To invert the shorttime Fourier transform, use theistft
(Signal Processing Toolbox) function.dlstft
(Signal Processing Toolbox) computes the deep learning shorttime Fourier transform. You must have Deep Learning Toolbox™ installed.stftLayer
(Signal Processing Toolbox) computes the shorttime Fourier transform within a deep learning network. You must have Deep Learning Toolbox installed.stftmag2sig
(Signal Processing Toolbox) reconstructs a signal starting with the magnitude of its timefrequency representation. You must have Deep Learning Toolbox installed.pspectrum
(Signal Processing Toolbox) orspectrogram
(Signal Processing Toolbox) computes the spectrogram.xspectrogram
(Signal Processing Toolbox) computes the crossspectrogram of two signals.You can also use the spectrogram view in Signal Analyzer (Signal Processing Toolbox) to view the spectrogram of a signal.
Use the persistence spectrum option in
pspectrum
(Signal Processing Toolbox) or Signal Analyzer (Signal Processing Toolbox) to identify signals hidden in other signals.
Example: Pulses and Oscillations
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))
Example: Audio Signal with Decreasing Chirps
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)
Example: Whale Song
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)
Example: Persistence Spectrum of Transient Signal
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)
Continuous Wavelet Transform (Scalogram)
Description
The wavelet transform is a linear timefrequency representation that preserves time shifts and time scalings.
The continuous wavelet transform (CWT) 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.
Potential Applications
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 shows how to classify ECG signals using scalograms and transfer learning.
How to Use
cwt
computes the continuous wavelet transform and displays the scalogram. Alternatively, create a CWT filter bank usingcwtfilterbank
and apply thewt
function. Use this method to run in parallel applications or when computing the transform for several functions in a loop. To invert the continuous wavelet transform, use theicwt
function.dlcwt
computes the deep learning CWT. You must have Deep Learning Toolbox installed.cwtLayer
computes the CWT within a deep learning network. You must have Deep Learning Toolbox installed.You can use Wavelet TimeFrequency Analyzer to compute and visualize the scalogram of a time series. You can also generate a script to reproduce the CWT analysis in your workspace.
Signal Analyzer (Signal Processing Toolbox) has a scalogram view to visualize the CWT of a time series.
Example: ECG Signal
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].
Example: Signal with Transient Events
Generate a signal sampled at 1 kHz for 1.5 seconds. The signal consists of two sinusoidal components that turn on and off abruptly, with additional "defects" at 222 milliseconds and 800 milliseconds.
Fs = 1e3; t = 0:1/Fs:1.51/Fs; addNoise = 0.025*randn(size(t)); x = cos(2*pi*150*t).*(t>=0.1 & t<0.5) + ... sin(2*pi*200*t).*(t>0.7 & t<=1.2); x = x+addNoise; x([222 800]) = x([222 800])+[2 2]; plot(t.*1000,x) xlabel("Time (ms)") ylabel("Amplitude")
Obtain the CWT of the signal. Use the analytic Morlet wavelet. Plot only the finestscale wavelet coefficients. Create a second axis and plot the original time data.
cfs = cwt(x,Fs,"amor"); plot(t,abs(cfs(1,:)),LineWidth=2) ylabel("Magnitude") set(gca,XTick=[0.1 0.222 0.5 0.7 0.8 1.2]) yyaxis right plot(t,x,"") ylabel("Amplitude") xlabel("Time(s)")
WignerVille Distribution
Description
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.
Potential Applications
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.
How to Use
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.
Example: Otoacoustic Emission
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.
Reassignment and Synchrosqueezing
Description
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.
Potential Applications
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 (Signal Processing Toolbox) shows how
fsst
(Signal Processing Toolbox) outputs can be fed into an LSTM network that classifies ECG signals.
How to Use
Use the
'reassigned'
option inspectrogram
(Signal Processing Toolbox), set the'Reassigned'
argument totrue
inpspectrum
(Signal Processing Toolbox), or check the Reassign box in the spectrogram view of Signal Analyzer (Signal Processing Toolbox) to compute reassigned spectrograms.fsst
(Signal Processing Toolbox) computes the Fourier synchrosqueezed transform. Use theifsst
(Signal Processing Toolbox) function to invert the Fourier synchrosqueezed transform. (See Fourier Synchrosqueezed Transform of Speech Signal (Signal Processing Toolbox) for reconstruction of speech signals usingifsst
(Signal Processing Toolbox).)wsst
computes the wavelet synchrosqueezed transform. Use theiwsst
function to invert the wavelet synchrosqueezed transform. (See Inverse Synchrosqueezed Transform of Chirp for reconstruction of a quadratic chirp usingiwsst
.)
Example: Echolocation Pulse
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.
tiledlayout(2,1) nexttile pspectrum(batsignal,Fs,"spectrogram", ... TimeResolution=280e6, ... OverlapPercent=85, ... MinThreshold=45, ... Leakage=0.9) nexttile 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].
Example: Speech Signals
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")
Example: Synthetic Seismic Data
Load the synthetic seismic data sampled at 100 Hz for 1 second.
load SyntheticSeismicData
The workspace now contains the synthetic data, x
, and a sample rate, Fs
.
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].
Example: Earthquake Vibration
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])
Example: Kobe Earthquake Data
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].
Example: Subsynchronous Oscillation in Power Systems
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].
ConstantQ Gabor Transform
Description
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.
Potential Applications
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.
How to Use
Example: Rock Music
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
The workspace now contains a vector of audio data, audio
, and a sample rate, fs
.
Set the frequency range over which the constantQ Gabor transform has a logarithmic frequency response to be the minimum allowable frequency to 2 kHz. Compute 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])
DataAdaptive Methods and Multiresolution Analysis
Description
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.
Potential Applications
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.
How to Use
emd
computes the empirical mode decomposition.vmd
computes the variational mode decomposition.ewt
computes the empirical wavelet transform.hht
computes the Hilbert Huang spectrum of an empirical mode decomposition.modwpt
computes the maximal overlap discrete wavelet packet transform. To invert the MODWPT transform, use theimodwpt
function. To obtain the MODWPT details, usemodwptdetails
. The MODWPT details provide zerophase filtering of the signal.modwt
computes the maximal overlap discrete wavelet transform. To invert the MODWT, use theimodwt
function. To obtain the MRA analysis, usemodwtmra
.tqwt
computes the tunable Qfactor wavelet transform. To invert the TQWT, use theitqwt
function. To obtain the MRA analysis, usetqwtmra
.You can use Signal Multiresolution Analyzer to visualize multilevel wavelet and data adaptivebased decompositions of realvalued 1D signals.
Wavelet Signal Analyzer enables visualization, analysis, and compression of 1D signals using decimated and nondecimated discrete wavelet and wavelet packet transforms.
Example: Bearing Vibration
Load the vibration signal from a defective bearing generated in the Compute Hilbert Spectrum of Vibration Signal (Signal Processing Toolbox) 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); tiledlayout(2,1) nexttile hht(imf(:,1),fs) nexttile hht(imf(:,3),fs, ... FrequencyLimits=[0 100])
References
[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.
See Also
Apps
 Signal Analyzer (Signal Processing Toolbox)  Signal Multiresolution Analyzer  Wavelet Signal Analyzer  Wavelet TimeFrequency Analyzer
Functions
cqt
cwt
cwtfilterbank
cwtLayer
dlcwt
dlstft
(Signal Processing Toolbox) emd
ewt
fsst
(Signal Processing Toolbox) hht
icqt
icwt
ifsst
(Signal Processing Toolbox) imodwpt
imodwt
istft
(Signal Processing Toolbox) itqwt
iwsst
kurtogram
(Signal Processing Toolbox) modwpt
modwptdetails
modwt
modwtmra
pkurtosis
(Signal Processing Toolbox) pspectrum
(Signal Processing Toolbox) spectrogram
(Signal Processing Toolbox) stft
(Signal Processing Toolbox) stftLayer
(Signal Processing Toolbox) tqwt
tqwtmra
vmd
wsst
wt
xspectrogram
(Signal Processing Toolbox) wvd
xwvd