Adhoc method to invert the reassigned time-frequency representation
Shristi RAJBAMSHI 1 ; Peter BALAZS 2 ; Nicki HOLIGHAUS 3
Acoustics Research Institute, Austria

Abstract
In this contribution, we present a heuristic method to invert the reassigned short time Fourier transform magnitude spectrum to allow the reconstruction of the original time domain signal. This is a simple method justinvolving an additional smearing step before phase retrieval. Finally, we provide some numerical evidence
that our method combined with existing phase retrieval methods shows excellent performance in synthesizing
pure tones and synthetic signals from the reassigned coefficients. Furthermore, we also illustrate the practicalimportance of our results by applying our method to a real-world audio signal.
Keywords: Time-frequency representation, time-frequency reassignment, invertibility

1

INTRODUCTION

Time-Frequency (TF) representation techniques such as the Short Time Fourier Transform (STFT) (1, 3), particularly its magnitude, have been long used to analyze non-stationary signals. Despite being an efficient tool for
analysis of complex signals, STFT often suffer from temporal-spectral resolution trade-off problem (2), thereby
resulting in a smeared visual representation. In order to tackle this problem, the concept of TF reassignment
was originated. This method, first introduced by Kodera (4) and later popularized by Auger and Flandrin (5),takes care of the aforementioned problem by utilizing the usually discarded phase information to produce a
sharper TF representation. In spite of this clear advantage, TF reassignment has not received much attention
because of its lack of invertibility, i.e., reconstruction of the original signal is not possible as the phase information is not available for synthesis. So far, there have been very few studies investigating the ways to reconstruct
the original signal from the reassigned coefficients. One such method based on additive sound modeling can beread in (6) and (7). In this paper, we present an approach to invert the reassigned STFT magnitude spectrum byapproximating the original STFT magnitude spectrum in combination with an existing phase retrieval method toallow the synthesis of the original signal.
1.1 Notation
In the following, we will denote matrices with uppercase letters, such as A, and column vectors as lowercase
letters, such as a. a[i, j] will denote an element of the matrix A where i and j denotes the row and column
index respectively. We will write variables as lowercase letters, e.g., s and constants as uppercase letters, e.g.,
S.
1
For a matrix A, the frobenius norm is defined as kAk f ro = (∑i ∑ j |a[i, j]|2 ) 2 . For a vector x = [x0 , x1 , · · · , xL−1 ]T ,the `1 norm is defined as kxk1 = ∑L−1
l=0 xl .
1 shristi.rajbamshi@oeaw.ac.at
2 peter.balazs@oeaw.ac.at
3 nicki.holighaus@oeaw.ac.at

2789

2

OVERVIEW OF TIME-FREQUENCY REASSIGNMENT

As per (5), TF reassignment technique can be applied to various TF representations. However, for simplicity
and clarity, we will only restrict to the STFT magnitude in this paper.
We will start by recalling a sub-sampled version of the discrete STFT called the Discrete Gabor Transform(DGT) (9). For any finite and real signal x ∈ RL , the DGT with respect to the analysis window w ∈ RL , time
step a ∈ N and number of modulations M ∈ N is given by,
L−1

Xw [n, m] = ∑ x[l + na]w[l]e−i2πml/M ,

(1)

l=0

where n ∈ {0, . . . , N−1} denotes the time-shift index, m ∈ {0, . . . , M−1} denotes the modulation index and [l+na]
is assumed to be evaluated modulo L according to circular indexing. The transform redundancy is given by theratio M/a.
If a = 1 and M = L, then (1) represents the full STFT, aka sliding FFT. In polar coordinate form, (1) canbe decomposed as,
Xw [n, m] = |Xw [n, m]|eiΦ[n,m] ,(2)where |Xw [n, m]| represents the magnitude and Φ[n, m] represents the phase.
Furthermore, the inverse DGT (IDGT) with respect to the synthesis window wd ∈ R is given as,
N−1 M−1x̄[l] = ∑ ∑ Xw [n, m]wd [l − na]ei2πm(l−na)/M ,(3)

n=0 m=0

where l ∈ {0, . . . , L−1}. Note that, if wd is a canonical dual window for w (9), then x̄ is equivalent to x up tosome numerical precision error.
Independent of the choice of window in (1), the STFT always suffers from uncertainty (2) manifesting ina visual representation either smeared in time, in frequency or in both directions, that cannot be arbitrarily reduced. In order to refine this blurred representation, reassignment technique is utilized. In general, reassignment
maps the TF data to the "true region of TF support" of the analyzed signal by means of the partial derivatives ofthe phase of the underlying TF representation. These partial derivatives are called instantaneous time-frequency
estimates. In our case, reassignment maps the magnitude at position (n, m) to a different TF position (ñ, m̃)based on the phase derivatives of the underlying STFT. These new TF positions are referred to as reassignment
operators and are computed as following (5, 10),

X 1 tw [n, m]TSR,ñ[n, m] = n + Re 
(4)
Xw [n, m]
XTSR dw [n, m]
L
Im
,
(5)
2π
Xw [n, m]where, Xw , X 1 tw and XTSR dw are the DGTs of the signal x corresponding to the analysis window w, timeTSRweighted version of the window w, i.e., T1SR tw [l] = T1SR t ·w(t)|t=lT and derivative of the windows with respectm̃[n, m] = m −SRto time, i.e., (TSR dw) [l] = TSR dw(t)dt |t=lTSR . Note that, TSR is the sampling period and we consider samples of acontinuous window function w(t), i.e. w[l] = w(t)|t=lT .
SR

Hence, the reassigned TFR for any point (n0 , m0 ) is then calculated as (5),M−1 −N−1

Rxw [n0 , m0 ] = ∑mwhere Sw [n, m] = |Xw [n, m]|.

2790

∑ Sw [n, m]δ (n0 − ñ[n, m])δ (m0 − m̃[n, m]).
n

(6)

3

PROPOSED METHOD TO OBTAIN THE ORIGINAL TIME-FREQUENCY REPRESENTATION

Although TF reassignment has many desirable properties such as TF shift invariance, energy conservation, perfect localization for pure sinusoids, chirps and impulses (5), etc, the major drawback of this method is that it is
not invertible which constrains its use. In order to address this problem, we propose a method to approximate
the original TF representation from the reassigned TF representation by implementing a two step post processing technique: The first step involves convolution of the reassigned coefficients with a kernel and the second
step involves recovery of the original STFT coefficients by employing existing phaseless retrieval methods.

x[l]|Xw [n, m]|

Xw [n, m]

Reassignment

|·|

DGT

Rxw [n0 , m0 ]Inversex[l]DGT

XPw [n, m]

Phase

CxW [n, m]

RetrievalConvolution

Figure 1. Block diagram of the proposed inverse reassignment technique
3.1 Convolution as Inverse Reassignment
In order to recover the original STFT magnitude coefficients from the reassigned coefficients Rxw [n0 , m0 ], we
perform two-dimensional convolution of the reassigned coefficients with a kernel. We will refer this procedure
as inverse reassignment. Mathematically, it can be formulated as,
CxW [n, m] = (Rxw ∗ W)[n, m] = ∑ ∑ Rxw [n0 , m0 ]W[n − n0 , m − m0 ]

(7)

m0 n0

where W[n, m] is the convolution kernel. In our work, we choose W[n, m] = w[n] ⊗ ŵ[m], where w is the same
window with respect to which the STFT was computed and ŵ is its discrete Fourier transform given by,
2π
1 K−1ŵ[m] = √ ∑ w[k]e−i K mk
K k=0

(8)

where m ∈ {0, . . . , M−1}. Further, we make an assumption that both w and ŵ are symmetric, localized around0 and kwk1 = kŵk1 = 1.
Proof in case of complex exponential
We now present a proof for the feasibility of our approach for the particular case of complex exponential
incorporating the concepts discussed in the previous sections. For that, we define a discrete complex exponential
2πsignal x ∈ RL of frequency f1 , i.e., x[l] = Aei L f1 l . Using (1), and choosing a = 1 and M = L, the STFT of x is,
2π

Xw [n, m] = Aei L f1 n ŵ[m − f1 ],(9)

Sw [n, m] =| Xw [n, m] |= Aŵ[m − f1 ],(10)

the magnitude is,2791

and the reassigned coefficients are,Rxw [n0 , m0 ] =A if m0 = f1.
0 otherwise

With the convolution kernel, W[n, m] = w[n] ⊗ ŵ[m], inverse reassignment yields,
Cxw [n, m] = (Rxw ∗ W)[n, m]= ∑ ∑ Rxw [n0 , m0 ]W[n − n0 , m − m0 )m0 n0

= ∑ ŵ[m − m0 ] ∑ Rxw [n0 , m0 ]w[n − n0 ]

(11)

n0

m0

Now, if m0 = f1 ,Cxw [n, m] = ŵ[m − f1 ] ∑ Rxw [n0 , f1 ]w[n − n0 ]n0= ŵ[m − f1 ] ∑ Aw[n − n0 ]n0= Aŵ[m − f1 ]∑ w[n − n0 ]∵ kwk1 = 1

n0

|

{z

=1

}

= Aŵ[m − f1 ].

(12)

Cxw [n, m] = |Xw (n, m)| .

(13)

Comparing (10) and (12), we obtain,

Hence, we can easily prove for the particular case of complex exponential that the operation of inverse
reassignment yields the original magnitude.
3.2 Phase Retrieval
The recovered STFT magnitude, i.e., Cxw [n, m], is then used to estimate the original STFT coefficients employing a Phase Retrieval (PR) method. PR can be defined as a method to recover the original signal from themagnitude only measurements. There are several studies focused on phase retrieval. In this section, we givea summary of three specific phase retrieval methods that we use in our work to synthesize the signal fromthe reassigned coefficients. We will denote them as PR-I, PR-II and PR-III. The first PR method is basedon an algorithm called Phase Gradient Heap Integration (PGHI) (12), second is based on an algorithm called
Griffin-Lim (GL) (11, 16) and third is based on the algorthim called Fast Griffin-Lim (13).
In PGHI, phase is retrieved by using the gradient theorem based on the assumption that the phase gradientand phase at one point is known. Phase gradient are obtained from the magnitude gradient. Since this a noniterative method, the algorithm is very fast and appropriate for long signals.
The GL algorithm is an iterative approach which converges to the original phase by repeatedly computing the
STFTs and inverse STFTs. First it performs an inverse STFT using the original magnitude and then computes
the STFT of the obtained times series. Next, it combines the original magnitude with the phase of the newly
computed STFT and then again the inverse STFT is calculated. These steps are repeated till the algorithm
converges. This algorithm is not very efficient as it converges very slowly.
The FGL algorithm is a phase retrieval method that treats the phase recovery problem as an optimization
problem. As the name suggests, this algorithm is a faster version of the traditional GL algorithm.

2792

4

SIMULATION RESULT

To evaluate the performance of our proposed method, we carry out experiments using four different test signals.
All the simulations are implemented in MATLAB using the LTFAT toolbox (14, 15). Our first test signal
is a complex exponential of 50 Hz generated as x[l] = ei2π50l , second is a synthetic signal generated as x[l] =cos(4800πl +350cos(10πl))+cos(20000πl −5400cos(2π(l −1))) +cos(750πl +80cos(4π(l −0.03))), third is the(in)famous audio signal called "glockenspiel" (14) and fourth is the speech signal called "linus" (14).
Details about all the test signals, such as their lengths (L), sampling frequencies ( f s) and parameters, a andM required for computing corresponding DGTs and IDGTs are enlisted in Table 1. For all the signals, the DGT
is computed with respect the sampled and periodized Gaussian window with variance 1 and IDGT is computedwith respect the canonical dual of the window used for DGT.
Table 1L(in samples)fs(in Hz)aMComplex Exponential

400

400

1

400

Synthetic

44100

44100

200

1000

Glockenspiel

262144

44100

200

2000

Linus

41461

8000

200

1000

For each simulation, DGT, reassignment, convolution and IDGT are computed using ‘dgt’,‘gabreassign’,‘pconv’ and ‘idgt’ fucntions available in the LTFAT toolbox. Additionally, the phase retrieval is carried out
using ‘constructPhase’ and ‘frsynabs’ routines which are also available in the LTFAT toolbox. Note that theroutines ‘constructPhase’ correspond to PR-I and ‘frsynabs’ with parameters griflim and fgriflim correspond to
PR-II and PR-III respectively.
4.1 Evaluation Metric
For each signal, we will reconstruct the original signal by employing various phase retrieval methods discussed
in Section 3.2 and calculate the consequent spectral convergence by including and excluding the convolution
step discussed in Section 3.1. We will use original phase as the benchmark to compare the performance ofvarious PR methods. The spectral convergence is given as (16),kX − X̄k f ro
20log10.
kXk f rowhere X̄ is the reconstructed STFT computed from the re-synthesized signal x̄ using 1.
4.2 Discussion
In Figure 2(a), it can be seen that our proposed method works perfectly for a pure tone as the difference
between the original and the reconstructed STFT magnitude is zero (up to numerical precision). Also for thesynthetic signal (see 2(b)), the reconstruction error is very small. Similar performance can be observed formore complex signals like glockenspiel and linus (see 2(c) and 2(d)). It should be noted that our method
performs well not only for the case of high redundancy, for example, in case of complex exponential when the
redundancy is 400, but also for the low redundancy case, for example in case of other three signals when the
redundancy is 5, 10 and 5 respectively. Remember that, redundancy is calculated as M/a.

2793

-40

-60

100
-80

50

150

-60

100

50

-100

0
0.2

0.4

0.6

0.8

-60

100

0

0.2

0.4

0.6

0.8

-60

100

0

0.2

0.4

0.6

-80

50

-100

0

-120

Time [s]

-40

150

-100

0

-120

Time [s]

-80

50

-100

0

-120

0

-80

150

-20

Frequency [Hz]

-20

-40

Frequency [Hz]

-20

-40

Frequency [Hz]Frequency [Hz]

150

-20

0.8

-120

0

0.2

0.4

Time [s]0.60.8Time [s](a) Complex Exponential

20

20

-70

10

5

0

-80

0

0.2

0.4

0.6

0.8

15

-70

10

-90

5

-100

0

1

-80

0

0.2

0.4

Time [s]

0.6

0.8

-60

Frequency [kHz]

15

20

-60

Frequency [kHz]Frequency [kHz]

-60

15

-70

10

-80

-60

Frequency [kHz]

20

15

10

-90

5

-90

5

-100

0

-100

0

1

0

0.2

0.4

Time [s]

0.6

0.8

-70

1

-80

-90

-100

0

0.2

0.4

Time [s]

0.6

0.8

1

Time [s](b) Synthetic Signal

20

20

-100

10
-120

5

15
-100

10
-120

-80

Frequency [kHz]

15

5

2

3

4

-100

10
-120

5

-80

0

1

2

3

4

-100

10
-120

-140

-140

0

Time [s]

15

5

-140

0
1

15

5

-140

0
0

20

-80

Frequency [kHz]Frequency [kHz]-80

Frequency [kHz]

20

5

0
0

1

Time [s]

2

3

4

5

0

1

Time [s]

2

3

4

5

Time [s](c) Glockenspiel
-60

-60

-60

-70

2

-80

1.5
1

-90

0.5

3

-70

2.5
2

-80

1.5
1

-90

0.5

0

-100

0

1

2

3

Time [s]

4

5

3.5

3

-70

2.5
2

-80

1.5
1

-90

0.5

0

-100

0

1

2

3

4

5

Frequency [kHz]

3
2.5

-60

3.5

Frequency [kHz]3.5

Frequency [kHz]Frequency [kHz]

3.5

3

-70

2.5
2

-80

1.5
1

-90

0.5

0

-100

0

1

Time [s]

2

3

4

5

0

-100

0

1

2

3

Time [s]Time [s]

(III)

(IV)

4

5

(d) Linus

(I)

(II)

Figure 2. (I) left: Original STFT magnitude, (II) left middle: Reassigned STFT magnitude, (III) right middle:
Reconstructed STFT magnitude (Inverse Reassignment), (IV) right: difference between (I) and (III)The Table 2 depicts the spectral convergence, using different phase retrieval methods while including andexcluding the convolution step denoted as "C" and "NC" respectively.

2794

Table 2. Spectral Convergence
Complex Exponential

Synthetic

Glockenspiel

Linus

NC

C

NC

C

NC

C

NC

C

Original Phase-7.6 dB-169.63 dB-8.67 dB

-28.1 dB

-8.79 dB

-31.06 dB-9.27 dB

-25.7 dB

PR-I

-7.65 dB

-169.63 dB-7.89 dB

-13.02 dB-6.64 dB

-14.95 dB-8.46 dB

-10.3 dB

PR-II

-7.65 dB

-21.19 dB-8.04 dB

-16.62 dB-6.72 dB

-14.39 dB-8.69 dB

-16.55 dB

PR-III-7.65 dB

-23.86 dB-8.10 dB

-19.20 dB-6.78 dB

-14.73 dB-8.70 dB

-17.46 dB

It can be deduced from Table 2 there is significant improvement in spectral convergence if we convolve thereassigned coefficient with the kernel before we perform the phase retrieval.

5

CONCLUSION

We have demonstrated that reassigned STFT magnitude can be inverted using the convolution technique basedon the assumption that the analysis window and its Fourier transform are symmetric and their `1 -norm is 1. The
proposed inverse reassignment method in combination with the existing PR techniques can be used to recover
the original time domain signal. Our method exhibited excellent performance not only for pure tones but also
for complex signals like audio and speech. Furthermore, pilot tests indicated that there was no audible difference between the original and reconstructed audio and speech signals. Therefore, we believe that this method
can be very useful in extending the usage of reassignment technique for present and future signal processing
applications which has been limited thus far due to the lack of invertibility.

