Estimates of the Reconstruction Error in Partially Redressed Warped Frames Expansions

TL;DR

This paper investigates the approximation errors involved in using compactly supported analysis-synthesis elements in warped frames, which are useful for non-uniform time-frequency audio analysis, especially in online processing scenarios.

Contribution

It provides a detailed analysis of the reconstruction error when using finite support warped frames, including examples and case studies.

Findings

Approximation errors are quantifiable and depend on support size.

Finite support introduces small but manageable reconstruction errors.

The study offers practical insights for real-time audio processing applications.

Abstract

In recent work, redressed warped frames have been introduced for the analysis and synthesis of audio signals with non-uniform frequency and time resolutions. In these frames, the allocation of frequency bands or time intervals of the elements of the representation can be uniquely described by means of a warping map. Inverse warping applied after time-frequency sampling provides the key to reduce or eliminate dispersion of the warped frame elements in the conjugate variable, making it possible, e.g., to construct frequency warped frames with synchronous time alignment through frequency. The redressing procedure is however exact only when the analysis and synthesis windows have compact support in the domain where warping is applied. This implies that frequency warped frames cannot have compact support in the time domain. This property is undesirable when online computation is required. Approximations in which the time support is finite are however possible, which lead to small reconstruction errors. In this paper we study the approximation error for compactly supported frequency warped analysis-synthesis elements, providing a few examples and case studies.

Figures

Equations

$$\begin{array}{c}\left[\mathcal{S}s\right](\tau ,\nu )=⟨s,h_{\tau ,\nu }⟩=⟨s,{\mathbf{T}}_{\tau }{\mathbf{M}}_{\nu }h⟩=\\ {\int }_{-\infty }^{+\infty }s(t)\overline{h_{0,0}(t-\tau )}{e}^{-j2\pi \nu (t-\tau )}dt,\end{array}$$

$$\mathcal{G}(h,a,b)=\{{\mathbf{T}}_{na}{\mathbf{M}}_{qb}h:n,q\in \mathbb{Z}\}$$

$$⟨s,{\mathbf{T}}_{na}{\mathbf{M}}_{qb}h⟩,n,q\in \mathbb{Z}$$

$$A\Vert s{\Vert }^2\le \sum _{l\in I}|⟨s,{\psi }_l⟩{|}^2\le B\Vert s{\Vert }^2\forall s\in \mathcal{H},\text{vspace}-2mm$$

$${\mathbf{W}}_{\tilde{\theta }}={\mathcal{F}}^{-1}{\mathbf{W}}_{\theta }\mathcal{F},$$

$${\hat{s}}_{fw}(\nu )=\left[\widehat{{\mathbf{U}}_{\tilde{\theta }}s}\right](\nu )=\sqrt{\left|\frac{d\theta }{d\nu }\right|}\hat{s}(\theta (\nu )),$$

$$\begin{array}{c}\left[{\mathcal{S}}_{\tilde{\theta }}s\right](\tau ,\nu )=\left[\mathcal{S}{\mathbf{W}}_{\tilde{\theta }}s\right](\tau ,\nu )=\\ ⟨{\mathbf{W}}_{\tilde{\theta }}s,h_{\tau ,\nu }⟩=⟨s,{\mathbf{W}}_{\tilde{\theta }}^{†}{h}_{\tau ,\nu }⟩,\end{array}$$

$$\begin{array}{c}{\hat{\tilde{h}}}_{\tau ,\nu }(f)=\left[\hat{{\mathbf{W}}_{{\tilde{\theta }}^{-1}}{h}_{\tau ,\nu }}\right](f)=\\ \sqrt{\frac{d{\theta }^{-1}}{df}}\hat{h}({\theta }^{-1}(f)-\nu )e^{-j2\pi {\theta }^{-1}(f)\tau },\end{array}$$

$$\begin{array}{c}{\phi }_{n,q}={\mathbf{T}}_{na}{\mathbf{M}}_{qb}h\text{hfill}\\ {\gamma }_{n,q}={\mathbf{T}}_{na}{\mathbf{M}}_{qb}g,\text{hfill}\end{array}$$

$$\begin{array}{c}{\tilde{\phi }}_{n,q}={\mathbf{U}}_{\tilde{\theta }}^{†}{\phi }_{n,q}={\mathbf{U}}_{{\tilde{\theta }}^{-1}}{\mathbf{T}}_{na}{\mathbf{M}}_{qb}h\text{hfill}\\ {\tilde{\gamma }}_{n,q}={\mathbf{U}}_{\tilde{\theta }}^{†}{\gamma }_{n,q}={\mathbf{U}}_{{\tilde{\theta }}^{-1}}{\mathbf{T}}_{na}{\mathbf{M}}_{qb}g,\text{hfill}\end{array}$$

$$s=\sum _{n,q\in Z}⟨s,{\tilde{\phi }}_{n,q}⟩{\tilde{\gamma }}_{n,q}.$$

