time-frequency gallery
this gallery provides you with an overview of the time-frequency 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 |
short-time fourier transform (spectrogram)
description
the short-time fourier transform is a linear time-frequency 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 .
the short-time fourier transform is invertible.
you can compute the cross-spectrogram of two signals to look for similarities in time-frequency space.
the persistence spectrum of a signal is a time-frequency view that shows the percentage of the time that a given frequency is present in a signal. the persistence spectrum is a histogram in power-frequency 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 time-frequency method include, but are not limited to:
audio signal processing: fundamental frequency estimation, cross synthesis, spectral envelope extraction, time-scale modification, time-stretching, and pitch shifting. (see 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
computes the short-time fourier transform. to invert the short-time fourier transform, use the function.
computes the deep learning short-time fourier transform. you must have deep learning toolbox™ installed.
or computes the spectrogram.
computes the cross-spectrogram 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 or signal analyzer 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:4-1/fs; x = 10*besselj(0,1000*(sin(2*pi*(t 2).^3/60).^5));
compute and plot the short-time fourier transform of the signal. window the signal with a 200-sample kaiser window with shape factor .
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 short-time 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 time-frequency 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 time-frequency plane with variable-sized windows. the window automatically widens in time, making it suitable for low-frequency phenomena, and narrows for high frequency phenomena.
potential applications
the applications of this time-frequency 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 signal-to-noise 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 time-frequency representations that can be used to train a convolutional neural network. (wavelet toolbox) shows how to classify ecg signals using scalograms and transfer learning.
how to use
(wavelet toolbox) computes the continuous wavelet transform and displays the scalogram. alternatively, create a cwt filter bank using (wavelet toolbox) and apply the (wavelet toolbox) function. use this method to run in parallel applications or when computing the transform for several functions in a loop.
(wavelet toolbox) inverts the continuous wavelet transform.
signal analyzer 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 mit-bih arrhythmia database [2].
wigner-ville distribution
description
the wigner-ville distribution (wvd) is a quadratic energy density computed by correlating the signal with a time and frequency translated and complex-conjugated version of itself.
the wigner-ville 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 first-order 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 time-frequency 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 many-body quantum systems.
how to use
computes the wigner-ville distribution.
computes the cross wigner-ville distribution of two signals. see 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 smoothed-pseudo 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 (wavelet toolbox).
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 short-time fourier transform and "squeezes" its values so that they concentrate around curves of instantaneous frequency in the time-frequency 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 time-frequency ridges.
potential applications
the applications of this time-frequency 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 deep-layer 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 time-frequency resolution of wsst improves the readability of the time-frequency view.
deep learning: synchrosqueezed transforms can be used to extract time-frequency features and fed into a network that classifies time-series data. shows how outputs can be fed into an lstm network that classifies ecg signals.
how to use
use the
'reassigned'
option in , set the'reassigned'
argument totrue
in , or check the reassign box in the spectrogram view of signal analyzer to compute reassigned spectrograms.computes the fourier synchrosqueezed transform. use the function to invert the fourier synchrosqueezed transform. (see for reconstruction of speech signals using .)
(wavelet toolbox) computes the wavelet synchrosqueezed transform. use the (wavelet toolbox) function to invert the wavelet synchrosqueezed transform. (see (wavelet toolbox) for reconstruction of a quadratic chirp using (wavelet toolbox).)
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.
subplot(2,1,1) pspectrum(batsignal,fs,'spectrogram','timeresolution',280e-6, ... 'overlappercent',85,'minthreshold',-45,'leakage',0.9) subplot(2,1,2) pspectrum(batsignal,fs,'spectrogram','timeresolution',280e-6, ... '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 .
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
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 "time-frequency 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].
constant-q gabor transform
description
the constant-q 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 constant-q 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 time-frequency 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 constant-q, making this technique appropriate for music signal processing.
how to use
(wavelet toolbox) computes the constant-q gabor transform.
(wavelet toolbox) inverts the constant-q gabor transform.
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
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])
data-adaptive 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 hilbert-huang 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 shift-invariant transform. you can obtain multiscale variance and correlation estimates, and invert the transform.
the tunable q-factor wavelet transform provides a parseval frame decomposition where energy is partitioned among components, as well as perfect reconstruction of the signal. the tunable q-factor wavelet transform is a technique that creates an mra with a user-specified q-factor. the q-factor 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 time-frequency 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
computes the empirical mode decomposition.
computes the variational mode decomposition.
(wavelet toolbox) computes the empirical wavelet transform.
computes the hilbert huang spectrum of an empirical mode decomposition.
(wavelet toolbox) computes the maximal overlap discrete wavelet transform. to obtain the mra analysis, use (wavelet toolbox).
(wavelet toolbox) computes the tunable q-factor wavelet transform. to obtain the mra analysis, use (wavelet toolbox).
example: bearing vibration
load the vibration signal from a defective bearing generated in the 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 high-frequency 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])
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 mit-bih arrhythmia database. ieee eng in med and biol 20(3):45-50 (may-june 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. time-frequency 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. time-frequency signal analysis and processing: a comprehensive reference elsevier, 2016.
see also
apps
- signal analyzer | (wavelet toolbox)
functions
- (wavelet toolbox) | (wavelet toolbox) | (wavelet toolbox) | | | (wavelet toolbox) | | | (wavelet toolbox) | (wavelet toolbox) | | | (wavelet toolbox) | | (wavelet toolbox) | (wavelet toolbox) | | | | | (wavelet toolbox) | (wavelet toolbox) | | (wavelet toolbox) | (wavelet toolbox) | | |