ACKNOWLEDGEMENT
This work was supported by the Austrian Science Fund (FWF) START-project FLAME (“Frames and Linear
Operators for Acoustical Modeling and Parameter Estimation”; Y 551-N13).

REFERENCES
1. Gröchenig K. Foundations of Time-Frequency Analysis, Applied Numerical Harmonic Analysis, Birkhäuser,
Boston, MA, 2001.
2. Gröcheni, K. Uncertainty Principles for Time-Frequency Representations, Chapter 2 in Advances in Gabor
Analysis (eds. H. G. Feichtinger and T. Strohmer), Birkhäuser Boston, 2003, pp 11-30.
3. Cohen, L. Time-Frequency Analysis: Theory and Applications, Prentice Hall Signal Processing Series, Prentice Hall, 1995
4. Kodera, K.; Gendrin R.; Villedary C. Analysis of time-varying signals with small bt values, IEEE Transactions on Acoustics, Speech and Signal Processing, Vol 26, Feb 1978, pp 64–76.
5. Auger F.; Flandrin P. Improving the readability of time-frequency and time-scale representation by the reassignment method, IEEE Transactions on Signal Processing, Vol 43(5), May 1995, pp 1068–1088.
6. Fitz, K.; Haken, L. On the use of time frequency reassignment in additive sound modeling, Journal of the
Audio Engineering Society, Vol 50(11), Nov 2002, pp 879-93.