$$\left[\hat{{\mathbf{W}}_{{\tilde{\theta }}^{-1}}\mathcal{S}{\mathbf{W}}_{\tilde{\theta }}s}\right](f,\nu )=\overline{{\hat{h}}_{0,0}({\theta }^{-1}(f)-\nu )}\hat{s}(f),$$

$${\widehat{\tilde{\stackrel{~}{h}}}}_{\tau ,\nu }(f)=\left[\widehat{{\mathbf{T}}_{\tau }{\mathbf{W}}_{\tilde{\theta }}{h}_{0,\nu }}\right](f)={\hat{h}}_{0,0}({\theta }^{-1}(f)-\nu )e^{-j2\pi f\tau },$$

$${\mu }_m(n)={\int }_{-\frac{1}{2}}^{+\frac{1}{2}}\sqrt{\frac{d\vartheta }{d\nu }}{e}^{j2\pi (n\vartheta (\nu )-m\nu )}d\nu ,$$

$$\tilde{x}(m)=\left[{\mathbf{D}}_{\tilde{\vartheta }}x\right](m)={⟨x,{\mu }_m⟩}_{{\ell }^2(\mathbb{Z})}.$$

$$\hat{\tilde{x}}(\nu )=\sqrt{\frac{d\vartheta }{d\nu }}\hat{x}(\vartheta (\nu )),$$

$$\begin{array}{c}{\widetilde{\stackrel{~}{\phi }}}_{n,q}={\mathbf{D}}_{{{\tilde{\vartheta }}_q}^{-1}}{\tilde{\phi }}_{\bullet ,q}=\sum _m{\eta }_{n,q}(m){\tilde{\phi }}_{m,q}\text{hfill}\\ {\widetilde{\stackrel{~}{\gamma }}}_{n,q}={\mathbf{D}}_{{{\tilde{\vartheta }}_q}^{-1}}{\tilde{\gamma }}_{\bullet ,q}=\sum _m{\eta }_{n,q}(m){\tilde{\gamma }}_{m,q},\text{hfill}\end{array}$$

$$s=\sum _{n,q\in (Z)}⟨s,{\widetilde{\stackrel{~}{\phi }}}_{n,q}⟩{\widetilde{\stackrel{~}{\gamma }}}_{n,q}.$$

$${\widehat{\tilde{\stackrel{~}{\phi }}}}_{n,q}(f)=A(f)\hat{h}({\theta }^{-1}(f)-qb)e^{-j2\pi n{\vartheta }_q(a{\theta }^{-1}(f))},$$

$$A(f)=\sqrt{\frac{d{\theta }^{-1}}{df}}{\sqrt{\frac{d{\vartheta }_q}{d\nu }}|}_{\nu =a{\theta }^{-1}(f)}.$$

$${\vartheta }_q(a{\theta }^{-1}(f))=d_{q}f$$

$${\widehat{\tilde{\stackrel{~}{\phi }}}}_{n,q}(f)=\sqrt{\frac{d_q}{a}}\hat{h}({\theta }^{-1}(f)-qb)e^{-j2\pi n{d}_{q}f}.$$

$$abK\le 1$$

$$d_q{\text{B}}_q\le 1$$

$${\text{B}}_q=\theta \left(qb+\frac{\text{B}}{2}\right)-\theta \left(qb-\frac{\text{B}}{2}\right)$$

$$a=\frac{T_h}{R}$$

$$b$$

$$q_{\sup }=\frac{{\theta }^{-1}(SR/2)}{b}.$$

$${\text{T}}_c=\frac{{\text{T}}_h}{{\theta }_{\inf }^{\prime }}{C}_{{\text{T}}_c}$$

$${\text{B}}_q=\theta \left(qb+\frac{Kb}{2}\right)-\theta \left(qb-\frac{Kb}{2}\right)\simeq {\theta }^{\prime }(qb)Kb$$

$$\Rightarrow d_q=\frac{1}{{\text{B}}_q}\simeq \frac{1}{{\theta }^{\prime }(qb)Kb}$$

Full text

1 Introduction

The availability of configurable time-frequency schemes is an asset for sound analysis and synthesis, where the elements of the representation can efficiently capture features of the signal. These include features of interpretation: e.g. glissandi and vibrati; human perception: e.g. non-uniform frequency sensitivity of the cochlea; music theory: e.g. scales (pentatonic, 12-tone, Indian scales with unequally spaced tones, etc.); physical effects: e.g. non-harmonic overtones in the low register of the piano or of percussion instruments and so forth. When used in the context of Music Information Retrieval, the adaptation of the representation to music scales is bound to improve. e.g. for instrument-, note-, chord- and overtone detection and recognition.

