## 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.

MethodFeaturesInvertibleExamples
• The short-time Fourier transform (STFT) has fixed time-frequency resolution.

• The spectrogram is the magnitude squared of the STFT.

• `stft`: Yes

• `spectrogram`: No

• The continuous wavelet transform (CWT) has a variable time-frequency resolution.

• The CWT preserves time shifts and time scalings.

Yes

• The Wigner-Ville distribution (WVD) is always real.

• Time and frequency marginal densities correspond to instantaneous power and spectral energy density, respectively.

• The time resolution of the WVD is equal to the number of input samples.

No

• Reassignment sharpens localization of spectral estimates.

• Synchrosqueezing "condenses" time-frequency maps around curves of instantaneous frequency.

• Both methods are especially suited to track and extract time-frequency ridges

• `pspectrum`: No

• `fsst`, `wsst`: Yes

• The constant-Q Gabor transform (CQT) tiles the time-frequency plane with variable-sized windows.

• The windows have adaptable bandwidth and sampling density.

• The ratio of center frequency to bandwidth (Q-factor) for all windows is constant.

Yes

• The empirical mode decomposition (EMD) decomposes signals into intrinsic mode functions.

• The variational mode decomposition (VMD) decomposes a signal into a small number of narrowband intrinsic mode functions.

• The empirical wavelet transform (EWT) decomposes signals into multiresolution analysis (MRA) components.

• The Hilbert-Huang transform (HHT) computes the instantaneous frequency of each empirical mode.

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 short-time Fourier transform is invertible.

• The spectrogram is the magnitude squared of the STFT.

• 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 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.

#### 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 $\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 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.

```whaleFile = fullfile(matlabroot,'examples','matlab','data','bluewhale.au'); [w,fs] = audioread(whaleFile);```

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 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. Classify Time Series Using Wavelet Analysis and Deep Learning (Wavelet Toolbox) shows how to classify ECG signals using scalograms and transfer learning.

#### How to Use

• `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.

#### 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.

#### 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 CWT-Based Time-Frequency Analysis (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. Waveform Segmentation Using Deep Learning shows how `fsst` outputs can be fed into an LSTM network that classifies ECG signals.

#### 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 $\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`

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.

#### 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])```

#### 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.

• 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.

#### Example: Bearing Vibration

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 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.