2795

7. Han, L., Sacchi, MD. Spectral decomposition and de-noising via time-frequency and space-wavenumber reassignment, Geophysical Prospecting, Vol 62(2), Mar 2014, pp 244-57.
8. Søndergaard, P.L. Gabor frames by Sampling and Periodization, Advanced in Computational Mathematics,
Vol. 27(4), Nov 2007, pp 355 –373.
9. Søndergaard, P.L.; Hansen P.C.; Christensen, O. Finite discrete Gabor analysis, Ph.D. dissertation, Technical
University of Denmark, 2007.
10. Fenet, S.; Badeau, R.; Richard, G. Reassigned time–frequency representations of discrete time signals andapplication to the Constant-Q Transform, Signal Processing, Elsevier,Vol 132, Mar 2017, pp 170-6.
11. Griffin, D.; Lim, J. Signal estimation from modified short-time Fourier transform, IEEE Transactions on
Acoustics, Speech and Signal Processing, Vol 32(2), 1984, pp 236-243.
12. Prusa, Z.; Balazs, P.; Søndergaard P.L. A non-iterative method for STFT phase (re)Construction, IEEE/ACM
Transactions on Audio, Speech, and Language Processing, Vol 25(5), May 2017, pp 1154-64.
13. Perraudin, N.; Balazs, P.; Søndergaard, P.L. A fast Griffin-Lim algorithm, IEEE Workshop on Applicationsof Signal Processing to Audio and Acoustics, October 2013, pp 1-4.
14. Søndergaard, P.L.; Torresani, B.; Balazs, P. The linear time frequency analysis toolbox, International Journalof Wavelets, Multiresolution and Information Processing, Vol 10(04), 2012, p. 1250032.
15. Prusa, Z.; Søndergaard, P.L., Holighaus, N.; Wiesmeyr, C.; Balazs, P. The Large Time-Frequency Analysis
Toolbox 2.0 in Sound, Music, and Motion, ser. LNCS. Springer International Publishing, 2014, pp 419–442.
16. Sturmel, N.; Daudet, L. Signal reconstruction from STFT magnitude: A state of the art. International conference on digital audio effects (DAFx), Sep 2011, pp 375-386.

2796