Traditionally, two extreme cases have been considered: Gabor expansions featuring uniform time and frequency resolutions and orthogonal wavelet expansions and frames featuring octave band allocation and constant uncertainty of the representation elements. In previous work [1, 2, 3, 4, 5, 6, 7], generalised Gabor frames have been constructed which allow for non-uniform time-frequency schemes with perfect reconstruction. In [3] the allocation of generalised Gabor atoms is specified according to a frequency or time warping map. In [8] the STFT redressing method is introduced, which, with the use of additional warping in time-frequency, shows under which conditions one can have generalised Gabor frames. These conditions stem from the interaction of sampling in time-frequency and frequency or time warping operators, which allow to incorporate the results in [3] in a more general context. It is shown that arbitrary allocation of the atoms is exactly possible in the so called painless case, i.e. in the case of finite time support of the windows for arbitrary time interval allocation and of finite frequency support of the windows for arbitrary frequency band allocation. Non-uniform frequency analysis by means of warping was introduced in [7]. Non-uniform Gabor frames with constant-Q were previously introduced in [1] , based on the theory developed in [2], where an ad hoc procedure was employed for their construction. In [3] we provided an alternative general method for their construction, using warping. In [8] the redressing method was introduced and applied to the general construction of non-stationary Gabor frames. In [9] the general theory was revisited with the real-time computational aspects in mind. There, the first approximate schemes were introduced without extensive testing.

Since online computation of the generalised Gabor analysis-synthesis is only possible with finite duration windows, the arbitrary frequency band allocation does not lead to an exact solution in applications that require real-time, while the arbitrary time interval allocation presents little or no problem. In [9] approximations leading to nearly exact representations were introduced.

In this paper we expand on the results in [9] and provide a study of the approximation error on a wide class of signals when finite duration windows are required in arbitrary frequency band allocation.

The paper is organised as follows. In Section 2 we review the concept of applying time and frequency warping to time-frequency representations, together with the redressing method, which involves a further warping operation in the time-frequency domain to reduce or eliminate dispersion. In Section 3 we introduce approximations suitable for the online computation of redressed frame expansions. In Section 4 the results of numerical experiments are shown, which provide estimations of the approximation error. In Section 5 we draw our conclusions.

2 Redressed Warped Gabor Frames

In this section we review the concepts leading to the definition of redressing in the context of frequency warped time-frequency representations. First we review the basic notions of STFT (Short-Time Fourier Transform) and Gabor frames. Then we move on to the definition of warped frames and then to the redressing procedure.

2.1 STFT and Gabor Frames

Gabor expansions can be considered as a form of sampling and exact reconstruction of the STFT. As is well-known, given a window $h$ and defining the time-shift operator ${\mathbf{T}}_{\tau}s(t)=s(t-\tau)$ and the modulation operator ${\mathbf{M}}_{\nu}s(t)=e^{j2\pi\nu t}s(t)$, the STFT is obtained by applying the operator $\mathcal{S}$ to the signal $s$:

[TABLE]

where the overbar denotes complex conjugation.

A Gabor system is generated by the kernel of $\mathcal{S}$ by sampling the time-frequency plane $(\tau,\nu)$:

[TABLE]

where $a,b>0$ are sampling parameters.

The scalar products of the signal with the members of the Gabor system

[TABLE]

provide evaluations of the STFT (1) of a signal $s$ with window $h$ at the time-frequency grid of points $(na,qb)$, with $n,q\in\mathbb{Z}$. The question whether the signal $s$ can be reconstructed from these evaluations can be addressed by introducing the concept of frame.

A sequence of functions $\left\{\psi_{l}\right\}_{l\in I}$ in the Hilbert space $\mathcal{H}$ is called a frame if there exist both positive constant lower and upper bounds $A$ and $B$, respectively, such that

[TABLE]

where $\|s\|^{2}=\langle s,s\rangle$ is the norm square or total energy of the signal. Frames generate signal expansions, i.e., the signal can be perfectly reconstructed from its projections over the frame.

A Gabor system that is a frame is called a Gabor frame. In this case, the signal can be reconstructed from the corresponding samples of the STFT (3). While not unique, reconstruction can be achieved with the help of a dual frame, which in turn is a Gabor frame generated by a dual window ${\tilde{h}}$. Perfect reconstruction essentially depends on the choice of the window and the sampling grid. One can show that there exist no Gabor frames when $ab>1$. See [10] for more informations about Gabor frames.

2.2 Warped STFT and Gabor Frames

The warped STFT can be obtained by warping the signal prior to applying the STFT operator. In this paper we focus on pure frequency warping.

A frequency warping operator $\mathbf{W}_{\tilde{\theta}}$ is completely characterised by a function composition operator $\mathbf{W}_{\theta}$, such that ${{\mathbf{W}}_{\theta}x}=x\circ\theta$, in the frequency domain:

[TABLE]

where $\mathcal{F}$ is the Fourier transform operator. The function $\theta$ is the frequency warping map, which transforms the Fourier transform $\hat{s}=\mathcal{F}s$ of a signal $s$ into the Fourier transform $\hat{s}_{fw}=\mathcal{F}s_{fw}$ of another signal $s_{fw}$. We affix the $\tilde{}$ symbol over the map $\theta$ as a reminder that the map operates in the frequency domain.

