# stftmag2sig

## Syntax

## Description

returns a reconstructed time-domain real signal, `x`

= stftmag2sig(`s`

,`nfft`

)`x`

, estimated from
the Short-Time Fourier Transform (STFT) magnitude,
`s`

, based on the Griffin-Lim algorithm. The function assumes
`s`

was computed using discrete Fourier transform (DFT) length
`nfft`

.

specifies additional options using name-value arguments. Options include, among others,
the FFT window and the method to specify initial phases. These arguments can be added to
any of the previous input syntaxes. For example,
`x`

= stftmag2sig(___,`Name=Value`

)`FrequencyRange="onesided",InitializePhaseMethod="random"`

specifies
that the signal is reconstructed from a one-sided STFT with random initial phases.

## Examples

### Reconstruct Sinusoid from STFT Magnitude

Consider 512 samples of a sinusoid with a normalized frequency of $\pi /60$ rad/sample and a DC value of 1. Compute the STFT of the signal.

n = 512; x = cos(pi/60*(0:n-1)')+1; S = stft(x);

Reconstruct the sinusoid from the magnitude of the STFT. Plot the original and reconstructed signals.

xr = stftmag2sig(abs(S),size(S,1)); plot(x) hold on plot(xr,"--",LineWidth=2) hold off legend("Original","Reconstructed")

Repeat the computation, but now pad the signal with zeros to decrease edge effects.

xz = circshift([x; zeros(n,1)],n/2); Sz = stft(xz); xr = stftmag2sig(abs(Sz),size(Sz,1)); xz = xz(n/2+(1:n)); xr = xr(n/2+(1:n)); plot(xz) hold on plot(xr,"--",LineWidth=2) hold off legend("Original","Reconstructed")

Repeat the computation, but now decrease edge effects by assuming that `x`

is a segment of a signal twice as long.

xx = cos(pi/60*(-n/2:n/2+n-1)')+1; Sx = stft(xx); xr = stftmag2sig(abs(Sx),size(Sx,1)); xx = xx(n/2+(1:n)); xr = xr(n/2+(1:n)); plot(xx) hold on plot(xr,"--",LineWidth=2) hold off legend("Original","Reconstructed")

### Reconstruct Audio Signal from STFT Magnitude

Load an audio signal that contains two decreasing chirps and a wideband splatter sound. The signal is sampled at 8192 Hz. Plot the STFT of the signal. Divide the waveform into 128-sample segments and window the segments using a Hamming window. Specify 64 samples of overlap between adjoining segments and 1024 DFT points.

load splat ty = (0:length(y)-1)/Fs; % To hear, type sound(y,Fs) wind = hamming(128); olen = 64; nfft = 1024; stft(y,Fs,Window=wind,OverlapLength=olen,FFTLength=nfft)

Compute the magnitude and phase of the STFT.

s = stft(y,Fs,Window=wind,OverlapLength=olen,FFTLength=nfft); smag = abs(s); sphs = angle(s);

Reconstruct the signal based on the magnitude of the STFT. Use the same parameters that you used to compute the STFT. By default, `stftmag2sig`

initializes the phases to zero and uses 100 optimization iterations.

[x,tx,info] = stftmag2sig(smag,nfft,Fs, ... Window=wind,OverlapLength=olen); % To hear the reconstruction, type sound(x,Fs)

Plot the original and reconstructed signals. For better comparison, shift the reconstructed signal vertically and to the right.

plot(ty,y,tx+500/Fs,x+1) legend("Original","Reconstructed",Location="best")

Output the relative improvement toward convergence between the last two iterations.

impr = info.Inconsistency

impr = 0.0579

Improve the reconstruction by doubling the number of optimization iterations and setting the initial phases to the actual phases from the STFT. Plot the original and reconstructed signals. For better comparison, plot the negative of the reconstructed signal and shift it vertically and to the right.

[x,tx,info] = stftmag2sig(smag,nfft,Fs, ... Window=wind,OverlapLength=olen, ... MaxIterations=200,InitialPhase=sphs); % To hear the reconstruction, type sound(x,Fs) plot(ty,y,tx+500/Fs,-x+1) legend("Original","Reconstructed",Location="best")

Output the relative improvement toward convergence between the last two iterations.

impr = info.Inconsistency

impr = 2.0848e-16

### Reconstruct Pulsating Signal from STFT Magnitude Using Gradient Descent

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. Divide the signal into 128-sample segments and window the segments using a Hann window. Specify 96 samples of overlap between adjoining segments and 128 DFT points.

win = hann(128); olen = 96; nfft = 128; stft(x,fs,Window=win,OverlapLength=olen,FFTLength=nfft)

Compute the magnitude and phase of the STFT.

s = stft(x,fs,Window=win,OverlapLength=olen,FFTLength=nfft); smag = abs(s); sphs = angle(s);

Reconstruct the signal based on the magnitude of the STFT. Use the same parameters that you used to compute the STFT. Instead of using the default Griffin-Lim algorithm, use the gradient descent method with the ADAM optimizer.

[xrec,trec,info] = stftmag2sig(smag,nfft,fs, ... Window=win,Method="gd",Optimizer="adam");

Plot the original and reconstructed signals. For better comparison, shift the reconstructed signal vertically so both are visible.

plot(t,x,trec,xrec+12) legend("Original","Reconstructed",Location="northwest") ylim([-5 27])

Compute the mean-squared error between the original signal and the reconstruction.

sum((x-xrec').^2)/length(x)

ans = 1.3620

### Reconstruct Multichannel Signal from STFT Magnitude

For purposes of reproducibility, set the random seed to the default value. Generate 160 samples of a two-channel sinusoid. Add noise to each channel.

The first channel has unit amplitude and a normalized sinusoid frequency of $$\pi /4$$ rad/sample.

The second channel has unit amplitude and a normalized sinusoid frequency of $$\pi /2$$ rad/sample.

Compute the STFT of the signal. Divide the signal into 32-sample segments and window the segments using a Hann window.

```
rng default
N = 160;
t = 0:N-1;
win = hann(32);
signal = cos(pi./[4;2]*t)'+randn(N,2)/5;
s = stft(signal,Window=win);
```

Reconstruct the signal from the magnitude of the STFT. Use the gradient descent method with the L-BFGS optimizer. The L-BFGS optimizer automatically selects the step size for each iteration. Track the normalized inconsistency of each channel during the optimization process.

[x,tx,info] = stftmag2sig(abs(s), ... size(s,1),Window=win, ... Method="gd",Optimizer="lbfgs", ... Display=true);

#Iteration | Normalized Inconsistency 1 | 1.2078e+00 1.2042e+00 20 | 3.2348e-02 9.6221e-03 40 | 1.9969e-02 6.8350e-03 60 | 2.5463e-03 1.2099e-03 80 | 2.4790e-03 8.4016e-04 100 | 5.4774e-03 7.1804e-04 Decomposition process stopped. The number of iterations reached the specified "MaxIterations" value of 100.

Plot the original and reconstructed signals.

tiledlayout(2,1) nexttile plot(t,signal(:,1)) hold on plot(tx,x(:,1),"--",LineWidth=2) hold off title("First Channel") legend("Original","Reconstructed", ... Location="northoutside") nexttile plot(t,signal(:,2)) hold on plot(tx,x(:,2),"--",LineWidth=2) hold off title("Second Channel")

Output the relative improvement toward convergence between the last two iterations.

info.Inconsistency

`ans = `*1×2*
0.0055 0.0007

## Input Arguments

`s`

— STFT magnitude

matrix | 3-D array

STFT magnitude, specified as a matrix or a 3-D array. `s`

must
correspond to a real-valued signal.

When

`s`

corresponds to a single-channel signal, specify`s`

as a matrix.When

`s`

corresponds to a multichannel signal, specify`s`

as a 3-D array, where the third dimension corresponds to the channels.

**Example: **

specifies the STFT magnitude of a sinusoid.`abs`

(`stft`

(`sin`

(pi/2*(0:255)),FFTLength=128))

**Example: **

specifies the STFT magnitude of a chirp sampled at 1 kHz.`abs`

(`stft`

(`chirp`

(0:1/1e3:1,25,1,50)))

**Example: **

specifies the STFT magnitude of a two-channel sinusoid with noise.`abs`

(`stft`

(`cos`

(pi./[4;2]*(0:159))'+`randn`

(160,2)/5))

**Data Types: **`single`

| `double`

`nfft`

— Number of DFT points

positive integer scalar

Number of DFT points, specified as a positive integer scalar. This argument is always required.

**Data Types: **`single`

| `double`

`fs`

— Sample rate

2*π* (default) | positive numeric scalar

Sample rate, specified as a positive numeric scalar.

`ts`

— Sample time

`duration`

scalar

Sample time, specified as a `duration`

scalar. Specifying `ts`

is equivalent to
setting a sample rate *f*_{s} =
1/`ts`

.

**Example: **`seconds(1)`

is a `duration`

scalar
representing a 1-second time difference between consecutive signal
samples.

**Data Types: **`duration`

### Name-Value Arguments

Specify optional pairs of arguments as
`Name1=Value1,...,NameN=ValueN`

, where `Name`

is
the argument name and `Value`

is the corresponding value.
Name-value arguments must appear after other arguments, but the order of the
pairs does not matter.

**Example: **`FrequencyRange="onesided",InitializePhaseMethod="random",Method="gd"`

specifies that the signal is reconstructed from a one-sided STFT with random initial phases
using the gradient descent method.

*
Before R2021a, use commas to separate each name and value, and enclose*
`Name`

*in quotes.*

**Example: **`"FrequencyRange","onesided","InitializePhaseMethod","random"`

specifies that the signal is reconstructed from a one-sided STFT with random initial
phases.

`Display`

— Inconsistency display option

`false`

or `0`

(default) | `true`

or `1`

Inconsistency display option, specified as a numeric or logical
`1`

(`true`

) or `0`

(`false`

). If this option is set to `true`

,
`stftmag2sig`

displays the normalized inconsistency after every
20 optimization iterations. The function also displays stopping information at the end
of the run.

**Data Types: **`logical`

`FrequencyRange`

— Frequency range of STFT magnitude

`"centered"`

(default) | `"twosided"`

| `"onesided"`

Frequency range of STFT magnitude, specified as one of these:

`"centered"`

— Treat`s`

as the magnitude of a two-sided, centered STFT. If`nfft`

is even, then`s`

is taken to have been computed over the interval (–*π*,*π*] rad/sample. If`nfft`

is odd, then`s`

is taken to have been computed over the interval (–*π*,*π*) rad/sample. If you specify time information, then the intervals are (–*f*_{s},*f*_{s}/2] cycles/unit time and (–*f*_{s},*f*_{s}/2) cycles/unit time, respectively, where*f*_{s}is the sample rate.`"twosided"`

— Treat`s`

as the magnitude of a two-sided STFT computed over the interval [0, 2*π*) rad/sample. If you specify time information, then the interval is [0,*f*_{s}) cycles/unit time.`"onesided"`

— Treat`s`

as the magnitude of a one-sided STFT. If`nfft`

is even, then`s`

is taken to have been computed over the interval [0,*π*] rad/sample. If`nfft`

is odd, then`s`

is taken to have been computed over the interval [0,*π*) rad/sample. If you specify time information, then the intervals are [0,*f*_{s}/2] cycles/unit time and [0,*f*_{s}/2) cycles/unit time, respectively, where*f*_{s}is the sample rate.

**Data Types: **`char`

| `string`

`InconsistencyTolerance`

— Inconsistency tolerance of reconstruction process

`1e-4`

(default) | positive scalar

Inconsistency tolerance of reconstruction process, specified as a positive scalar. The reconstruction process stops when the Normalized Inconsistency is lower than the tolerance.

**Data Types: **`single`

| `double`

`InitializePhaseMethod`

— Phase initialization

`"zeros"`

(default) | `"random"`

Phase initialization, specified as `"zeros"`

or
`"random"`

. Specify only one of
`InitializePhaseMethod`

or `InitialPhase`

.

`"zeros"`

— The function initializes the phases as zeros.`"random"`

— The function initializes the phases as random numbers distributed uniformly in the interval [–*π*,*π*].

**Data Types: **`char`

| `string`

`InitialPhase`

— Initial phases

matrix | 3-D array

Initial phases, specified as a real numeric matrix or 3-D array. Elements of
`InitialPhase`

must be in the range [–*π*, *π*]. `InitialPhase`

must have the same size as
`s`

. Specify only one of
`InitializePhaseMethod`

or
`InitialPhase`

.

**Example: **`angle(`

specifies the phases of the short-time Fourier transform of a random
signal.`stft`

(`randn`

(1000,1)))

**Example: **`2*pi*(rand(size(stft(randn(1000,1))))-1/2)`

specifies a
matrix of random phases distributed uniformly in the interval [–*π*, *π*]. The matrix has the same size as the short-time Fourier transform of
a random signal.

.

**Data Types: **`single`

| `double`

`InputTimeDimension`

— Input time dimension

`"acrosscolumns"`

(default) | `"downrows"`

Input time dimension, specified as one of these:

`"acrosscolumns"`

— The function assumes that the time dimension of`s`

is across the columns and the frequency dimension is down the rows.`"downrows"`

— The function assumes that the time dimension of`s`

is down the rows and the frequency dimension is across the columns.

**Data Types: **`char`

| `string`

`MaxIterations`

— Maximum number of optimization iterations

`100`

(default) | positive integer scalar

Maximum number of optimization iterations, specified as a positive integer scalar.
The reconstruction process stops when the number of iterations is greater than
`MaxIterations`

.

**Data Types: **`single`

| `double`

`Method`

— Signal reconstruction algorithm

`"gla"`

(default) | `"fgla"`

| `"legla"`

| `"gd"`

Signal reconstruction algorithm, specified as one of these:

`"gla"`

— The original reconstruction algorithm, proposed by Griffin and Lim and described in [1].`"fgla"`

— A fast Griffin-Lim algorithm proposed by Perraudin, Balazs, and Søndergaard and described in [2].`"legla"`

— A fast algorithm proposed by Le Roux, Kameoka, Ono, and Sagayama and described in [3].`"gd"`

— A gradient descent method, proposed by Ji and Tie and described in [4].*(since R2023b)*

**Note**

To specify the gradient descent method, you must have a Deep Learning Toolbox™ license.

**Data Types: **`char`

| `string`

`Optimizer`

— Optimizer used in the gradient descent method

`"sgdm"`

(default) | `"adam"`

| `"rmsprop"`

| `"lbfgs"`

*Since R2023b*

Optimizer to use in the gradient descent method, specified as one of these:

`"sgdm"`

— Stochastic gradient descent with momentum (SGDM) optimizer`"adam"`

— Adaptive moment estimation (ADAM) optimizer`"rmsprop"`

— Root-mean-square propagation (RMSProp) optimizer`"lbfgs"`

— Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimizer

This argument applies only when `Method`

is set to
`"gd"`

.

**Data Types: **`char`

| `string`

`OverlapLength`

— Number of overlapped samples

`75%`

of window length (default) | nonnegative integer

Number of overlapped samples between adjoining segments, specified as a positive
integer smaller than the length of `Window`

. Successful signal
reconstruction requires `OverlapLength`

to match the number of
overlapped segments used to generate the STFT magnitude. If you omit
`OverlapLength`

or specify it as empty, it is set to the largest
integer less than or equal to 75% of the window length, which is 96 samples for the
default Hann window.

**Data Types: **`double`

| `single`

`StepSize`

— Step size for iterative updates

`0.01`

(default) | positive scalar

`TruncationOrder`

— Truncation order for `"legla"`

update rule

positive integer

Truncation order for `"legla"`

update rule, specified as a
positive integer. This argument applies only when `Method`

is set
to `"legla"`

and controls the number of phase values updated in each
iteration of that method. If not specified, `TruncationOrder"`

is
determined using an adaptive algorithm.

**Data Types: **`single`

| `double`

`UpdateParameter`

— Update parameter for fast Griffin-Lim algorithm

`0.99`

(default) | positive scalar

Update parameter for the fast Griffin-Lim algorithm, specified as a positive
scalar. This argument applies only when `Method`

is set to
`"fgla"`

and specifies the parameter for that method's update
rule.

**Data Types: **`single`

| `double`

`Window`

— Spectral window

`hann(128,"periodic")`

(default) | vector

Spectral window, specified as a vector. Successful signal reconstruction requires
`Window`

to match the window used to generate the STFT magnitude.
If you do not specify the window or specify it as empty, the function uses a periodic
Hann window of length 128. The length of `Window`

must be greater
than or equal to 2.

For a list of available windows, see Windows.

**Example: **`hann(128,"periodic")`

and
`(1-cos(2*pi*(128:-1:1)'/128))/2`

both specify the default window
used by `stftmag2sig`

.

**Data Types: **`double`

| `single`

## Output Arguments

`x`

— Reconstructed time-domain signal

vector | matrix

Reconstructed time-domain signal, returned as a vector or a matrix.

**Note**

If you want `x`

and `s`

to be the same
length, the value of
`(length(`

must be an integer. Use `x`

)-noverlap)/(length(window)-noverlap)`Window`

to specify the length of `window`

and `OverlapLength`

to specify `noverlap`

.

`t`

— Time instants

vector

Time instants at which the signal is reconstructed, returned as a vector.

`info`

— Reconstruction process information

structure

Reconstruction process information, returned as a structure containing these fields:

`ExitFlag`

— Termination flag.A value of

`0`

indicates the algorithm stopped when it reached the maximum number of iterations.A value of

`1`

indicates the algorithm stopped when it met the relative tolerance.

`NumIterations`

— Total number of iterations.`Inconsistency`

— Average relative improvement toward convergence between the last two iterations.`ReconstructedPhase`

— Reconstructed phase at the last iteration.`ReconstructedSTFT`

— Reconstructed short-time Fourier transform at the last iteration.

## More About

### Short-Time Fourier Transform

The short-time Fourier transform (STFT) is used to analyze how the frequency
content of a nonstationary signal changes over time. The magnitude squared of the STFT is
known as the *spectrogram* time-frequency representation of the signal.
For more information about the spectrogram and how to compute it using Signal Processing Toolbox™ functions, see Spectrogram Computation with Signal Processing Toolbox.

The STFT of a signal is computed by sliding an *analysis window*
*g*(*n*) of length *M* over the signal and calculating the
discrete Fourier transform (DFT) of each segment of windowed data. The window hops over the
original signal at intervals of *R* samples, equivalent to *L* = *M* –
*R* samples of overlap between adjoining segments. Most window functions taper
off at the edges to avoid spectral ringing. The DFT of each windowed segment is added to a
complex-valued matrix that contains the magnitude and phase for each point in time and
frequency. The STFT matrix has

$$k=\lfloor \frac{{N}_{x}-L}{M-L}\rfloor $$

columns, where *N _{x}* is the length
of the signal

*x*(

*n*) and the ⌊⌋ symbols denote the floor function. The number of rows in the matrix equals

*N*

_{DFT}, the number of DFT points, for centered and two-sided transforms and an odd number close to

*N*

_{DFT}/2 for one-sided transforms of real-valued signals.

The *m*th column of the STFT matrix $$X(f)=\left[\begin{array}{ccccc}{X}_{1}(f)& {X}_{2}(f)& {X}_{3}(f)& \cdots & {X}_{k}(f)\end{array}\right]$$ contains the DFT of the windowed data centered about time *mR*:

$${X}_{m}(f)={\displaystyle \sum _{n=-\infty}^{\infty}x(n)\text{\hspace{0.17em}}g(n-mR)\text{\hspace{0.17em}}{e}^{-j2\pi fn}}.$$

### Normalized Inconsistency

*Normalized inconsistency* measures the
improvement toward convergence of the reconstruction process in successive optimization
iterations.

The definition of normalized inconsistency depends on the signal reconstruction
algorithm `Method`

.

When

`Method`

is`"gla"`

,`"fgla"`

, or`"legla"`

, the inconsistency is defined as$$\text{Inconsistency}=\frac{\Vert \mathrm{STFT}\left(\mathrm{ISTFT}\left({s}_{\text{est}}\right)\right)-{s}_{\text{est}}\Vert}{\Vert {s}_{\text{est}}\Vert},$$

where

*s*_{est}is the complex short-time Fourier transform estimated at each iteration, the brackets denote the matrix norm, STFT denotes the short-time Fourier transform, and ISTFT denotes its inverse.When

`Method`

is`"gd"`

, the inconsistency is defined as$$\text{Inconsistency}=\frac{\Vert {s}_{\text{est},\text{\hspace{0.17em}}i}-{s}_{\text{est},\text{\hspace{0.17em}}i-1}\Vert}{\Vert s\Vert},$$

where

*s*_{est,i}is the estimated complex short-time Fourier transform estimated at the*i*th iteration, and*s*is the input short-time Fourier transform magnitude.

`stftmag2sig`

uses the MATLAB^{®} function `norm`

to compute matrix norms. For more
information about the STFT and its inverse, see Short-Time Fourier Transform and Inverse Short-Time Fourier Transform.

## Tips

If you are using the gradient descent method and the reconstruction is not satisfactory, set

`Display`

to`true`

. Observe the inconsistency during iterations. If the inconsistency does not decrease, reduce`StepSize`

for better reconstruction.If you are using the gradient descent method, the L-BFGS optimizer usually provides the best results. This optimizer automatically selects the step size for each iteration. However, the L-BFGS optimizer may require more computation time than other optimizers to run the same number of iterations.

## References

[1] Griffin, Daniel W., and Jae S.
Lim. "Signal Estimation from Modified Short-Time Fourier Transform." *IEEE
Transactions on Acoustics, Speech, and Signal Processing*. Vol. 32, Number 2,
April 1984, pp. 236–243. https://doi.org/10.1109/TASSP.1984.1164317.

[2] Perraudin, Nathanaël, Peter
Balazs, and Peter L. Søndergaard. "A Fast Griffin-Lim Algorithm." In *2013 IEEE
Workshop on Applications of Signal Processing to Audio and Acoustics*, New Paltz,
NY, October 20–23, 2013. https://doi.org/10.1109/WASPAA.2013.6701851.

[3] Le Roux, Jonathan, Hirokazu
Kameoka, Nobutaka Ono, and Shigeki Sagayama. "Fast Signal Reconstruction from Magnitude STFT
Spectrogram Based on Spectrogram Consistency." In *Proceedings of the 13th
International Conference on Digital Audio Effects (DAFx-10)*, Graz, Austria,
September 6–10, 2010.

[4] Ji, Li, and Zhou Tie. “On Gradient
Descent Algorithm for Generalized Phase Retrieval Problem.” In *2016 IEEE 13th
International Conference on Signal Processing (ICSP)*, 320–25. Chengdu, China:
IEEE, 2016. https://doi.org/10.1109/ICSP.2016.7877848.

## Extended Capabilities

### C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

`"gd"`

method is not supported.

### GPU Arrays

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

The `stftmag2sig`

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

`"legla"`

method is not supported.

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

## Version History

**Introduced in R2020b**

### R2023b: Reconstruct signal using gradient descent algorithm

The `stftmag2sig`

function can reconstruct a signal using the
gradient descent algorithm. You must have a Deep Learning Toolbox license to use this functionality.

## MATLAB Command

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)