If the warping map is one-to-one and almost everywhere differentiable then a unitary form $\mathbf{U}_{\tilde{\theta}}$ of the warping operator can be defined by the following frequency domain action

[TABLE]

where $\nu$ denotes frequency. We assume henceforth that all warping maps are almost everywhere increasing so that the magnitude sign can be dropped from the derivative under the square root.

Given a frequency warping operator $\mathbf{W}_{\tilde{\theta}}$, the warped STFT is defined through the operator $\mathcal{S}_{\tilde{\theta}}$ as follows

[TABLE]

which is indeed a warped version of (1), where $\mathbf{W}^{\dagger}_{\tilde{\theta}}$ is the adjoint of the warping operator. If the warping operator is unitary then we have $\mathbf{W}^{\dagger}_{\tilde{\theta}}=\mathbf{W}_{\tilde{\theta}}^{-1}=\mathbf{W}_{\tilde{\theta}^{-1}}$. In that case, warping the signal prior to STFT is perfectly equivalent to perform STFT analysis with inversely frequency warped windows. The warped STFT is unitarily equivalent to the STFT so that a number of properties concerning conditioning and reconstruction hold [11].

The Fourier transforms of the frequency warped STFT analysis elements are

[TABLE]

i.e., the warped STFT analysis elements are obtained from frequency warped modulated windows centred at frequencies $f=\theta(\nu)$. The windows are time-shifted with dispersive delay, where the group delay is $\tau\tfrac{{d\theta^{-1}}}{{df}}$.

Frequency warping generally disrupts the time organisation of signals due to the fact that the time-shift operator ${\mathbf{T}}_{\tau}$ does not commute with the frequency warping operator [9].

From (4) it is easy to see that any unitary operation, in particular unitary warping, on a frame results in a new frame with the same frame bounds $A$ and $B$ [11]. Since the atoms are not generated by shifting and modulating a single window function, the resulting frames are not necessarily of the Gabor type. However, warping prior to conventional Gabor analysis and unwarping after Gabor synthesis always leads to perfect reconstruction.

Starting from a Gabor frame (analysis) $\left\{{{\varphi}_{n,q}}\right\}_{q,n\in{\mathbb{Z}}}$ and dual frame (synthesis) $\left\{{{\gamma}_{n,q}}\right\}_{n,q\in{\mathbb{Z}}}$:

[TABLE]

where $h$ and $g$ are dual windows, warped frames can be generated, following (8), by unwarping the analysis and synthesis frames. In the case of non-unitary warping, a frequency domain scaling operation is necessary in order to reconstruct the original signal. For the case of unitary warping we simply have:

[TABLE]

where $\left\{{\tilde{\varphi}_{n,q}}\right\}_{q,n\in{\mathbb{Z}}}$ is the frequency warped analysis frame and $\left\{{\tilde{\gamma}_{n,q}}\right\}_{n,q\in{\mathbb{Z}}}$ is the dual warped frame for the synthesis. With these definitions, one obtains the signal expansion

[TABLE]

Warped Gabor frames suffer from the same problem as the warped STFT. Indeed the Fourier transforms of the warped Gabor frame elements bear frequency dispersive delays so that dispersive time samples are produced by the direct application of the frequency warped frame analysis.

2.3 Redressing Methods

As shown in [8, 9], the dispersive delays intrinsic to the warped STFT can be redressed, i.e. made into constant delays in each analysis band if frequency unwarping is performed in the transformed time domain, i.e. with respect to time shift. In other words, instead of (7) we consider the similarity transformation ${{\mathbf{W}_{\tilde{\theta}}^{\dagger}\mathcal{S}\mathbf{W}_{\tilde{\theta}}}}$ on the STFT operator, which is time-shift covariant. In fact, one has:

[TABLE]

which is in the form of a time-invariant filtering operation, corresponding to convolution in time domain. The filters are frequency warped versions of the modulated windows in the traditional STFT. The Fourier transform of the redressed analysis elements are

[TABLE]

which shows that the dispersive delays in the analysis elements (8) are brought back to non-dispersive delays.

In redressing warped Gabor frames one faces a further difficulty due to time-frequency sampling. In this case, inverse frequency warping can only be applied to sequences (with the respect to the time shift index) and may not perfectly reverse the dispersive effect of the original map on delays.

Unitary frequency warping in discrete time can be realised with the help of an orthonormal basis of $\ell^{2}(\mathbb{Z})$ constructed from an almost everywhere differentiable warping map $\vartheta$ that is one-to-one and onto $[-\tfrac{1}{2},+\tfrac{1}{2}[$, as follows:

[TABLE]

where $n,m\in\mathbb{Z}$ (see [12, 13, 14, 15, 16]). The map $\vartheta$ can be extended over the entire real axis as congruent modulo $1$ to a $1$-periodic function.

Given any sequence $\{x(n)\}$ in $\ell^{2}(\mathbb{Z})$, the action of the discrete-time unitary warping operator $\mathbf{D}_{{\tilde{\vartheta}}}$ is defined as follows:

[TABLE]

In fact, the sequence $\{\tilde{x}(m)\}$ in $\ell^{2}(\mathbb{Z})$ satisfies

[TABLE]

where the $\hat{}$ symbol, when applied to sequences, denotes discrete-time Fourier transform. The sequences $\overline{\eta_{m}(n)}$ define the nucleus of the inverse unitary frequency warping $\ell^{2}(\mathbb{Z})$ operator $\mathbf{D}_{{\tilde{\vartheta}}^{-1}}=\mathbf{D}_{{\tilde{\vartheta}}}^{\dagger}$. where $\eta_{m}(n)={\mu_{n}(m)}$.

In order to limit or eliminate time dispersion in the frequency warped Gabor expansion, the discrete-time frequency warping operator $\mathbf{D}_{{\tilde{\vartheta}}^{-1}}$ is applied to the time sequence of expansion coefficients over the warped Gabor frames. Since the operator is applied only on the time index, for generality, one can include dependency of the map and of the sequences $\eta_{n}$ on the frequency index $q$. The process can be equivalently described by defining the redressed frequency warped Gabor analysis and synthesis frames as follows:

[TABLE]

obtaining:

[TABLE]

One can show [9] that the Fourier transforms of the redressed frame are

[TABLE]

where

[TABLE]

Hence, dispersion is completely eliminated if

[TABLE]

for any $f\in\mathbb{R}$, where $d_{q}$ are positive constants controlling the time scale in each frequency band. In this case, the Fourier transforms of the redressed frame elements become:

[TABLE]

When all $d_{q}$ are identical all the time samples are aligned to a uniform time scale throughout frequencies. If the $d_{q}$ are distinct, time realignment when displaying the non-uniform spectrogram is a simple matter of different time base or time scale for each frequency band.

Unfortunately, due to the discrete nature of the redressing warping operation, each map $\vartheta_{q}$ is constrained to be congruent modulo $1$ to a $1$-periodic function, while the global warping map $\theta$ can be arbitrarily selected. Moreover, the functions $\vartheta_{q}$ must be one-to-one in each unit interval, therefore they can have at most an increment of $1$ there.

In Fig. 1 we illustrate the phase linearisation problem. There, the black curve is the amplitude scaled warping map $d_{q}\theta(\nu)$ and the grey curve represents the map $\vartheta_{q}(a\nu)$, which is $1/a$-periodic. Both maps are plotted in the abscissa $\nu=\theta^{-1}(f)$. By amplitude scaling the warping map $\theta$ one can allow the values of the map to lie in the range of the discrete-time warping map $\vartheta_{q}$. The amplitude scaling factors happen to be the new time sampling intervals $d_{q}$ of the redressed warped Gabor expansion.

In the “painless” case, which was hand picked in [3], the window $h$ is chosen to have compact support in the frequency domain. Through equation (19) and condition (21), the redressing method shows that, given any continuous and almost anywhere differentiable and increasing warping map, only in the painless case one can exactly eliminate the dispersive delays with the help of (17). In fact, in this case linearisation of the phase is only required within a finite frequency range given by the frequency support of the frame elements [8, 9], which is compatible with the periodicity constraint of the redressing maps $\vartheta_{q}$.

In the general case, a perfect time realignment of the components is not guaranteed. Notice, however, that, by construction, the redressed warped Gabor systems are guaranteed to be frames for any choice of the maps $\vartheta_{q}$ satisfying the stated periodicity conditions, even when the phase is not completely linearised. Locally, within a certain band of the warped modulated windows it is possible to linearise the phase of the complex exponentials in (19). In the sequel we will refer to this band as the essential bandwidth since, hopefully, the magnitude of the Fourier transform of the window is negligible outside it, at least as a design goal.

Unfortunately, in general both the painless case and the partially redressed cases lead to infinite duration windows, which are undesirable when online computation is required. In the next we are going to propose some approximations and study the reconstruction error also through numerical experiments.

3 Real-time computation of the Warped Gabor Expansion

For real-time computation one needs to make assumptions on the signal, as well on the window functions $h$ in order to keep the computational load as low as possible.

By requiring the window $h$ to be real-valued and $\theta$ to be odd, we obtain $\hat{\tilde{\tilde{\varphi}}}_{n,-q}(-f)$ = $\overline{\hat{\tilde{\tilde{\varphi}}}_{n,q}(f)}$. If, additionally the input signal is real-valued, then the coefficients $c_{n,q}=\langle s,\tilde{\tilde{\varphi}}_{n,q}\rangle$ fulfil $c_{n,-q}=\overline{c_{n,q}}$ and thus we only need to compute about half the coefficients and frame elements. By enforcing shift-invariance of the warped frame elements (i.e. $\tilde{\tilde{\varphi}}_{n,q}(t)=\tilde{\tilde{\varphi}}_{0,q}(t-na)$), which we must do since otherwise the pre-computation of the frame elements and hence the real-time computation of the expansion would be impossible, it is sufficient to only store $\tilde{\tilde{\varphi}}_{0,q}$ for non-negative values of q. It is advisable that the warping map is selected as a continuously differentiable function, since the error in the resynthesised signal in the frequency domain at the points of discontinuity of $\theta^{\prime}$ is very high. A strictly positive derivative helps the atoms not get too stretched in the time domain, which is undesirable in real time applications. To avoid aliasing we further require that $\theta(qb)=\text{SR}/2$ and $\theta^{\prime\prime}(qb)=0$ for some $q\in\mathbb{N}$, where SR denotes the sampling rate. This ensures that the high frequencies are smoothly mapped back to the negative frequencies and avoids extra aliasing introduced by the approximation. By requiring the Fourier transform of the window $h$ to have an analytic expression, we can avoid numerical errors in the computation of the warped windows. We further only worked with windows which are dual to itself, i.e. $\varphi_{n,q}=\gamma_{n,q}$ which is actually not a necessary restriction.

3.1 Implementation Details

We implemented the approximate warped redressed Gabor analysis and synthesis as C externals interfaced to Pure-Data (32 bit and 64 bit). It runs both under Windows and Linux and can use multiple cores. The warped windows are computed using equation (22) and applying an IDFT of that data. The inner products $c_{n,q}$ are directly computed with a loop. Also for the synthesizing part, the sum $\sum_{n,q}c_{n,q}\tilde{\tilde{\varphi}}_{n,q}$ is directly computed. It turned out that this is sufficiently fast for real-time computation and whence we made no use of fast convolution algorithm. Since the data rate of the analysis part is not uniform in time, PD’s signal connections are not suitable to transfer the coefficients. Therefore we use the signal connections to transfer pointers rather than signals.

3.2 Interface Details

For simplicity we require that the essential bandwidth B is a multiple $K\in\mathbb{N}$ of the frequency shift parameter $b$, i.e. $\text{B}=Kb$. Then, in order to obtain a frame the following conditions must be fulfilled as can be seen easily in Figure 1 [9, eq. (34) to (36)]:

[TABLE]

where $\text{B}_{q}$ is the ess. bandwidth of the warped modulated window.

[TABLE]

In the case of an exponential increasing warping map, it makes sense to set the frequency shifting parameter in a way that adjacent notes fall away from the windows main lobe in the frequency domain (If there is such a main lobe, as in the case of the raised cosine window). This on the other hand determines the window length $\text{T}_{h}$, a minimal value for $R$ and the time shift parameter

[TABLE]

with $R\geq K$ due to equation (23). We set $R=K$ for simplicity. Equations (23) are fulfilled by setting

[TABLE]

where $C_{b},C_{d}$ can be chosen inside the PD-patch. $C_{b}$ together with $R$ controls the oversampling, where $C_{d}$ controls the bandwidth in which the phase is linearised and also influences the oversampling. The $d_{q}$ are different for each band which results in a non-uniform data rate for each band. We remark that the numbers $d_{q}$ have to be chosen, such that $d_{q}/\text{SR}\in\mathbb{N}$, where SR is the sampling rate.

The number of bands $q_{\sup}$ we need for a given SR is

[TABLE]

The assumptions on the warping map ensure that this is a natural number.

Due to the warping, the supports of the windows are in general unbounded. Thus we compute the windows with a zero padded array of length of $\text{T}_{c}$ ($c$ for compute), which is defined as

[TABLE]

where $C_{\text{T}_{c}}\geq 1$ can be chosen inside the PD-patch and $\theta^{\prime}_{\inf}=\inf_{f\in\mathbb{R}}\theta^{\prime}(f)$. Due to the same reason it is indispensable to cut the windows after warping. To define a sensible atom length $T_{q}$ after warping, we set the atoms to zero after the point from which they were smaller than $\left|\tilde{\tilde{\varphi}}_{q}(t)\right|<1/C_{\text{cut}}\max_{t}\left|\tilde{\tilde{\varphi}}_{q}(t)\right|$, with $C_{\text{cut}}$ a constant which can be chosen inside the PD-patch. Afterwards the windows are truncated accordingly and with respect to the parameter $\text{T}_{\max}$, which defines the maximal desired window length. It is important that the windows are cut after aligning them to the time origin. The second approach suggested in [9] to obtain windows with finite length by computing only a linearised version of the warped windows, leads to very bad reconstruction.

3.3 Computational Costs

3.3.1 Analysis

A rough estimation yields the following. Let $\text{T}_{q}$ denote the length of the $q^{th}$-window . It will be clear from the context if $\text{T}_{q}$ denotes seconds or samples. To compute the inner product of that window with the signal one needs $2\text{T}_{q}$ many real multiplications and summations. This has to be done every $d_{q}$ samples. Hence, per sample we have $4\text{T}_{q}/d_{q}$ many floating point operations (flops) in average. If the essential bandwidth is not too large, then $\text{T}_{q}$ is proportional to $\text{T}_{0}/\theta^{\prime}(qb)$ and $d_{q}$ is proportional to $1/(Kb\theta^{\prime}(qb))$ since, linearising $\theta$ around $qb$ yields $\theta(qb+\frac{Kb}{2})\simeq\theta(qb)+\theta^{\prime}(qb)\frac{Kb}{2}$ and hence

[TABLE]

Summing up over all windows we get the average number of operations per sample $N_{\text{avg}}$:

[TABLE]

Where $q_{\sup}$, defined above, denotes the number of bands and depends on $\theta^{-1}$. That means the complexity of the analysis part is proportional to $\text{T}_{0}$, $K$ and $b$ and $q_{\sup}\sim\theta^{-1}(SR)$.

If there is a lower bound for the $d_{q}$’s, one can choose identical numbers for all bands. However, this increases the computational load tremendously, making real time computation impossible.

The above is an estimation of the average computational cost. In the worst case all inner products of the windows with the signal have to be computed starting at one frame. The number of operations for that frame is

[TABLE]

However, if one processes the audio stream block wise, the worst case cannot arise for all samples in the block at once, because two worst case scenarios have a distance of at least $M=\max_{q}{d_{q}}$ (actually $M=\text{lcm}(\{d_{q}\})$), the next $M$ samples have for sure lower computational cost. Furthermore, the next $m=\min_{q}{d_{q}}$ samples have zero computational cost. Therefore the average costs are a suitable measure for the analysis process.

3.3.2 Synthesis Part

The complexity of the synthesis part is the same as that of the analysis part. Furthermore, in the synthesis part the worst case scenario can be avoided, because only the parts of the next frame buffer have to be computed in real-time, the rest can be computed later. Nevertheless, our given implementation is not optimised in this direction.

3.3.3 Memory Costs

Our algorithm for precomputing the windows needs $2q\text{T}_{0}$ cells. With slight modifications it only needs $\text{T}_{0}+\sum_{q}\text{T}_{q}$. The smallest possible number of cells being needed is $\sum_{q}\text{T}_{q}$. The analysis and synthesis algorithm both need at least a buffer of size $\textrm{\emph{audiobuffer}}+\text{T}_{\max}$, where $\text{T}_{\max}$ denotes the window length of the longest windows used in the analysis and synthesis, and audiobuffer the buffer length in which the audio is processed.

The frame elements for our tests below needed between 50 and 400 MB.

4 Computational Error

The measured $err$-numbers are the difference between the averaged RMS-amplitude in dB of the input signal and the analysed-synthesised output signal (i.e. negative signal to noise ratio). For comparison: 16 bit quantization (which is CD-quality) has $err\simeq-96$ dB, 8 bit quantization has $err\simeq-50$ dB. The tests were conducted with the PD-patch available at tommsch.com/science.php in real-time over a time of about 20 s with double precision floating point numbers and a sample rate of $44.1$ kHz. We used the following stationary and non-stationary test signals

• •

white: White noise

• •

sine X: A pure sine tone with X Hz

• •

const: A constant signal

• •

clicks: Clicks with a spacing of 1 s

• •

beet: Beethoven - Piano Sonata op 31.2, length 25 s.

• •

speech: A man counting from one to twenty, length 20 s.

• •

fire: A firework, length 23 s.

• •

atom: A sample of sparse synthesised warped Gabor atoms which were also used for that specific test run.

Since our method would lead to perfect reconstruction if the windows were time-shift invariant, the behaviour of the algorithm for stationary signals over a long time is of greatest interest. The constant signal is of interest since it is the one with the lowest possible frequency and our algorithm may bear problems with low frequencies due to the necessary cutting. The clicks represent the other extreme point of signals. The atoms are of interest since they shows whether our algorithm has the ability to reproduce its atoms with high quality.

For the warping map, we used a function which is exponentially increasing between two frequencies $f_{in}$ and $f_{out}$, namely $\theta(f)=f_{0}2^{-f/k}$ where $f_{0},k\in\mathbb{R}$ are constant parameters which can be set inside the PD-patch. For all the tests we set $f_{0}=12$ and $k=36$. Outside of $\pm f_{out}$ and between $\pm f_{in}$ the map is linear. The function attains exactly Nyquist frequency, i.e. at $\theta((q_{\sup}-1)b)=\text{SR}/2$. The frequencies $f_{in}$, $f_{out}$ are chosen such, that the resulting map is $C^{1}$. See Figure 4 for a plot of $\theta^{-1}$.

4.1 Raised Cosine Window

As proposed in [9] we first performed tests with a raised cosine window.

[TABLE]

with Fourier-Transform

[TABLE]

where T is the total duration of the window, $R\geq 3$ is an integer, $f_{0}=\frac{3}{2\text{T}}$, $a=\frac{3}{2Rf_{0}}$, $b=\frac{2f_{0}}{R}$. See [9, Section 4] for an ansatz about how to compute these parameters. This window has a very slow decay in the frequency domain after warping:

[TABLE]

This means that either the warped windows are not bounded anymore (i.e.: $\tilde{\tilde{\varphi}}_{q}\notin L^{\infty}(\mathbb{R})$) or that they are discontinuous, which is clearly visible in the warped windows, a plot of one window can be seen in Figure 5. Hence the computation of the warped windows with the IDFT bears numerical errors. Also the test results with this window were suboptimal, yielding an error between $-50$ dB and $-60$ dB, depending on the chosen parameters and input signal.

Changing the parameter $R$ or the parameter $C_{d}$ had a significant effect on the error. In Figure 6 on can see the influence of $R$. This can be expected, since the essential bandwidth, in which the phase is linearised, depends on that parameter. Since oversampling is directly proportional to the computational load (as computed in  (30)), for real-time applications there is a natural upper bound. $C_{d}\approx 2,5$ provided good results. Too small and too big numbers for $C_{d}$ both gave worst results.

On the contrary, changing the parameter $b$ while preserving the oversampling factor $R$ and the parameter $a$ had no effect, as long $C_{b}\leq 2$.

The relation between the window length after cutting and the error was nearly proportional to the value of $\sum_{q}\text{T}_{q}$ in a certain range, see Figure 8. From this observation one can determine a suitable cutting parameter $C_{\text{cut}}$. Changing the parameter $\text{T}_{\max}$ has clearly only an influence on the $err$ for low frequencies.

The parameter $C_{\text{T}_{c}}$ had no big influence as long it was about $\text{T}_{c}$ was about twice the window length $\text{T}_{h}$, see Figure 7

The warping map has a very big influence on the error. At points where the map is not smooth the error for these frequencies is order of magnitudes higher. This can be seen in Figures 9 and 10 for the sine with 30 Hz signal. At this point, our used warping map is only $C^{1}$. For a $C^{0}$ warping map the error was again 10 dB higher.

The table in Figure 9 shows selected numerical results with well chosen parameters for the raised cosine window. All values in dB, values above $-60$ dB and below $-70$ dB are coloured. The number $4\sum\text{T}_{q}/d_{q}$ denotes the average computational average complexity (see above).

4.2 Gaussian Window

In order to obtain a window with proper decay in the frequency domain after warping, we used a Gaussian window. This window does not overlap-add to one. Hence higher overlap is necessary to minimize the deviation from one. The warped windows were still fast decaying in the time domain, resulting in the possibility to cut them much shorter then the warped raised-cosine windows which compensated the high computational load due to the high overlap.

Our used Gaussian window and its Fourier Transform is

[TABLE]

where T is the total duration of the window, $R\geq 2$ the overlap factor is an integer, $f_{0}=\frac{R}{2\text{T}}$, $a=\frac{\text{T}}{R}$, $b=\frac{1}{aR}$ and $C\simeq 0.893249\cdots$ is a constant factor used to approximate the overlap-add to one condition.

This window has a fast enough decay to ensure that the warped windows in the frequency domain still are in $L^{1}(\mathbb{R})$ and hence their Inverse Fourier transformed (i.e. the warped windows in the time domain) are bounded and continuous. This window leads to significantly better results down to $-96$ dB which is the limit when using PD’s single precision numbers. The influence of the parameters on the error was the same as for the raised-cosine window.

The tables in Figure 10 show selected numerical results.

4.3 Comparison of these two windows

For the raised cosine a high number of bands must be used to achieve a small error. In Figure 8 the tests of the red dots are conducted using 92 bands, the tests with yellow dots were with 229 bands. Since the windows have a bad decay in the frequency domain, a high essential bandwidth has to be used too, which entails big overlap in time. and hence a very high average computational load, in our example $33.6$k floating point operations per sample. If one uses similar parameters to the ones used for the Gaussian window in Figure 10, the error is in the range of $-36$ dB.

Gaussian windows on the other hand do not overlap add to one and hence are not dual to themselves. Hence we at first used a very high overlap to minimize the deviation from 1, which turned out to be not necessary later. The fast decay in the time domain allows to cut the windows much shorter than the raised cosine window. This decreases the computational load, in our example only 7.3k flops per sample in average. In Figure 8 one can see that for the Gaussian window, with proper parameters, the error is in the range of $-90$ dB.

5 Conclusions

We have introduced a novel, flexible, easy way to construct frames starting from a classical Gabor frame. Tests show, that even in the painful case, where perfect time realignment of the components is not guaranteed and hence our method does not lead to perfect reconstruction, the error can be made as small as the accuracy of single precision floating point numbers for a wide range of signals. This does not prove that the error is small for all signals, but gives a good estimation of the error to be expected. The transform is working in real time.

We are going to implement maps that can be arbitrarily defined, e.g., by means of interpolation of a selected number of points. We will also implement Gabor Multipliers [17] in the redressed warped Gabor expansion (partially done already with the external winarym̃issing).

Since we already implemented this method as a Pure Data external, it is ready to use for audio-applications. The whole external explained in detail as well as the source code can be found online at  tommsch.com/science.php and in [18].

6 Acknowledgements

Many thanks to the great number of suggestions by the great number of anonymous reviewers on how to improve the paper.