format: text

text:
# Joint Time–Frequency Analysis: a tool for exploratory analysis and filtering of non-stationary time series

## 1 Introduction

The analysis of time series has always been an issue of central interest to astronomers. Indeed, the most natural way to get some insights on the characteristics of a physical system consists in studying its behavior over time. Unfortunately, maybe because until few years ago the most interesting objects known to change their luminosity were periodic/semiperiodic stars and multiple star systems, time series analysis became synonymous of research of periodicities and power-spectrum became the tool for signal analysis. Only in the eighties people began to realize that various astrophysical objects, principally extragalactic sources, can show very complex and unpredictable time evolutions (see, for example, Miller & Wiita 1991; Duschl et al. 1991). In spite of that, although inadequate, power-spectrum continued to be the preferred approach by most of researchers. Such a predilection can be explained with the substantial failure of the methods proposed as alternatives (Vio et al. 1992). The main limitations of such methods are two:

1. The implicit assumption of stationarity of the processes underlying the signals. In other words, an observed signal is assumed to be characterized by statistical properties that do not change with time. However, real signals in their vast majority are nonstationary.
2. The implicit assumption that it is possible to obtain the equations describing the evolution of a given system solely through the analysis of an experimental signal. In reality, if one takes into account that a given time series is only the one-dimensional output of a system which very often can be described only in terms of a set of ordinary and/or partial differential equations (in other words a time series is the one-dimensional projection of a more complex dynamics), it is easy to understand that this position is not well founded.

The consequence of these points is that if not carried out in a well defined physical context, signal analysis can be expected to be useful only for explorative purposes (e.g. feature detection) and/or for filtering tasks (e.g. removal of unwanted components as, for example, noise or spurious contributions). Furtehrmore, in order to obtain meaningful results, the nonstationarity of signals must be explicitly considered in the analysis.

## 2 Data exploration

### 2.1 Limits of the classical approach

It is well known that the power-spectrum $`\mathrm{PS}(\nu )`$ of a continuous and integrable (possibly complex) signal $`x(t)`$ is defined as

$$\mathrm{PS}(\nu )=|\mathrm{FT}(\nu )|^2,$$
(1)
where the frequency range is the full real line, i.e. $`\mathrm{}<\nu <+\mathrm{}`$, and

$$\mathrm{FT}(\nu )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x(t)e^{j2\pi \nu t}𝑑t$$
(2)
is the Fourier transform of $`x(t)`$. $`\mathrm{PS}(\nu )`$ can be considered as a measure of the similarity between signal $`x(t)`$ and the complex sinusoid $`\mathrm{exp}(j2\pi \nu t)`$. Indeed, if $`\mathrm{FT}(\nu )`$ is again integrable, the Fourier inversion formula reads

$$x(t)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\mathrm{FT}(\nu )e^{j2\pi \nu t}𝑑\nu ,$$
(3)
i.e., the signal can be written (uniquely) as a kind of superposition of pure frequencies. Therefore, the larger $`\mathrm{PS}(\nu )`$, the more important is the contribution of the pure frequency $`\nu `$ within $`x(t)`$.

This is the reason why $`\mathrm{PS}(\nu )`$ is so useful in the search of periodic components. However, this is also its more serious limit in the analysis of non-periodic signals (Feichtinger & Strohmer 1998). Fig.1a shows a signal constituted by two sinusoids with amplitudes different from zero only within a short time interval:

$$x(t)=a_1(t)\mathrm{sin}(2\pi \nu _1t)+a_2(t)\mathrm{sin}(2\pi \nu _2t)$$
(4)
with $`\nu _1=0.1\text{Hz}`$, $`\nu _2=0.4\text{Hz}`$ and

$$a_1(t),a_2(t)=\{\begin{array}{cc}1\hfill & \text{if }10t50\text{ and }70t110\text{ };\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}$$
(5)
with $`t`$ expressed in seconds. The two key facts are:

* this signal clearly is not periodic and
* it does not maintain its characteristics constant over time (in poor words, it is non-stationary).

From Fig.1b it is evident that $`\mathrm{PS}(\nu )`$ is not able to extract the salient features of $`x(t)`$. In particular, the power spectrum clearly contains no information about the time intervals on which these two frequencies occur. Moreover, since $`\mathrm{PS}(\nu )0`$ also at frequencies $`\nu \nu _1,\nu _2`$, one could erroneously infer that $`x(t)`$ contains more than two sinusoidal components. The main problem is that, if a signal changes its characteristic over time, a representation limited exclusively to the frequency domain is absolutely inadequate. Indeed, the function $`\mathrm{exp}(j2\pi \nu t)`$ does not have any temporal localization: it has neither start nor end. This is an useful characteristic only in case of signals with constant statistical properties (i.e. periodic or stationary signals). Unfortunately, in nature this kind of signals is not so common. Therefore, for a good analysis of most “natural” signals it is necessary to use a representation that is able to consider also the temporal dimension. A possible solution is a representation, that we call ideal power-spectrum, $`\mathrm{IPS}(t,\nu )`$, defined as the squared amplitude of the sinusoids that at a given time instant $`t`$ constitute the signal of interest.

That is shown in Fig.2 where it is possible to realize that with $`\mathrm{IPS}(t,\nu )`$, contrary to $`\mathrm{PS}(\nu )`$, the characteristics of the signal in Fig.1a are very well defined; it is clear that such a signal consists of two pure sinusoids having different frequencies but amplitudes which evolve in the same way.

### 2.2 Time–frequency representation of continuous signals

In the previous section we have seen that in case of non-stationary signals a joint time-frequency representation (JTFR) might be a good analysis tool. However, the JTFR in Fig.2 has been obtained directly from the model with which the time series has been simulated. Unfortunately, in an experimental context the dynamical model underlying a given time series will usually be unknown. Consequently, the problem consists in the way such a JTFR can be estimated from a set of available data. The simplest solution is given by the so called Fourier Spectrogram, $`\mathrm{FS}(t,\nu )`$,

$$\mathrm{FS}(t,\nu )=|\mathrm{STFT}(t,\nu )|^2,$$
(6)
where $`\mathrm{}<t<+\mathrm{},\mathrm{}<\nu <+\mathrm{}`$, and

$$\mathrm{STFT}(t,\nu )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}w^{}(t\tau )x(\tau )e^{j2\pi \nu \tau }𝑑\tau ,$$
(7)
(symbol “” stands for complex conjugation), is the so called short time fourier transform. Here

$$\begin{array}{ccc}w(t)\hfill & 0\hfill & \text{for }T/2tT/2;\hfill \\ & 0\hfill & \text{otherwise}\hfill \end{array}$$
(8)
is a window function usually indicated with the term of analysis function. The rationale is that, for signals $`x(t)`$ which maintain their harmonic components unchanged at least approximately within a time interval of length $`T`$ and centered at $`t_o`$, $`\mathrm{FS}(t_o,\nu )`$ can be considered an estimate of $`\mathrm{IPS}(t_0,\nu )`$. Although effective, this methods is characterized by two drawbacks:

* since $`\mathrm{FS}(t,\nu )`$ is estimated via the Fourier transform of only a segment of the available signal, it suffers a degradation of the frequency resolution inversely proportional to the length of $`w(t)`$. For example, if $`w(t)`$ is the rectangular window

$$w(t)=\{\begin{array}{cc}1\hfill & \text{for }T/2tT/2;\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}$$
(9)
it is easy to show that

$$\mathrm{STFT}(t_o,\nu )=\mathrm{FT}(\nu )\mathrm{sinc}_T(t_o,\nu ),$$
(10)
where the symbol “$``$” stands for convolution, and

$$\mathrm{sinc}_T(t_o,\nu )=\frac{\mathrm{sin}(T\pi \nu )}{\pi \nu }e^{j2\pi \nu t_o}.$$
(11)
This means that, in $`\mathrm{FS}(t_o,\nu )`$, the narrowest structure along the frequency direction is charaterized by a frequency support equal to $`2/T`$. Such value corresponds to the full width at zero intensity of the main lobe in the graphical representation of $`|\mathrm{sinc}_T(t_o,\nu )|^2`$;
* since $`\mathrm{FS}(t,\nu )`$ is calculated by using all the points in a interval of length $`T`$ centered at $`t`$, it suffers also a degradation of the resolution along the time dimension which is proportional to the length of $`w(t)`$. Indeed, in case of a signal constituted by an impulse at time $`t_0`$, and again with a rectangular $`w(t)`$, $`\mathrm{FS}(t,\nu )0`$ when $`t_oT/2tt_o+T/2`$. In other words, the narrowest structure along the time direction is charaterized by a time support equal to $`T`$.

From these two points it is easy to understand that the main limitation of $`\mathrm{FS}(t,\nu )`$ lies in the impossibility to increase simultaneously the time and the frequency resolutions. In other words, there is a trade-off: if $`w(t)`$ is chosen to have good time resolution (i.e., smaller $`T`$), then its frequency resolution must be deteriorated and vice versa. In principle, a possible way to mitigate this occurrence is the choice of a window $`w(t)`$ for which the degradation of the time-frequency resolution is reduce to a minimum. In this regard, it has been proved (e.g., Qian & Chen 1996) that, if the time-frequency resolution of a representation is measured by its energy concentration in the time-frequency domain, then the Gaussian function

$$w(t)=\left(\frac{\alpha }{\pi }\right)^{\frac{1}{4}}e^{\frac{\alpha }{2}t^2}$$
(12)
represents an “optimal” choice. Without any a priori information on the spectral characteristics of the signal, $`\alpha `$ is a free parameter with which it is possible to balance between an increased time or frequency resolution according to the current necessities. In reality, if $`w(t)`$ is a smooth function, well localized around zero, it will serve well as a window for the above purpose, independently from the precise analytic form. In this regard note that, because of its lack of smoothness, the rectangular window, used above only for ease of discussion, does not represent a good choice because of its very poor frequency resolution.

### 2.3 Time-frequency representation of discrete signals

In practical situations a signal is available only in discrete form, say $`x[k]`$, obtained by sampling $`x(t)`$ on a regular grid of $`L_s`$ time instants $`k=0,1,2,\mathrm{},L_s1`$. Consequently, $`\mathrm{FS}(t,\nu )`$ and $`\mathrm{STFT}(t,\nu )`$ have to be modified to their discrete variants:

$`\mathrm{DFS}[m,n]`$ $`=`$ $`|\mathrm{DSTFT}[m,n]|^2;`$ (13)
$`\mathrm{DSTFT}[m,n]`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{L_s1}{}}}w^{}[km]x[k]e^{j2\pi nk/L_s},`$ (14)
with $`m`$, $`n`$, $`0m,nL_s1`$, representing the discrete time and frequency indices, respectively.

Fig.3 shows the representation $`\mathrm{DFS}[m,n]`$ of a time series obtained by sampling the signal in Fig.1 with a time step of one second: in spite the loss of resolution described above, the salient characteristics of the signal remain well defined.

### 2.4 Some practicalities

#### 2.4.1 Digital implementation

The calculation of $`\mathrm{DFS}[m,n]`$ requires some technical adjustments. The first one regards the edge effects, due to the finite length of the signal, that do not permit the calculation of this representation for $`m<L_w/2`$ and $`m>L_s1L_w/2`$, with $`L_w`$ the length of $`w[k]`$. Usually this problem is solved through one of the following three methods:

* by assuming $`x[k]`$ and $`w[k]`$ to be periodic with period $`L_s`$ (NB. since usually $`L_w<L_s`$, $`w[k]`$ has to be zero padded);
* by zero padding one of the ends of $`x[k]`$ with $`L_w`$ data and assuming the new time series to be periodic with period $`L_s+L_w`$ (for an efficient implementation of this method see Qian & Chen 1996);
* by imposing that $`0mM_T`$, with $`M_T=L_sL_w+1`$ the number of possible shifts of $`w[k]`$ inside the signal length (for an efficient implementation of this method see Munk 1997, 2001);

If not stated otherwise, in the following the discussion will be developed according to the first approach.

A more serious drawback is represented by the fact that, by increasing the length of $`x[k]`$, the corresponding $`\mathrm{DFS}[m,n]`$ becomes quickly an huge matrix. Indeed, already $`L_s=1000`$ will result in a matrix with $`10^6`$ elements. A viable solution is to cut the signal in shorter segments and then to analyze them separately. An alternative comes out from the observation that $`\mathrm{DFS}[m,n]`$ is an highly redundant representation; it is composed by $`L_s^2`$ elements that, however, have been obtained from only $`L_s`$ data. Since $`L_s`$ is the dimension of the signal space, at most $`L_s`$ elements can be linear independent. This means that, in principle, it is possible to recover arbitrary signals if samples of $`\mathrm{DFS}[m,n]`$ for at least $`L_s`$ coordinates are available. According to this observation, Eqs.(13),(14) could be modified in

$$\mathrm{SDFS}[m,n]=|a_{mn}|^2;$$
(15)
$`a_{mn}`$ $`=`$ $`\mathrm{DSTFT}[m\mathrm{\Delta }M,n\mathrm{\Delta }N]`$
$`=`$ $`{\displaystyle \underset{k=0}{\overset{L_s1}{}}}w^{}[km\mathrm{\Delta }M]x[k]e^{j2\pi nk\mathrm{\Delta }N/L_s}`$
where $`0m<M=L_s/\mathrm{\Delta }M`$, $`0n<N=L_s/\mathrm{\Delta }N`$, with $`\mathrm{\Delta }M`$ and $`\mathrm{\Delta }N`$ integer divisors of $`L_s`$ such that $`L_s/(\mathrm{\Delta }M\times \mathrm{\Delta }N)1`$. In such a signal representation the coefficients $`a_{mn}`$ are called the Gabor coefficients and the equality (2.4.1) goes by the name Gabor representation (GR).

Given the above arguments, in order to avoid “unnecessary” redundancy, one could be tempted to use a $`\mathrm{SDFS}[m,n]`$ containing only $`L_s`$ elements. In reality, numerical stability of the expansion and easiness of the analysis (Daubechies 1992, and see below), make it necessary to work with redundant representations. The degree of redundancy is expressed through the so called oversampling rate, say $`\alpha `$,

$$\alpha =\frac{\text{number of elements of }\mathrm{SDFS}[m,n]}{\text{length of }x[k]},$$
(17)
that will take values in the range $`[1,L_s]`$. When $`\alpha =1`$, $`\mathrm{SDFS}[m,n]`$ is said to be sampled at a critical rate, whereas when $`\alpha >1`$ it is said oversampled. Using a full short-time Fourier transform corresponds to the case $`\alpha =L_s`$.

#### 2.4.2 Which $`L_w`$ for $`w[k]`$?

One important step in the construction of JTFR’s, such as $`\mathrm{SDFS}[m,n]`$, is the choice of the analysis function $`w[k]`$. Although the exact shape of this function is not so critical, the value of $`L_w`$ can severely influence the final result. In fact, in order to be meaningful, a $`\mathrm{SDFS}[m,n]`$ has to reflect both the instantaneous frequency content and the (possible) non-stationarity of the signal. However, very often these two requirements conflicts one with the other. The reason is that a large $`L_w`$ is preferable when the frequency content is of interest and shows little changes, whereas the time evolution of a signal can be better tracked with a small $`L_w`$. Again, we are forced to decide within a trade-off.

Despite its importance and its long-time use in signal processing, surprisingly this subject has found little attention in current literature (however, see Stankovic & Katkovnik 1999; Stankovic 2001, and references therein). A possible explanation is that this problem is not so well defined. Indeed, the concept of “optimal” $`L_w`$ is strictly linked to the problem at hand. For example, if one is interested in tracking the “long term” evolution of a signal rather than possible short “bursts”, then larger $`L_w`$’s represent the “best” choice. The contrary is true if the short term evolution is of interest. However, in certain applications (e.g. automatic signal detection) it should be useful to have an “objective” method able to provide a “good” $`L_w`$ without the intervention of an operator. Of course, that requires the adoption of an optimality criterion. At this regards, a possible choice is that $`L_w`$ provides the representation with the “best” global time-frequency resolution. In this case the following procedure is useful:

* a set, say $`\mathrm{\Lambda }`$, of integer values for $`L_w`$ is fixed;
* the $`\mathrm{SDFS}[m,n]`$’s corresponding to the $`L_w\mathrm{\Lambda }`$ are calculated;
* for each of these $`\mathrm{SDFS}[m,n]`$ the two-dimensional autocorrelation function $`\mathrm{ACF}[k,l]`$ is computed:

$$\begin{array}{c}\mathrm{ACF}[k,l]=\text{ }\hfill \\ \hfill =_{m=0}^{M1}_{n=0}^{N1}\mathrm{SDFS}[m,n]\mathrm{SDFS}[m+k,n+l]\end{array}$$
(18)
with $`M<k<M`$ and $`N<l<N`$.
* the value of $`L_w`$ to which corresponds the “sharpest” $`\mathrm{ACF}[m,n]`$ is selected as the searched length.

The idea behind such a procedure is that the concept of “best” global resolution implies a representation maximally concentrated in the time-frequency plane. Indeed, a value of $`L_w`$ too large or too small will results in a bad resolution along the time or frequency dimension and therefore in a representation that tends to spread in the time-frequency domain. The use of $`\mathrm{ACF}[k,l]`$ stems by the observation that the larger the spread of the representation, the more similar (correlated) the values are of the adjacent elements in $`\mathrm{SDFS}[m,n]`$. This will reflect in $`\mathrm{ACF}[k,l]`$’s with significant values also for lags $`k,l`$ different from zero.

The only open question in this procedure is the measure of the “sharpness” of $`\mathrm{ACF}[k,l]`$. At this regards there are various possibilities. In our numerical experiments we have seen that the simple count of the elements of $`\mathrm{ACF}[k,l]`$ with values larger than a given threshold, say $`50\%`$ of the maximum, is able to provide already satisfactory results. This method for the estimation of $`L_w`$ works very well even in the situation of very noisy data (see Figs.4 and 5). The only concern regards the calculation of $`\mathrm{ACF}[k,l]`$ that, if carried out on the basis of its definition (18), can be computationally very expensive. Indeed, for saving computing time, this operation has to be carried out in the Fourier domain via the inverse double Fourier transform of the bidimensional power-spectrum of $`\mathrm{SDFS}[m,n]`$. Although very fast, with this approach it is implicitly assumed that $`\mathrm{SDFS}[m,n]`$ is periodic with regards to both indices; this fact introduces edge effects. Fortunately, they can be avoided via zero padding; if interested in the correlation for lags as large as $`\pm K`$ along both dimensions, then both ends of $`\mathrm{SDFS}[m,n]`$ must be zero padded in such a way to obtain a matrix with dimensions $`(M+K)\times (N+K)`$.

### 2.5 How to improve the resolution of the time-frequency representations?

The main limitation of $`\mathrm{SDFS}[m,n]`$ is its bad resolution. In many practical situations that is not a so serious problem. However, in certain applications (e.g. tracking of quasi-periodic sources) more precise results can be required. This is a subject that has been extensively treated in current literature. Here three methods are shortly presented that are able to provide interesting results. We warn, however, that this list is by no means complete and other important approaches are available in current literature as, for instance, the reduced interference distributions (Cohen 1995).

#### 2.5.1 Wigner-Ville distribution

The loss of resolution suffered by the $`\mathrm{FS}(t,\nu )`$ is due to the fact that, for each time instant $`t`$, $`\mathrm{STFT}(t,\nu )`$ is calculated by using only a segment of the signal. In principle, if such windowing were avoided, it could be possible to obtain a JTFR with a much better resolution. Indeed such a JTFR does exist, it is named Wigner-Ville Distribution, and comes out from the problem of defining a representation, say $`\mathrm{WVD}(t,\nu )`$, with the property that:

$`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{WVD}(t,\nu )𝑑\nu `$ $`=`$ $`|x(t)|^2;`$ (19)
$`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{WVD}(t,\nu )𝑑t`$ $`=`$ $`\mathrm{PS}(\nu ).`$ (20)
In other words, the marginals with respect to frequency and time have to provide, respectively, the squared signal and the corresponding power-spectrum. The form of this function, that however is not the only one sharing properties (19) and (20), is given by (Flandrin 1999; Groc̈henig 2001)

$$\begin{array}{c}\mathrm{WVD}(t,\nu )=\hfill \\ \hfill =\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x(t+\tau /2)x^{}(t\tau /2)e^{j2\pi \nu \tau }d\tau .\end{array}$$
(21)
In reality, the use of this JTFR has been rather limited because of two drawbacks:

* $`\mathrm{WVD}(t,\nu )`$ can go negative. In fact, except very particular cases, e.g., $`x(t)`$ given by a Gaussian function, it will go negative for all signals. The non-positive property prohibits an energy interpretation of this representation;
* the quadratic nature of $`\mathrm{WVD}(t,\nu )`$ determines the occurrence of non-linear effects when dealing with multicomponents signals.

This last point is particularly troublesome. In fact, if $`x(t)=g(t)+r(t)`$, it happens that

$$\begin{array}{c}\mathrm{WVD}_x(t,\nu )=\hfill \\ \hfill =\underset{\mathrm{Auto}\mathrm{Terms}}{\underset{}{\mathrm{WVD}_g(t,\nu )+\mathrm{WVD}_r(t,\nu )}}+2\underset{\mathrm{Cross}\mathrm{Terms}}{\underset{}{Re\{\mathrm{WVD}_{g,r}(t,\nu )\}}}\end{array}$$
(22)
where

$$\begin{array}{c}\mathrm{WVD}_{g,r}(t,\nu )=\hfill \\ \hfill =\underset{\mathrm{}}{\overset{+\mathrm{}}{}}g(t+\tau /2)r^{}(t\tau /2)e^{j2\pi \nu t}d\tau .\end{array}$$
(23)
In general, the $`\mathrm{WVD}(t,\nu )`$ corresponding to a signal containing $`N`$ components will be characterized by $`N`$ auto-terms, that constitutes the desired information, and $`0.5\times N\times (N1)`$ real cross-terms which can be considered a sort of unwanted interference disturb. The worst point is that the intensity of such cross-terms is larger than those of the corresponding auto-terms. As shown in Fig.6, that can mean unreadable results. In spite of these problems, $`\mathrm{WVD}(t,\nu )`$ has very interesting properties not shared by most of the other representations. For this reason, in these last years, one of the most important subject in the field of time-frequency analysis has been the research of methods able to reduce the importance of the interference terms with the minimum deterioration of such properties.

#### 2.5.2 The reassignment method

It is possible to show (Flandrin 1999) that the spectrogram $`\mathrm{FS}(t,\nu )`$ is the 2D-convolution of the Wigner-Ville distribution of the signal $`x(t)`$, say $`\mathrm{WVD}_x(t,\nu )`$ with the Wigner-Ville distribution, say $`\mathrm{WVD}_w(t,\nu )`$ of the window function $`w(t)`$:

$$\begin{array}{c}\mathrm{FS}(t,\nu )=\hfill \\ \hfill =\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\mathrm{WVD}_x(s,\xi )\mathrm{WVD}_w(ts,\nu \xi )dsd\xi .\end{array}$$
(24)
From this equation it is possible to see that the reason why $`\mathrm{FS}(t,\nu )`$, contrary to $`\mathrm{WVD}(t,\nu )`$, does not present the interference terms is due to the fact that such terms are smoothed out by the convolution. Such operation, however, deteriorates the time-frequency resolution.

A closer look at Eq.(24) shows that $`W_w(ts,\nu \xi )`$ delimits a time-frequency domain at the vicinity of the $`(t,\nu )`$ point, inside which a weighted average of the $`\mathrm{WVD}_x`$ values is performed. This is the reason of the bad time-frequency resolution of $`\mathrm{FS}(t,\nu )`$. Indeed, the key point of the reassignment principle is that these values have no reason to be symmetrical distributed around $`(t,\nu )`$, which is the geometrical center of this domain; their average should not be assigned at this point, but rather at the center of gravity of this domain, which is much more representative of the local energetic distribution of the signal. In practice, with the reassignment approach one moves each value of the spectrogram computed at any point $`(t,\nu )`$ to another point $`(\widehat{t},\widehat{\nu })`$ which is the center of gravity of the signal energy distribution around $`(t,\nu )`$:

$$\widehat{t}(t,\nu )=\frac{{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}s\mathrm{WVD}_w(ts,\nu \xi )\mathrm{WVD}_x(s,\xi )𝑑s𝑑\xi }{{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{WVD}_w(ts,\nu \xi )\mathrm{WVD}_x(s,\xi )𝑑s𝑑\xi }$$
(25)
$$\widehat{\nu }(t,\nu )=\frac{{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\xi \mathrm{WVD}_w(ts,\nu \xi )\mathrm{WVD}_x(s,\xi )𝑑s𝑑\xi }{{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\mathrm{WVD}_w(ts,\nu \xi )\mathrm{WVD}_x(s,\xi )𝑑s𝑑\xi }$$
(26)
and thus leads to a reassigned spectrogram, whose value at any point $`(t^{},\nu ^{})`$ is the sum of all the spectrogram values reassigned to this point:

$$\begin{array}{c}\mathrm{FS}_r(t^{},\nu ^{})=\hfill \\ \hfill =\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\mathrm{FS}(t,\nu )\delta (t^{}\widehat{t}(t,\nu ))\delta (\nu ^{}\widehat{\nu }(t,\nu ))𝑑t𝑑\nu \end{array}$$
(27)
Fig.7 show $`\mathrm{FS}_r[k,l]`$ for the data shown in Fig.1. A detailed description of the algorithm and of its digital implementation can be found in Chassande-Mottin (1998).

#### 2.5.3 Gabor spectrogram

Recently, Qian & Chen (1996) have presented an interesting technique for improving the time-frequency resolution of $`\mathrm{FS}(t,\nu )`$. They start from the consideration that each cross-term in $`\mathrm{WVD}(t,\nu )`$ is located midway the pair of corresponding auto-terms and tends to highly oscillate in both time and frequency directions. On the other hand, useful properties, such as the time marginal, frequency marginal and instantaneous frequency, are obtained by averaging $`\mathrm{WVD}(t,\nu )`$. For example, the mean instantaneous frequency $`\nu _{\mathrm{inst}}`$, a quantity that is commonly used for characterizing non-stationary signals, can be evaluated via:

$$\nu _{\mathrm{inst}}=\frac{{\displaystyle \nu \mathrm{WVD}(t,\nu )𝑑\nu }}{{\displaystyle \mathrm{WVD}(t,\nu )𝑑\nu }}.$$
(28)
These observations suggest that, if $`\mathrm{WVD}(t,\tau )`$ can be decomposed as the sum of two-dimensional (time and frequency) localized harmonic functions, then we could use the low-order harmonic terms to delineate the time-dependent spectrum with reduced interference. High-order harmonics that introduce high oscillation have relatively small averages and thereby have negligible influence on the useful properties. The signal energy and useful properties are mainly determined by a few low-order harmonic terms. On these bases, they propose the so called Gabor spectrogram $`\mathrm{GS}_D`$ representation

$$\mathrm{GS}_D[k,l]=\underset{d=0}{\overset{D}{}}P_d[k,l]$$
(29)
where

$$\begin{array}{c}P_d[k,l]=2_{A_d}a_{mn}a_{mn}^{}\times \hfill \\ \hfill \mathrm{exp}\{\frac{\left[k\mathrm{\Delta }M(m+m^{})/2\right]^2}{2\sigma _k^2}\frac{\left[l\mathrm{\Delta }N(n+n^{})/2\right]^2}{2\sigma _l^2}\}\times \\ \hfill \mathrm{exp}\{\frac{j2\pi }{L_w}[k(nn^{})\mathrm{\Delta }N+\\ \hfill +(l\frac{n+n^{}}{2}\mathrm{\Delta }N)\mathrm{\Delta }M(m^{}m)]\},\end{array}$$
(30)
where

$$A_d=\{(m,m^{}),(n,n^{});|mm^{}|+|nn^{}|=d\},$$
(31)
$$\sigma _k=\frac{T_1}{2\sqrt{\pi }},\sigma _l=\frac{L_w}{T_12\sqrt{\pi }},$$
(32)
with $`0m,m^{}M_{tot}1`$, $`N/2n,n^{}N/21`$, $`0k<L_s`$, and $`L_w/2l<L_w/21`$. Here, $`T_1`$ is a parameter whose value can be chosen $`\sqrt{\mathrm{\Delta }MN}`$ and $`a_{mn}`$ are the Gabor coefficients calculated via Eq.(2.4.1) that here we rewrite in a slightly modified form

$$\begin{array}{c}a_{mn}=e^{j2\pi mn\mathrm{\Delta }M/N}\times \hfill \\ \hfill _{k=0}^{L_w1}w^{}[k]x[k+m\mathrm{\Delta }M]e^{j2\pi kn/N},\end{array}$$
(33)
$`0mM_T=[(L_sL_w)/\mathrm{\Delta }M+1]`$, $`0n<N`$. It is possible to show that the larger the value of $`D`$ the better is the resolution of the representation but, at the same time, the higher is the contribution of the high-order harmonics in the decomposition of $`\mathrm{WVD}(t,\nu )`$. That means that there is a trade-off between time-frequency resolution and the importance of the interference terms. Usually $`D=3`$,$`4`$ is sufficient for many practical applications. For the time series of Fig.3, $`\mathrm{GS}_4[k,l]`$ provides results very similar to those of Fig.7. A detailed description of the algorithm and of its digital implementation can be found in Munk (1997, 2001).

## 3 Time-variant filtering of discrete signals

Although data exploration is an interesting application of time-frequency analysis, the issue where such an approach shows all its potency is signal filtering (Matz & Hlawatsch 2002; Hlawatsch & Matz 2002). Indeed, given the capability of a JTFR in characterizing a non-stationary signal, we may expect that a filtering technique, developed within this framework, be able to adapt its action to the instantaneous properties of the signal itself. Since it is much easier to treat this subject in the context of the linear representations ($`\mathrm{WVD}`$, $`\mathrm{GS}_D`$ and $`\mathrm{FS}_r`$ are bilinear representations carrying no phase information in a form easy to use), we shall concentrate in our discussion on the use of the Gabor transform.

### 3.1 Some preliminary notes

We have seen that the GR of a discrete signal corresponds to map a one-dimensional array onto a two dimensional matrix. Of course, in order to benefit of this kind of representation for signal filtering purposes, it is necessary to be able to do the reverse job. It can be shown that the relationship between a discrete signal $`x[k]`$ and the corresponding GR is given by (Groc̈henig 2001)

$$x[k]=\underset{m=0}{\overset{M1}{}}\underset{n=0}{\overset{N1}{}}a_{mn}g[km\mathrm{\Delta }M]e^{j2\pi kn/N}$$
(34)
with $`0k<L_s`$ and

$$a_{mn}=\underset{k=0}{\overset{L_s1}{}}w^{}[km\mathrm{\Delta }M]x[k]e^{j2\pi kn/N}$$
(35)
with $`0m<M,0n<N`$. Here, the synthesis function $`g[k]`$ and the analysis function $`w[k]`$ have to satisfy the so called biorthogonality condition (or Wexler–Raz identity).

$$\underset{k=0}{\overset{L_s1}{}}w^{}[k]g[k+pN]e^{j2\pi kq/\mathrm{\Delta }M}=\frac{\mathrm{\Delta }M}{N}\delta _p\delta _q$$
(36)
with $`\mathrm{\Delta }N<q<\mathrm{\Delta }N,0p<\mathrm{\Delta }M`$. Pairs of functions satisfying Eq.(36) are termed dual functions or windows. Since this formula is symmetric with respect to $`g[k]`$ resp. $`w[k]`$, duality is a symmetric relation.

For oversampled GR, with $`\alpha >2L_w/(L_w+\mathrm{\Delta }M)`$, the system of linear Eq.(36) is underdetermined. In other words, $`g[k]`$ and $`w[k]`$ do not determine each other in a unique way. Hence, for a given $`g[k]`$ there is a certain freedom in the choice of the corresponding dual complementary $`w[k]`$. The solution suggested by Qian & Chen (1996) is to solve system (36) with the constraint, say $`\mathrm{\Gamma }`$, that the shapes of $`g[k]`$ and $`w[k]`$ be as close as possible (in least squares sense):

$$\mathrm{\Gamma }=\mathrm{min}\underset{k=0}{\overset{L_s1}{}}\left(\frac{w[k]}{w}\frac{g[k]}{g}\right)^2.$$
(37)
This choice turns out to correspond to the solution of the system (36) via a pseudo-inverse method. For example, by supposing $`g[k]`$ fixed and rewriting system (36) in a matrix notation,

$$\mathrm{𝐇𝐰}^{}=𝐮,$$
(38)
the standard synthesis function $`𝐰`$, i.e., the window satisfying Eq.(36) with minimal norm, is given by

$$𝐰^{}=𝐇^{}(\mathrm{𝐇𝐇}^{})^1𝐮.$$
(39)
It can be shown that the larger $`\alpha `$ the more similar are $`g[k]`$ and $`w[k]`$ (see Fig.8). In particular, when $`\alpha =L_s`$ (i.e., $`\mathrm{\Delta }M=\mathrm{\Delta }N=1`$) then the two functions are identical.

In case of a critical sampled GR (i.e. $`\alpha =1`$,) $`w[k]`$ and $`g[k]`$ are univocally determined. However, according to the Balian-Low theorem, both the functions cannot be simultaneously well localized in time and frequency (see, for example, Mallat 1998). In the practical application, the consequence is that $`w[k]`$ and $`g[k]`$ have very different shapes. This is the main reasons why in the filtering problems it is preferable to work with redundant representations (see below).

### 3.2 Time–variant signal denoising

As shown in Fig.9, the choice of $`\gamma _a`$ is an important step in a filtering procedure. Unfortunately, without substantial a priori information (typical situation in many applications), there are no so many ways to select the “optimal” value for this threshold.

In many practical situations the signal of interest is contaminated by a broad–band random noise (white–noise). The key point is that such a noise tends to spread evenly over the entire joint time-frequency domain, whereas the signal component is usually concentrated in a relatively small region (e.g., see Fig.4). This fact suggests a simple procedure, similar to the hard thresholding method introduced by Donoho for wavelets (see, for example, Donoho 1995), to filter out the noise contribution:

* determination of those indices for which the corresponding coefficients $`a_{mn}`$ may be considered as contributed by the noise. This step can be carried out in several ways according to the problem at hand. A typical solution consists in setting a threshold value, say $`\gamma _a^2`$, for the elements $`|a_{mn}|^2`$ of $`\mathrm{SDFS}[m,n]`$. The idea is that the coefficients $`a_{mn}`$ of GR, whose absolute value is smaller than $`\gamma _a`$, can be assumed unaffected by the signal contribution;
* masking (zeroing) of such coefficients in order to obtain a “cleaned” version of the coefficients, denoted by $`\stackrel{~}{a}_{mn}`$, equal to zero in areas only occupied by noise;
* synthesis of the filtered signal $`\widehat{x}[k]`$ via Eq.(34).

Although this approach is fairly simple, its application requires some care in order to avoid unpleasant side effects. In particular, one point is important:

> once the analysis function $`g[k]`$ is fixed, it is advisable that the synthesis function $`w[k]`$ be calculated via Eq.(39). In other words, the shapes of $`g[k]`$ and $`w[k]`$ have to be as close as possible.

The good reason for this is that $`g[k]`$ is often chosen to have a good localization both in time and frequency. If $`w[k]`$ would not be forced to share this property it might be badly localized along one or both of these two dimensions. The consequence is that the synthesis of $`\widehat{x}[k]`$, via the coefficient $`\stackrel{~}{a}_{mn}`$, could provide different results from what expected;

Here it is necessary to stress that, in general, the set of modified coefficients $`\stackrel{~}{a}_{mn}`$ will not constitute a valid GR representation of any signal whatever. Indeed, it happens that

$$\widehat{x}[k]\underset{m=0}{\overset{M1}{}}\underset{n=0}{\overset{N1}{}}\stackrel{~}{a}_{mn}g[km\mathrm{\Delta }M]e^{j2\pi kn/N}.$$
(40)
In particular the coefficients $`\widehat{a}_{mn}`$, obtained from $`\widehat{x}[k]`$ via transform (35), can be different from zero even for those indices $`m_0,n_0`$, for which $`\stackrel{~}{a}_{m_0n_0}=0`$. The reason is that coefficients $`a_{mn}`$ are obtained via a convolution operation and therefore no physical signal can correspond to a representation $`\stackrel{~}{a}_{mn}`$ characterized by the sharp cutoffs created by the masking procedure. In any case, it can be shown (see appendix A) that, if $`w[k]`$ and $`g[k]`$ satisfy the constraint (37), $`\widehat{x}[k]`$ is a signal whose Gabor coefficients $`\widehat{a}_{mn}`$ are as close as possible (in least squares sense) to $`\stackrel{~}{a}_{mn}`$.

A last note regards the fact that, in certain situations, the zeroing of coefficients $`a_{mn}`$ can constitute a too drastic operation. In fact, it can introduce some spurious structures in the reconstructed signals. In this case, a possible alternative is represented by the weighting of coefficients $`a_{mn}`$, according to the expected level of noise contamination, in a way similar to that suggested by Donoho (1995) (soft thresholding) in the context of wavelets.

#### 3.2.1 Practical determination of an “optimal” $`\gamma _a`$

Apart the very common assumption of Gaussianity, the typical information available about the noise is an estimate of its standard deviation. Conversely to the case of classical Fourier analysis, the statistical characterization of a pure random Gaussian process in the time-frequency domain is much more complex. The problem is that coefficients $`a_{mn}`$ are not independent. This situation suggests that the simplest (and safest) way to determine the threshold $`\gamma _a`$ is the following:

* simulation of a noise process (not necessarily Gaussian) via a pseudo-random generator;
* calculation of the corresponding $`\mathrm{SDFS}[m,n]`$ with the same modalities used for the calculation of the GR of the signal;
* estimation of a simple statistical index (e.g. empirical percentiles), describing the characteristics of such a representation, that provides the desired threshold for the $`\mathrm{SDFS}[m,n]`$ of the signal.

Although very simple, in our numerical simulations this method has proved to be rather solid (see Fig.10 to compare with Fig.9). An alternative to the last step could be the calculation of a set of filtered time series $`\widehat{x}[k]`$ corresponding to different choices of $`\gamma _a`$. After that, among these solutions it is chosen the time series that satisfies some statistical requirements as, for example, that the standard deviation of the residuals between the original and the filtered signal is equal to the nominal value of noise level. Although appealing, from our numerical experiments this approach appears to be not so stable.

In the case one lacks also the information regarding the standard deviation of the noise, the simplest solution consists in considering directly the $`\mathrm{SDFS}[m,n]`$ of the signal and then in determining the threshold interactively by delimiting, via some graphical facility, the region where the signal is supposed to give its contribution. Another possibility is the use of some statistical indices of the values of $`\mathrm{SDFS}[m,n]`$ (e.g. percentiles).

#### 3.2.2 Could the SNR be further improved?

The synthesis of $`\widehat{x}[k]`$, via the coefficients $`\stackrel{~}{a}_{mn}`$, permits to filter out the noise contribution corresponding to the coefficients $`a_{mn}`$ that have been zeroed. However, what about the contamination of the remaining coefficients? In their book, Qian & Chen (1996) suggest a method that, according to their opinion, should further improve the SNR of the final results. Their idea is very simple and consists in an iterative application of the procedure described in the previous section. In other words, once obtained the filtered $`\widehat{x}[k]`$, to this signal is applied the same procedure used for $`x[k]`$ and so on.

The most serious problem of this method consists in the lack of theoretical argumentation that guarantee its reliability. Indeed, the repeated application of the masking procedure corresponds to a contraction, namely to an operation for which at the n-th iteration

$$x_n[k]^2<x_{n1}[k]^2.$$
(41)
In other words, the total energy of the filtered signal is continuously reduced with the number of iterations. At first sight, one could accept this fact as useful since effectively the noise tends to increase the power of an observed time series. Unfortunately, such an increase regards the global characteristics of a signal whereas the procedure suggested by Qian & Chen is applied to the coefficient $`a_{mn}`$ that reflect its local properties. In other words, given an observed Gabor coefficient, the only thing one is allowed to claim is that, because of the noise contamination, it will be different from the true one. For a correct filtering, however, we should be able to decide at least whether the contaminated coefficient is larger or smaller than the uncontaminated one.

The only situation where the Qian & Chen procedure can be expected to work is when $`\gamma _a`$ is chosen too low. In situations like this one, a large fraction of the unmasked coefficient $`a_{mn}`$ could be erroneously assumed linked to the signal. In general, these coefficients are smaller, in absolute value, than the coefficients effectively reflecting the signal contribution. Consequently, the contraction connected with the Qian’s & Chen’s technique can be able to smooth out them before the signal component be too negatively affected. However, we stress that, contrary to what claimed by the two authors, the improvement is due to the removal of Gabor coefficients still related to the noise rather that a better estimate of the coefficients representing the contribution of the signal. Better results are to be expected with a more appropriate choice of $`\gamma _a`$. All that is shown in the Fig.11 to compare with Fig.9.

### 3.3 Signal separation

In certain practical situations it can happen that a given signal is the mixture of two or more components. Often it can be of interest to separate the various contributions. In this respect, it can be shown (Hopgood 2000, and appendix B) that, within the context of linear filtering, this operation is possible (to a given degree of accuracy) only when a representation of signal is available where the various components are disjoint. For example, in the Fourier domain, two components are “perfectly” separable if their power-Spectra do not overlap. Therefore, in the context of the joint time-frequency representations, it can be expected that separability of signal mixtures is possible if the components do not overlap in the time-frequency domain (see Figs.12 and 13). A simple and yet effective procedure is the same as that presented in section 3.2: each component is reconstructed by zeroing the Gabor coefficients corresponding to the remaining ones.

In order to understand the potential but also the limits of such an approach, it is necessary to understand what does it mean the fact that two signal are disjoint in the time-frequency domain. To discuss this problem, let suppose to have a signal $`x(t)`$ given by the sum of two components, say $`x_1(t)`$ and $`x_2(t)`$. If $`\mathrm{SDFS}[m,n]`$ is interpreted as an energy distribution, no overlap in the time-frequency plane implies that

$$|x(t)|^2𝑑t=|x_1(t)|^2𝑑t+|x_2(t)|^2𝑑t.$$
(42)
In other words, the energy of the observed signal is equal to the sum of the energies corresponding to the single components. This point is very important since indicates that $`x_1(t)`$ and $`x_2(t)`$ can be separated only if they do not interference and the total energy is conserved. From the statistical point of view this condition means that the two components are uncorrelated. In fact, if $`x_1(t)`$ and $`x_2(t)`$ are supposed, without loss of generality, zero mean processes, then the integrals in Eq.(42) provide the variance of the signals. However, it is well known that

$$\sigma _{x(t)}^2=\sigma _{x_1(t)}^2+\sigma _{x_2(t)}^2+2\sigma _{x_1(t)x_2(t)}$$
(43)
where the last term provides the covariance between $`x_1(t)`$ and $`x_2(t)`$. Therefore, Eq.(42) implies that $`\sigma _{x_1(t)x_2(t)}=0`$.

From the above discussion it is evident that, in the separation of signals overlapping in the time-frequency domain, things have to be expected much more complex. For example, Figs.14 shows the $`\mathrm{SDFS}[m,n]`$ corresponding to three signals, sampled with a frequency of one Hertz, that have been obtained through the following models:

$$x(t)=\mathrm{sin}[2\pi t\nu _a(t)]+\mathrm{sin}[2\pi t\nu _b(t)+\varphi _i]$$
(44)
with $`i=1,2,3`$, $`\nu _a(t)=4.0\times 10^3\times t+0.10`$, $`\nu _b(t)=3.5\times 10^3\times t+0.45`$, $`0t100`$ sec., and $`\varphi _1=0`$, $`\varphi _2=\pi /2`$, and $`\varphi _3=\pi `$ rad. In practice, the only difference lies in the phase $`\varphi `$ of the second component. The key point to note is that, although the $`\mathrm{SDFS}[m,n]`$’s are different, these signals are characterized by the same $`\mathrm{IPS}(t,\nu )`$ (see Fig.14). In other words, in this specific case the equality (42) does not hold. Of course, that is the consequence of the interference of the components present in $`x(t)`$. Unfortunately, the details of such a interference will depend not only on the amplitudes of the components at a given time instant, but also on the corresponding phases. Therefore, their separation should require that each coefficient $`a_{mn}`$ be able to provide information about four parameters. At least of some a priori information (e.g. the value of the phases), this is clearly an underdetermined problem. The situation becomes worst in case the number of overlapping signals is larger than two.

One time more, these facts suggest that in order to obtain some insights on the characteristics of a complex system, in general, the analysis of a signal, if not developed within a given physical context, is absolutely insufficient.

## 4 JTFR better than wavelets?

Joint time-frequency analysis does not constitute the only possible tool for the analysis of non-stationary signals. Nowadays the so-called Wavelet Analysis arises as a new alternative to JTFA. In the continuous case the (Scalogram) of a signal $`x(t)`$ is given by:

$$\mathrm{SCL}(t,a)=|\mathrm{WT}(t,a)|^2$$
(45)
with $`W(t,a)`$ being the continuous wavelet transform given by

$$\mathrm{WT}(t,a)=\frac{1}{\sqrt{|a|}}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x(s)\mathrm{\Psi }^{}\left(\frac{st}{a}\right)𝑑s.$$
(46)
Here $`\mathrm{\Psi }(t)`$ is a function, named wavelet, with the property that $`\mathrm{\Psi }(t)=0`$. The idea behind such an approach is that, if $`\mathrm{\Psi }(t)`$ is centered in the neighborhood of $`t=0`$ and has a Fourier transform centered at a frequency $`\nu _00`$, then in the time-frequency plane $`\mathrm{\Psi }\left[(tb)/a\right]`$ will be centered at $`t=b`$ and $`\nu =\nu _0/a`$, respectively. Consequently, $`\mathrm{SCL}(b,a)`$ will reflect the signal behavior in the vicinity of such time instant and frequency. The main difference with $`\mathrm{FS}(t,\nu )`$ is that the variable $`a`$ corresponds to a scale factor, in the sense that taking $`|a|>1`$ dilates $`\mathrm{\Psi }`$ and taking $`|a|<1`$ compress it. Therefore, in contrast to the analysis function $`w(t)`$ in JTFR which maintains a fixed duration (or width) but changes its shape due to frequency modulation, in wavelet analysis when the scale factor $`a`$ is varied $`\mathrm{\Psi }(t,a)`$ maintains its shape but changes its duration and consequently its bandwidth. The practical consequence is that, while $`\mathrm{FS}(t,\nu )`$ is characterized by a time-frequency resolution constant over the entire time-frequency domain, $`\mathrm{SCL}(t,a)`$ presents a better time resolution at high frequencies but a lower resolution at low frequencies (see Fig.15 to compare with Fig.3). At this point, the natural question arises concerning which choice, between the time-frequency or wavelet approach, is the most appropriate one in the context of the analysis of astrophysical signals.

Of course, the answer strictly depends on the problem at hand. In particular, it is useful to have a representation characterized by a non-constant time-frequency resolution when the signal of interest presents a non-stationary evolution on very different time scales as, for example, in case of non-stationary self similar (fractal) signals. On the other hand, especially in problems of exploratory data analysis, when the time behavior of a signal is dominated by one or more non-stationary processes which evolve on similar time scales (e.g. a periodic process that changes slowly its period), then it is by far better to use a representation with a constant time-frequency resolution since it permits an easier tracking of the changing characteristics of the signal. In the astrophysical applications, very often this is the situation of interest.

It is necessary to mention that, especially in the field of discrete wavelets, the above indicated limitations have been overcome in the sense that some approaches have been devised that permit a large flexibility in the choice of the resolution across the time-frequency plan (e.g., wavelet packets, local cosine bases, cf. Mallat 1998). In spite of that, in the exploratory analysis of the typical astrophysical signals the joint-time frequency approach maintains its supremacy since it is more intuitive and with better developed theoretical backgrounds as far as the interpretation is concerned. Indeed, while JTFA essentially makes use of localized sinusoids, the shape of the wavelets is much more difficult to explain from the physical point of view. Consequently, the results provided by the time-frequency approach are more amenable for a direct physical interpretation.

## 5 Practical applications

To illustrate the analysis strengths of JFTA in astronomical applications, we have applied some simple JFTA algorithms to a sequence of X-ray observations of low mass X-ray binary Sco X-1 (van der Klis, priv. comm.) made with the Proportional Counter Array (PCA) on the Rossi-XTE spacecraft (Bradt et al. 1993). The time series used, shown in Fig.16 (lower panel), is a 3 seconds subset of the available sequence that is 1500 seconds long. Sampling is regular with a time step of 125 microsecond. As Sco X-1 shows quasi-periodic oscillations (QPO) in the X-ray light emission, its light curve is characterized a periodic component that changes its period over time. Therefore, an important experimental activity consists in tracking such changes. Figs.16,17 show that even a very simple JTFR as DFS is able to provide very interesting results, whereas Figs.18, 19 show how it is possible to extract the component of interest from the observed signal. Although still preliminary and partial, this example represents a good illustration of the possibilities of the application of the technique to astrophysical (and other) time series.

## 6 Summary and conclusions

In this paper we have considered the joint time-frequency domain as a possible framework for the analysis of the astrophysical signals. The methodologies based on such a domain appear very interesting since they allow reliable analyses and filterings of non-stationary time series that are very difficult, if not impossible, to carry out in the standard Fourier domain. These methodologies appear also useful for the separation of signal mixtures.

###### Acknowledgements.

We thank Prof. M. van der Klis for the availability of the unpublished data on Sco X-1 and Prof. H. Feicthinger for his useful comments.

## Appendix A Masking of the Gabor Coefficients as Least Squares Filtering

Transforms (34) and (35) can be rewritten as matrix multiplications. In particular, the Gabor transform (35) described as a linear mapping from the signal space to the coefficient space $`𝒞^{MN}`$ takes the form

$$𝐚=𝐆_𝐰^{}𝐱$$
(47)
where $`𝐚=\mathrm{vec}\{a_{mn}\}`$ and $`𝐆_𝐠`$ is an $`(MN)\times L_s`$ matrix

$$𝐆_𝐰=\left\{\begin{array}{ccccccccccccccc}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{T}_{\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{T}_{2\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{}& & & & & & & & & & & & & & \\ \mathrm{T}_{M\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{M}_{\mathrm{\Delta }N}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{M}_{\mathrm{\Delta }N}\mathrm{T}_{\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{M}_{\mathrm{\Delta }N}\mathrm{T}_{2\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{}& & & & & & & & & & & & & & \\ \mathrm{M}_{\mathrm{\Delta }N}\mathrm{T}_{M\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{M}_{2\mathrm{\Delta }N}g& & & & & & & & & & & & & & \\ \mathrm{M}_{2\mathrm{\Delta }N}\mathrm{T}_{\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{}& & & & & & & & & & & & & & \\ \mathrm{M}_{N\mathrm{\Delta }N}\mathrm{T}_{(M1)\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \\ \mathrm{M}_{N\mathrm{\Delta }N}\mathrm{T}_{M\mathrm{\Delta }M}𝐰^{}& & & & & & & & & & & & & & \end{array}\right\}$$
(48)
where $`\mathrm{M}_q`$ and $`\mathrm{T}_p`$ are, respectively, the modulation and the translation operators

$$\begin{array}{cc}(\mathrm{M}_qw)[k]\hfill & =w[k]e^{i2\pi qk/L_s}\hfill \\ (\mathrm{T}_pw)[k]\hfill & =w[kp],\hfill \end{array}$$
(49)
and

$$(\mathrm{M}_q\mathrm{T}_pw)[k]=w[kp]e^{i2\pi qk/L_s}.$$
(50)
Since $`𝐆_𝐰`$ is a full rank matrix, from Eq.(47) it easy to show that the transform (34) can be expressed in the form

$$𝐱=𝐆_𝐰^{}𝐚,$$
(51)
where

$$𝐆_𝐰^{}=(𝐆_𝐰^{}𝐆_𝐰)^1𝐆_𝐰^{}$$
(52)
is the pseudo-inverse matrix of $`𝐆_𝐰`$ with dimensions $`L_s\times (MN)`$.

Here the key point is that, if in Eq.(47) we suppose $`𝐚`$ and $`𝐆_𝐰`$ fixed and $`𝐱`$ unknown, then we are dealing with a system of $`M\times N`$ equations in $`L_s`$ unknowns. For oversampled GR the number $`M\times N`$ is larger than $`L_s`$ and therefore the system of equations (47) is overdetermined. In this case, it can be shown that Eq.(51) delivers the corresponding least squares solution (Bjorck 1990). The same still holds also when $`𝐚`$ is not in the range of $`𝐆_𝐰`$. That happens exactly, for example, when the set of $`a_{mn}`$ does not constitutes a valid GR for $`x[k]`$. In this case Eq.(47) has to be rewritten as

$$\stackrel{~}{𝐚}\stackrel{LS}{=}𝐆_𝐰\widehat{𝐱}$$
(53)
in order to highlight that the equality has to be intended in least squares sense. However, $`\widehat{x}[k]`$ will be characterized by a GR $`\widehat{a}_{mn}=𝐆_𝐰\widehat{𝐱}`$. Consequently, Eq.(53) can be rewritten as

$$\stackrel{~}{𝐚}\stackrel{LS}{=}\widehat{𝐚}.$$
(54)
The meaning of this last expression is that, although $`\stackrel{~}{a}_{mn}`$ is not a valid GR, $`\widehat{x}[k]`$ synthesized via Eq.(51) corresponds to a signal whose GR is as close as possible to $`\stackrel{~}{a}_{mn}`$ in least squares sense.

## Appendix B Signal Separability

In this section we briefly present the conditions according to which, given an observed signal $`x(t)=d(t)+n(t)`$, the component of interest $`d(t)`$ and the unwanted component $`n(t)`$ can be separated.

A signal $`x(t)`$ may be represented on an arbitrary domain $`\lambda `$ by the following general integral transform:

$$X(\lambda )=_Tx(t)K(t,\lambda )𝑑t\lambda \mathrm{\Lambda }$$
(55)
$$x(t)=_\mathrm{\Lambda }X(\lambda )k(\lambda ,t)𝑑\lambda tT$$
(56)
where $`\mathrm{\Lambda }`$ is the region in the signal space that the representation $`X(\lambda )`$ lies in <sup>1</sup><sup>1</sup>1In case of stochstic signals, integrals are to be interpreted in a mean-square limit.. Here the functions $`K(t,\lambda )`$, called direct transform kernel, and $`k(\lambda ,t)`$, called the inverse transform kernel, are linked by the relationships:

$$_\mathrm{\Lambda }k(\lambda ,t)K(\tau ,\lambda )𝑑\lambda =\delta (t\tau )$$
(57)
$$_Tk(\lambda ,t)K(t,\widehat{\lambda })𝑑t=\delta (\lambda \widehat{\lambda })$$
(58)
For example in case of the Fourier transform $`\lambda `$ is equal to the Fourier frequency $`\nu `$, $`K(t,\nu )=\mathrm{exp}(i2\pi \nu t)`$ and $`k(\nu ,t)=K^{}(t,\nu )`$.

If signal $`d(t)`$ has to be reconstructed by linear filtering of the observed signal $`x(t)`$ (i.e., by making use of linear operators), then it has to hold (see, for example, Hopgood 2000)

$$\widehat{d}(t)=_Th(t,\tau )x(\tau )𝑑\tau ,$$
(59)
where $`h(t,\tau )`$ is the response at time $`t`$ given an impulse occurred at the filter input at time $`\tau `$. A sufficient condition for the existence of $`h(t,\tau )`$ is that it can be written in the form

$$\begin{array}{cc}h(t,\tau )\hfill & =_{\mathrm{\Lambda }_d}k(\lambda ,t)K(\tau ,\lambda )𝑑\lambda \hfill \\ & +_{\mathrm{\Lambda }_0}_{\mathrm{\Lambda }_0}H_0(\lambda ,\widehat{\lambda })k(\lambda ,t)K(\tau ,\widehat{\lambda })𝑑\lambda 𝑑\widehat{\lambda }.\hfill \end{array}$$
(60)
Here $`\mathrm{\Lambda }_d`$ is the region of the $`\mathrm{\Lambda }`$ space over which the spectral component of $`d(t)`$ does not overlap the spectral component of $`n(t)`$, whereas $`\mathrm{\Lambda }_0`$ is the spectral region where these two components overlap. $`H_0(\lambda ,\widehat{\lambda })`$ is the solution of

$$P_{dx}(\lambda ,\widehat{\lambda })=_{\mathrm{\Lambda }_0}H_0(\lambda ,\stackrel{~}{\lambda })P_{xx}(\stackrel{~}{\lambda },\widehat{\lambda })𝑑\stackrel{~}{\lambda },\lambda ,\widehat{\lambda }\mathrm{\Lambda }_0$$
(61)
where

$$P_{yz}(\lambda ,\widehat{\lambda })=\{\begin{array}{cc}Y(\lambda )Z^{}(\widehat{\lambda })\hfill & \text{if }y(t)\text{}z(t)\text{ deterministic};\hfill \\ E[Y(\lambda )Z^{}(\widehat{\lambda })]\hfill & \text{if }y(t)\text{}z(t)\text{ stochastic};\hfill \end{array}$$
(62)
is the generalized cross power-spectrum concerning to the processes $`y(t)`$ and $`z(t)`$. It can be shown (Hopgood 2000) that the reconstruction $`\widehat{d}(t)`$, as provided by Eqs.(59) and (60), is characterized by a variance $`\sigma ^2(t)=E[|\widehat{d}(t)d(t)|^2]`$ equal to

$$\sigma ^2(t)=_{\mathrm{\Lambda }_0}_{\mathrm{\Lambda }_0}P_{\sigma \sigma }(\lambda ,\widehat{\lambda })k(\lambda ,t)k^{}(\widehat{\lambda },t)𝑑\lambda 𝑑\widehat{\lambda }$$
(63)
with $`P_{\sigma \sigma }`$ the generalized power-spectrum of $`\sigma ^2(t)`$.

Here the important points are three:

* the first term in the rhs of Eq.(60) is independent of any signal. Moreover, it constitutes a so called ideal filter for $`\mathrm{\Lambda }_d`$, namely a filter which passes without distorsion all generalized frequencies components falling in $`\mathrm{\Lambda }_d`$ and rejects all others;
* the second term in the rhs of Eq.(60) depends, via Eq.(61), on the crosscorrelation function between $`x(t)`$ and $`d(t)`$ and therefore, since $`x(t)=d(t)+n(t)`$, on the crosscorrelation between $`d(t)`$ and $`n(t)`$;
* the magnitude of $`\sigma ^2(t)`$, given by Eq.(63), depends on the extension of the overlap between $`d(t)`$ and $`n(t)`$.

This last point is particularly interesting since it indicates that signals $`d(t)`$ and $`n(t)`$ can be perfectly separated only when they have disjoint supports in $`\mathrm{\Lambda }`$. Indeed, in such a case, $`\mathrm{\Lambda }_0`$ is an empty set and the three points above imply, respectively, that:

* Eq.(60) simplifies in

$$h(t,\tau )=_{\mathrm{\Lambda }_d}k(\lambda ,t)K(\tau ,\lambda )𝑑\lambda ;$$
* $`d(t)`$ and $`n(t)`$ are uncorrelated signals;
* $`\sigma ^2(t)=0`$.

In other words, the perfect separation of two signals, via a linear filtering, requires the following condition be satisfied:

> there exists some domain where the generalized spectral representation of the desired components are disjoint. Thus, the desired components must be uncorrelated processes.

All that holds also in the context of joint time-frequency analysis. The only different thing is that the function $`h(t,\tau )`$, with arguments in the time domain, has to be changed, via a Fourier transform with respect to $`\tau `$, in $`H(t,\nu )`$ which has arguments in the mixed time-frequency domain.


text:
# Universal quantum computing based on monodromy representations

## Logical gates for a quantum computer

Quantum mechanical processes give rise to new types of computation. Computational networks built out of quantum-mechanical gates provide a natural framework for constructing quantum computers. The computing capacity of such computers drastically exceeds that of traditional computer , .

The quantum analogue of the classical bit is the quantum bit or *qubit*. Just as the classical bit is represented by a system which can adopt one of two distinct states ‘$`0`$’ and ‘$`1`$’, one can define a quantum bit as follows:

Definition 1. A *qubit* is a quantum system whose state can be fully described by a superposition of two orthogonal eigenstates labeled $`|0`$ and $`|1`$.

The space of quantum states is a Hilbert space, which we denote by $``$. Therefore, a qubit is a normalized state in the two-dimensional Hilbert space $`𝐂^2`$. A general state $`|\mathrm{\Psi }`$ of the qubit is given by $`|\mathrm{\Psi }=\alpha |0+\beta |1`$ with $`|\alpha ^2|+|\beta ^2|=1`$.

The value of the qubit is the observable $`N`$ with the Hermitian operator $`N|i=i|i`$ on the Hilbert space $`𝐂^2`$, or, in the matrix representation

$$N=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right).$$

The expectation value of $`N`$ is given by

$$N=\mathrm{\Psi }|N|\mathrm{\Psi }=(\alpha ^{}\beta ^{})\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \end{array}\right)=|\beta |^2$$

Thus, $`N`$ gives the probability to find the system in state $`|1`$ if a measurement is performed on the qubit.

While the state of a classical computer can be given as the collection of distinct states of all bits in the memory and processor registers, the “state of a qubit” is a meaningless term, if the machine state is the combined state of more than one system.

Definition 2. The machine state $`|\mathrm{\Psi }`$ of an $`n`$-qubit quantum computer is given by $`|\mathrm{\Psi }=_{(d_0,\mathrm{},d_{n1})}c_{d_0,\mathrm{},d_{n1}}|d_0\mathrm{}d_{n1}`$ with $`|c_{d_0,\mathrm{}.d_{n1}}|^2=1`$.

The quantum state space $``$ is thus the tensor product of $`n`$ single qubit Hilbert spaces $`_j𝐂^2`$, i. e.

(1)
$$_1_2\mathrm{}_n(𝐂^2)^{(n)}$$

One of the fundamental facts of the classical theory of computing is that there exist finite sets of simple functions, called *fundamental* or *universal gates*, such that any function $`f:𝐁^n𝐁^m`$, where $`𝐁`$ is the Boolean algebra $`\{0,1\}`$, can be constructed in a simple explicit way from them. For example the set {NOT, OR, AND} is such a basis of classical computation.

Similarly, one considers *quantum gates* on $`k`$ qubits (or $`k`$-gates) — unitary $`2^k\times 2^k`$-matrices acting on the quantum state space of $`k`$ qubits. A fundamental problem in quantum computing is to find basis of gates, which is “universal”. (More precisely, one distinguishes between universal bases and more restrictive *exactly universal* ones; a basis $`B`$ is called exactly universal if, for each $`k2`$, every unitary $`k`$-qubit operator can be obtained exactly by a circuit made up of the $`k`$-qubit gates produced from the elements of the basis $`B`$. See )

In real world computation, important rôle play devices independent from the environment noise. Thus two of the main requirements on error-free operations are to have a set of gates that is both universal for quantum computing and that can operate in noise-producing environment, i. e. is fault-tolerant (see ). For gates involving irrational multiples of $`\pi `$, called non-elementary, fault-tolerant realization is impossible. Thus, presence of this property makes an elementary gate inappropriate for physical realization.

It is feasible that any quantum system which one would consider potentially useful for quantum computing should contain a generating set of gates, which forms a basis in the above sense.

Two bases $`A`$ and $`B`$ are called equivalent, if the gates in the basis $`A`$ can be exactly realized using gates in the basis $`B`$ and vice versa.

Let us introduce the following operators:

$$\sigma _x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _z=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\sigma _y=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _z^\alpha =\left(\begin{array}{cc}1& 0\\ 0& e^{i\pi \alpha }\end{array}\right),$$

where $`\sigma _x,\sigma _y,\sigma _z`$ are Pauli matrices and have the property: every traceless and Hermitian $`2\times 2`$ unitary matrix $`U`$ can be expressed as $`U=x\sigma _x+y\sigma _y+z\sigma _z,`$ where $`x,y,z𝐑`$ and $`x^2+y^2+z^2=1`$.

For every unitary operator $`U\mathrm{U}(2)`$ we define the controlled $`U`$-operator as

$$\mathrm{\Lambda }_k(U)|x_1,\mathrm{},x_k|\xi =\{\begin{array}{cc}|x_1,\mathrm{},x_kU|\xi \hfill & \text{if }_{j=1}^kx_j=1\hfill \\ |x_1,\mathrm{},x_k|\xi \hfill & \text{if }_{j=1}^kx_j=0\hfill \end{array}$$

for all $`x_1,\mathrm{},x_k𝐁`$. Here $`_{j=1}^kx_j`$ denotes action of the operator AND on the boolean variables $`\{0,1\}`$. For example, if $`k=0`$, then $`\mathrm{\Lambda }_0(\sigma _x)=\mathrm{NOT}_q:HH`$, which is called the *quantum NOT* operator and by definition one has $`\sigma _x=\mathrm{NOT}_q`$, which acts on the state $`\alpha |0+\beta |1`$ via $`\alpha |0+\beta |1\beta |0+\alpha |1`$. If the quantum state $`\alpha |0+\beta |1`$ is written in vector form as $`(\alpha ,\beta )^{}`$, then NOT$`{}_{q}{}^{}(\alpha ,\beta )=(\beta ,\alpha )`$. Second important example is the so called *controlled-not* operator — cNOT: $``$, which acts on the pairs of qubits and carries out bitwise summation cNOT$`|u,v=|u,vu`$, where $``$ is addition mod 2. Let us consider one more operator — the *controlled controlled NOT* operator ccNOT: $``$. The classical one and two bit operators OR and AND can be expressed via ccNOT as follows:

$`\mathrm{ccNOT}|1,1,a`$ $`=|1,1,(\mathrm{NOT}a),`$
$`\mathrm{ccNOT}|a,b,0`$ $`=|a,b,(a\mathrm{AND}b),`$
for $`a,b𝐁`$.

On the other hand, every Boolean function can be expressed in the basis {NOT, OR, AND} and therefore there exists the unitary operator

$$U_f:\mathrm{}\mathrm{}$$

such that $`U_f|x,0=|x,f(x)`$.

Two and three qubit operations, respectively, cNOT and ccNOT, are represented by one qubit operations as follows:

$`\mathrm{cNOT}`$ $`=|01|\mathrm{𝟏}+|11|\mathrm{NOT},\mathrm{ccNOT}`$ $`=|00|\mathrm{𝟏}\mathrm{𝟏}+|11|\mathrm{cNOT}.`$
This means that if it is possible to find a finite subset of U(2) which generates a dense subset of U(2), then we obtain all operators for quantum computing.

The universal sets of gates for computation have been extensively studied. To emphasize importance of such works let us recall some fundamental facts which concern relations between different bases.

The set of gates $`\{\mathrm{\Lambda }_2(\sigma _x),\sigma _z^{1/2},H\}`$, where

$$H=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),$$

the so called Hadamard operator, is a fault-tolerant basis . On the other hand the set of gates $`\{\mathrm{\Lambda }_1(\sigma _z^{1/2}),H\}`$ is universal  and this basis is equivalent to the basis $`\{\mathrm{\Lambda }_2(\sigma _x),\sigma _z^{1/2},H\}`$. In  is considered the basis $`\{\mathrm{\Lambda }_1(\sigma _x),\sigma _z^{1/4},H\}`$ and proved that it is not equivalent to the basis $`\{\mathrm{\Lambda }_2(\sigma _x),\sigma _z^{1/2},H\}`$.

## Connection with singularity as generator of quantum gates

In  is considered the model of quantum computing, called by authors holonomic, based on the well known Berry phase. Main ingredient in holonomic quantum computation is a smooth vector bundle $`EX`$ with fibre $`𝐂^2\mathrm{}𝐂^2`$ and unitary connection $`\mathrm{\Omega }`$. The encoding space of information is in this case the fibre of the bundle and processing of information is represented by the holonomy operator $`P\mathrm{exp}(_\gamma \mathrm{\Omega })`$, which acts on the encoding space as $`vP\mathrm{exp}(_\gamma \mathrm{\Omega })v`$, where $`\gamma :[0,1]X`$ is a smooth path and $`P`$ denotes path-ordered exponential.

In the present section is given the construction of a universal set of gates from the monodromy representation of any system of differential equations of Fuchs type, which can be considered as dynamic equation of a quantum system. By our opinion, sets of gates obtained in this way are expected to be fault-tolerant by reasons similar to ones given in : action of a monodromy operator corresponding to a loop is unchanged by small fluctuations of the loop, provided they do not result in crossing any singularities.

At first let us review a simple example, the differential equation with regular singular points (for more details see ) $`\frac{df}{dz}=\frac{a}{z}f`$, $`a𝐂`$, on $`𝐂\{0\}`$. The solution of this equation is the many valued function $`f(z)=z^a`$ which by definition is $`z^a=e^{a\mathrm{log}z}`$, which under analytic continuation along a path $`\gamma `$ looping once counterclockwise around the origin transforms into $`e^{2\pi ia}z^a`$. The fundamental group of $`𝐂\{0\}`$ is isomorphic to $`𝐙`$, with generator $`\gamma `$, and the monodromy representation in GL$`{}_{1}{}^{}(𝐂)`$ is given by $`\gamma e^{2\pi ia}`$. Analogously, let $`A`$ be a Hermitian $`2\times 2`$ constant matrix, then the matrix function $`(zs)^A`$ is a many valued function on the domain $`|zs|>0`$ and is the solution of the matrix differential equation $`\frac{dF}{dz}=\frac{A}{z}F`$. The monodromy representation of this system is given by $`\gamma e^{2\pi iA}`$, and as $`A`$ is Hermitian, $`e^{2\pi iA}`$ is unitary and therefore we obtain the monodromy representation in the unitary group. In particular, if we can obtain any one-qubit gate, for example $`\tau `$, it will suffice to choose a Hermitian $`2\times 2`$ matrix $`\alpha `$, such that $`e^{2\pi i\alpha }=\tau `$, then $`e^{2\pi i\alpha }`$ acts as locally unitary operator (processing of information) on the fibre (encoding space of information) $`𝐂^2`$ of the trivial vector bundle $`𝐂^2\times 𝐂^{}𝐂^{}`$.

In  we consider the case, when the holomorphic vector bundle is given on a punctured compact Riemann surface. In this case, fundamental group of the base of the bundle is a free group and obtaining a universal set of logical gates does not present a difficult problem after application of solution methods of the Riemann-Hilbert problem.

Let $`\mathrm{𝐂𝐏}^n`$ be the $`n`$-dimensional complex projective space, and let $`s_1,s_2,s_3,s_4\mathrm{𝐂𝐏}^1`$ and $`s_j\mathrm{}`$. Denote by $`X_4=\mathrm{𝐂𝐏}^1\{s_1,s_2,s_3,s_4\}`$ and by $`\gamma _1,\gamma _2,\gamma _3,\gamma _4`$ the generators of $`\pi _1(X_4,z_0)`$ with relations $`\gamma _1\gamma _2\gamma _3\gamma _4=1`$. Let $`M_1,M_2,M_3,M_4`$ be matrices such that $`E_j=\frac{1}{2\pi i}\mathrm{ln}M_j,j=1,2,3,4`$ and $`M_1,M_2,M_3`$ generate a basis of the Lie algebra of SU(2) and satisfy the condition $`M_1M_2M_3M_4=1`$. Consider the representation $`\rho :\pi _1(X_4,z_0)\mathrm{SU}(2)`$ defined by the map $`\gamma _jM_j`$. Then for $`\rho `$ there exits a system of differential equations of Fuchs type

$$dF(z)=(\underset{j=1}{\overset{4}{}}\frac{A_j}{zs_j}dz)F(z),$$

whose monodromy representation coincides with $`\rho `$ . The monodromy representation $`\rho `$ induces two-dimensional vector bundle $`E\mathrm{𝐂𝐏}^1`$ with meromorphic connection form

$$\omega =\underset{j=1}{\overset{4}{}}\frac{A_j}{zs_j}dz$$

and Chern number $`c_1(E)=_{j=1}^4E_j`$. The solution space $``$ of the system is a two-dimensional vector space. Moreover $`\pi _1(X_4,z_0)`$ acts on $``$ and any unitary operator can be obtained in this way.

Therefore, we have proved the following proposition:

Proposition 1 . The connection form $`\omega `$ generates the basis for the computation.

Consider the general case. Let $`X`$ be a compact Riemann surface of genus $`g2`$ with marked point $`z_0`$. The fundamental group of $`X^{}=X\{z_0\}`$ is generated by $`\alpha _1,\beta _1,\mathrm{},\alpha _g,\beta _g,\gamma `$ with relations $`\mathrm{\Pi }_{j=1}^g[\alpha _j,\beta _j]=\gamma `$, where $`\alpha _1,\beta _1,\mathrm{},\alpha _g,\beta _g`$ are generators of the fundamental group $`\pi _1(X)`$ and $`\gamma `$ is a loop going around $`z_0`$. Consider such homomorphism $`\rho :\pi _1(X)`$ SU(2) that Im$`\rho `$ is a dense subgroup of SU(2). The homomorphism $`\rho `$ defines a two-dimensional holomorphic vector bundle $`E^{}X^{}`$ and there exists a SU(2)-system of differential equations $`Df=\omega f`$ on $`X`$ which has one regular singular point $`z_0`$ and the monodromy representation of this system coincides with $`\rho `$. The pair $`(E^{},\omega )`$ can be extended to a possibly holomorphically nontrivial bundle $`EX`$, for which $`\omega `$ can be a meromorphic connection. Therefore we proved the following proposition

Proposition 2 . The connection $`\omega `$ of the bundle $`EX`$ densely generates all unitary operators $``$.

Below we consider the case when the fundamental group of the base is not free.

Suppose $`D=_{j=1}^{m+1}D_j`$ is a divisor, where $`D_j`$, $`j=1,\mathrm{},m+1`$ are hyperplanes in $`\mathrm{𝐂𝐏}^n`$. The main result of this paper is the following theorem.

Theorem. There exists a Fuchs type Pfaff system

$$df=\omega f$$

on $`\mathrm{𝐂𝐏}^nD`$ whose monodromy representation gives a universal set of quantum gates.

Let us choose a line $`L`$ which intersects the divisor $`D=_{j=1}^{m+1}D_j`$ at nonsingular points $`a_j=LD_j`$, $`j=1,2,\mathrm{},m+1`$. Consider the fundamental group $`\pi _1(L\{a_1,\mathrm{},a_{m+1}\},z_0)`$. Let $`[\gamma _1],\mathrm{},[\gamma _m]`$ be the generators of $`\pi _1(L\{a_1,\mathrm{},a_{m+1}\},z_0)`$, where the loops have form $`\gamma _j=\sigma _j[\alpha _j]\sigma _j^1`$, $`\alpha `$ is a path from $`z_0`$ to a neighborhood $`V_{a_j}`$ and $`\sigma _j`$ is small a loop in the the neighborhood $`V_{a_j}`$ which generates $`\pi _1(V_{a_j}\{a_j\})𝐙`$.

It is known that $`[\gamma _1],\mathrm{},[\gamma _m]`$ too are generators of $`\pi _1(\mathrm{𝐂𝐏}^nD,z_0)`$ under some conditions.

Suppose we have the family of representations

(2)
$$\rho _\lambda :\pi _1(\mathrm{𝐂𝐏}^nD,z_0)\mathrm{GL}_m(𝐂)$$

such that

(3)
$$\rho _\lambda ([\gamma _j])=1+\lambda M_1^j+\lambda ^2M_2^j+\mathrm{}+\lambda ^kM_k^j+\mathrm{},$$

where $`M_k^j`$ are $`m\times m`$-matrices. The family of representations (2) satisfying the condition (3) is called analytic.

Suppose $`D_j=\{h_j=0\}`$ and consider the 1-form $`\mathrm{\Omega }_k=_{j=1}^{m+1}U_k^j\frac{dh_j}{h_j}`$, where $`U_k^j`$ are constant matrices and $`_{j=1}^{m+1}U_k^j=0`$. Consider the family of meromorphic 1-forms

(4)
$$\mathrm{\Omega }(\lambda )=\lambda \mathrm{\Omega }_1+\lambda ^2\mathrm{\Omega }_2+\mathrm{}+\lambda ^k\mathrm{\Omega }_k+\mathrm{}.$$

If $`U^j(\lambda )=_{k=1}^{\mathrm{}}\lambda ^kU_k^j`$ is a converging power series then the family $`\mathrm{\Omega }(\lambda )`$ is called analytic family of Fuchs systems.

It is known that for every $`\lambda `$ this system satisfies the condition

$$d_z\mathrm{\Omega }(\lambda )=0,zC.$$

If the analytic family of Fuchs systems satisfies the condition

$$\mathrm{\Omega }(\lambda )\mathrm{\Omega }(\lambda )=0,$$

then the family of Fuchs systems is called integrable.

Similar terminology will be used for a Pfaff system $`df=\mathrm{\Omega }(\lambda )f`$.

The monodromy representation of the integrable family of Pfaff systems

(5)
$$df=\mathrm{\Omega }(\lambda )f$$

is an analytic family of representations of the fundamental group $`\pi _1(\mathrm{𝐂𝐏}^nD,z_0)`$ .

Let $`\mathrm{\Phi }`$ be the fundamental matrix of the system (5) represented by the Peano series

(6)
$$\mathrm{\Phi }=1+\mathrm{\Omega }(\gamma )+\mathrm{\Omega }(\gamma )\mathrm{\Omega }(\gamma )+\mathrm{}+\mathrm{\Omega }(\gamma )\mathrm{}\mathrm{\Omega }(\gamma )+\mathrm{},$$

where $`\mathrm{\Omega }(\gamma )\mathrm{}\mathrm{\Omega }(\gamma )`$ is the Chen iterated integral. If such representation of the fundamental matrix is chosen, then for every $`\gamma \pi _1(\text{CP}^nD,z_0)`$ one will have

$$\rho ([\gamma ])=1+_\gamma \mathrm{\Omega }(\gamma )+_\gamma \mathrm{\Omega }(\gamma )\mathrm{\Omega }(\gamma )+\mathrm{}+_\gamma \mathrm{\Omega }(\gamma )\mathrm{}\mathrm{\Omega }(\gamma )+\mathrm{}.$$

Theorem 1 . For every analytic family of representations (2), when the parameter is small, there exists a family of Pfaff systems (5) whose monodromy coincides with (2).

Sketch of proof. We will show, that there exists a family of Fuchs type Pfaff systems (5), where $`\mathrm{\Omega }(\lambda )`$ has the form (4), $`\mathrm{\Omega }_k=_{j=1}^mU_k^j\omega _j`$, $`\omega _j=\frac{dh_j}{h_j}\frac{dh_{m+1}}{h_{m+1}}`$, so that its family of representations coincides with (2).

Step 1. By (2) we find $`U_k^j`$, $`k=1,2,\mathrm{}`$, $`j=1,2,\mathrm{},m`$. Indeed, let $`\rho _\lambda ([\gamma _j])`$ be represented as Peano series like (6), and rewrite it in Lappo-Danilevski form :

(7)
$$\rho ([\gamma ])=1+\underset{j=1}{\overset{m}{}}_\gamma \omega _jU^j(\lambda )+\underset{j,k=1}{\overset{m}{}}_\gamma \omega _j\omega _kU^j(\lambda )U^k(\lambda )+\mathrm{}+$$

$$+\underset{j_1,\mathrm{},j_k=1}{\overset{m}{}}_\gamma \omega _{j_1}\mathrm{}\omega _{j_k}U^{j_1}(\lambda )U^{j_k}(\lambda );$$

for the generators we have (2) and therefore we obtain:

$$\lambda M_1^j+\lambda ^2M_2^j+\mathrm{}+\lambda ^kM_k^j+\mathrm{}=$$

$$=\lambda _{\gamma _j}\mathrm{\Omega }_1+\lambda ^2(_{\gamma _j}\mathrm{\Omega }_1\mathrm{\Omega }_1+_{\gamma _j}\mathrm{\Omega }_2)+\mathrm{}$$

$$\mathrm{}+\lambda ^k(_{\gamma _j}\mathrm{\Omega }_k+\underset{q=2}{\overset{m}{}}\underset{k_1+\mathrm{}+k_q=k}{}_{\gamma _j}\mathrm{\Omega }_{k_1}\mathrm{}\mathrm{\Omega }_{k_q})+\mathrm{}$$

This implies

$$_{\gamma _j}\mathrm{\Omega }_1=M_1^j,j=1,\mathrm{},m$$

$$_{\gamma _j}\mathrm{\Omega }_1\mathrm{\Omega }_1+_{\gamma _j}\mathrm{\Omega }_2=M_2^j,j=1,\mathrm{},m,$$

$$_{\gamma _j}\mathrm{\Omega }_k+\underset{q=2}{\overset{m}{}}\underset{k_1+\mathrm{}+k_q=k}{}_{\gamma _j}\mathrm{\Omega }_{k_1}\mathrm{}\mathrm{\Omega }_{k_q}=M_k^j,j=1,\mathrm{},m$$

and so on. As $`_{\gamma _j}\mathrm{\Omega }_k=2\pi iU_k^j`$, for $`U_k^j`$, $`j=1,\mathrm{},m,k=1,2,\mathrm{}`$ we obtain

$$U_1^j=\frac{1}{2\pi i}M_1^j,$$

$$U_2^j=\frac{1}{2\pi i}(M_2^j_{\gamma _j}\mathrm{\Omega }_1\mathrm{\Omega }_1)=\frac{1}{2\pi i}(M_2^j\underset{k_1,k_2=1}{\overset{m}{}}_{\gamma _j}\omega _{k_1}\omega _{k_2}U_1^{k_1}U_1^{k_2}),$$

$$\mathrm{}\mathrm{}\mathrm{}.$$

$$U_k^j=\frac{1}{2\pi i}(M_k^j\underset{q=2}{\overset{k}{}}\underset{k_1+k_2+\mathrm{}+k_q=k}{}_{\gamma _j}\mathrm{\Omega }_{k_1}\mathrm{}\mathrm{\Omega }_{k_q}),$$

$`j=1,2,\mathrm{},m.`$ Therefore we obtain a formal family of 1-forms

(8)
$$\mathrm{\Omega }(\lambda )=\lambda \mathrm{\Omega }_1+\lambda ^2\mathrm{\Omega }_2+\mathrm{}+\lambda ^k\mathrm{\Omega }_k+\mathrm{}=\underset{j=1}{\overset{m}{}}U^j(\lambda )\omega _j,$$

$$\mathrm{\Omega }_k=\underset{j=1}{\overset{m}{}}U_k^j\omega _j,U^j(\lambda )=\underset{k=1}{\overset{\mathrm{}}{}}\lambda ^kU_k^j.$$

Step 2. The formal series (8) is convergent for small $`\lambda `$ .

Step 3. The family of 1-forms (8) is integrable, i. e. the identity $`\mathrm{\Omega }(\lambda )\mathrm{\Omega }(\lambda )=0`$ is satisfied .

Theorem 2. Let

$$\rho :\pi _1(\mathrm{𝐂𝐏}^nD,z_0)\mathrm{GL}_m(𝐂)$$

be a representation such that $`\rho ([\gamma _j])`$, $`j=1,\mathrm{},m`$, are close to identity. Then $`\rho `$ is realizable as monodromy representation of an integrable Fuchs system $`df=f\mathrm{\Omega }`$, $`\mathrm{\Omega }=_{j=1}^mU^j\omega _j`$, where $`U^j`$ are close to the zero matrix.

We will apply results of this subsection to a special divisor. In particular, suppose $`D=_{i<j}D_{ij}D_0`$, where $`D_{ij}=\{(z_0,\mathrm{},z_n)\mathrm{𝐂𝐏}^n|z_i=z_j,i,j0\}`$, $`D_0=\{z\mathrm{𝐂𝐏}^n|z_0=0\}`$. Let $`X_n=\mathrm{𝐂𝐏}^nD=\{(z_1,\mathrm{},z_n)𝐂^n|z_iz_j,ij\}`$. The fundamental group of $`X_n`$ is called the pure braid group on $`n`$ strings, which we will denote by $`P_n`$. The symmetric group $`S_n`$ acts on $`X_n`$ by $`g.(z_1,\mathrm{},z_n)=(z_{g(1)},\mathrm{},z_{g(n)})`$, $`gS_n`$. The fundamental group of the quotient space $`X_n/S_n`$ is called the braid group on $`n`$ strings and denoted by $`B_n`$. The braid group has $`n1`$ generators $`\sigma _1,\mathrm{},\sigma _{n1}`$ with relations

(9)
$$\sigma _i\sigma _{i+1}\sigma _i=\sigma _{i+1}\sigma _i\sigma _{i+1},1in2,$$

(10)
$$\sigma _i\sigma _j=\sigma _j\sigma _i,|ij|2.$$

Choose a point $`z_0=(1,2,\mathrm{},n)X_n`$ and denote by $`pr:X_nX_n/S_n`$ the natural projection. We have exact sequence of groups:

$$1P_nB_nS_n1$$

and generators of the pure braid group $`P_n`$ are $`\tau _{ij}`$, $`1i<jn`$, where

(11)
$$\tau _{ij}=\sigma _i\sigma _{i+1}\mathrm{}\sigma _{j1}\sigma _j^2\sigma _{j1}^1\mathrm{}\sigma _i^1,1i<jn$$

satisfying the relations

$$\tau _{rs}\tau _{ik}=\tau _{ik}\tau _{rs},\text{ for }s<i\text{ or }k<s,$$

$$\tau _{ks}\tau _{ik}\tau _{ks}^1=\tau _{is}^1\tau _{ik}\tau _{is}\text{ for }i<k<s,$$

$$\tau _{rk}\tau _{ik}\tau _{rk}^1=\tau _{ik}^1\tau _{ir}^1\tau _{ik}\tau _{ir}\tau _{ik}\text{ for }i<r<k,$$

$$\tau _{rs}\tau _{ik}\tau _{rs}^1=\tau _{is}^1\tau _{is}\tau _{ir}\tau _{ik}\tau _{ir}^1\tau _{is}^1\tau _{ir}\tau _{is}\text{ for }i<r<k<s.$$

Consider the matrix valued 1-form

(12)
$$\mathrm{\Omega }=\underset{1i<jn}{}\mathrm{\Omega }_{ij}dlog(z_iz_j)$$

holomorphic on $`X_n`$, where $`\mathrm{\Omega }_{ij}M_m(C),1i<jn`$. Let $`EX_n`$ be a holomorphic rank $`m`$ trivial vector bundle with connection $``$ for which $`\mathrm{\Omega }`$ is connection form. Holomorphic sections $`f=(f^1,\mathrm{},f^m)`$ of this bundle are solutions of the Fuchs system

(13)
$$df=\mathrm{\Omega }f,$$

where $`\mathrm{\Omega }`$ has the form (12).

Proposition 3. The system (13) is integrable if and only if the following condition is satisfied

(14)
$$[\mathrm{\Omega }_{ij},\mathrm{\Omega }_{ik}+\mathrm{\Omega }_{jk}]=[\mathrm{\Omega }_{ij}+\mathrm{\Omega }_{ik},\mathrm{\Omega }_{ji}],i<j<k,$$

(15)
$$[\mathrm{\Omega }_{ij},\mathrm{\Omega }_{kl}]=0,\text{ for distinct }i,j,k,l.$$

By theorem 2 we have :

Theorem 3. If $`\rho :P_n\mathrm{GL}_m(𝐂)`$ is such representation that $`\rho (\tau _{ij})1`$ is sufficiently small for each $`1i<jn`$, then there exist matrices $`\mathrm{\Omega }_{ij}M_m(𝐂)`$, $`1i<jn`$ close to $`0`$, which satisfy the conditions (14)-(15) and monodromy representation of the Fuchs system

$$df=\underset{1i<jn}{}\mathrm{\Omega }_{ij}\frac{d(z_iz_j)}{z_iz_j}f$$

coincides with $`\rho `$.

Let us apply the above results for designing logical gates for quantum computing. Choose such $`n`$ and representation $`\rho :P_n\mathrm{SU}(m)`$ that, the image Im$`\rho (P_n)`$ is dense in SU(m) and $`\rho (\tau _{ij})1`$ is sufficiently small for each $`1i<jn`$. Then the conditions of theorem 3 are satisfied and therefore there exists the unitary connection

(16)
$$=\underset{1i<jn}{}\mathrm{\Omega }_{ij}\frac{d(z_iz_j)}{z_iz_j}$$

on the trivial $`m`$-dimensional vector bundle $`𝐂^m\times X_nX_n`$ which gives the action on $`𝐂^m`$ as follows: $`𝐂^mv\rho (\tau )v`$, for every $`\tau P_n`$.

In this manner we obtain universal quantum logical gates which are based on holomorphic vector bundles over $`X_n`$ with integrable connection (16).

## Conformal field theory

Consider a 1-form $`\mathrm{\Omega }`$ of special type defined in (12). Let $`V_1,\mathrm{},V_m`$ be $`𝔰l_2`$-modules. Put $`V=V_1\mathrm{}V_m`$. Let $`\{I_j\}`$ be an orthogonal basis of $`𝔰𝔩_2`$ and $`𝐜=_jI_jI_jU(𝔰𝔩_2)`$ the Casimir element in the universal enveloping algebra $`U(𝔰𝔩_2)`$. Let $`\mathrm{\Delta }:U(𝔰𝔩_2)U(𝔰𝔩_2)U(𝔰𝔩_2)`$ be the diagonal homomorphism determined by $`\mathrm{\Delta }(x)=x1+1x`$, $`x𝔰𝔩_2`$. Set $`\mathrm{\Omega }=\frac{1}{2}(\mathrm{\Delta }𝐜𝐜11𝐜)`$. Consider a family $`\rho _i:𝔰𝔩_2\mathrm{End}(V_i)`$, $`i=1,\mathrm{},n`$ of irreducible representations and define the representations

(17)
$$\rho _i\rho _j:𝔰𝔩_2\mathrm{End}(V_1\mathrm{}V_n)$$

by formulæ

$$(\rho _i\rho _j)(x)=I_1\mathrm{}\rho _i(x)\mathrm{}I_j\mathrm{}I_n+I_1\mathrm{}I_i\mathrm{}\rho _j(x)\mathrm{}I_n,$$

where $`I_k`$ denotes the identity operator acting on $`V_k`$.

The representations (17) extend to the universal enveloping algebra $`U(𝔰𝔩_2)`$, these representations we denote again by $`\rho _i\rho _j`$; thus we have the representations $`\rho _i\rho _j:U(𝔰𝔩_2)\mathrm{End}(V_1\mathrm{}V_n)`$. Let $`\mathrm{\Omega }_{ij}=(\rho _i\rho _j)\mathrm{\Omega }`$. The linear operators $`\mathrm{\Omega }_{ij}:VV`$, $`i<j`$, act as $`\mathrm{\Omega }1+1\mathrm{\Omega }`$ on $`V_iV_j`$ and trivially on all of the other factors. The Fuchs type Pfaff system

(18)
$$\frac{\mathrm{\Psi }(z_1,\mathrm{},z_n)}{z_i}=\frac{1}{\lambda }\underset{j=1,ij}{\overset{n}{}}\frac{\mathrm{\Omega }_{ij}}{z_iz_j}\mathrm{\Psi },i=1,\mathrm{},n,$$

where $`\mathrm{\Psi }(z_1,\mathrm{},z_n)`$ is a $`V`$-valued function on $`X_n`$, is called the Knizhnik-Zamolodchikov equation. Here $`\lambda `$ is a complex parameter. Solutions of (18) are covariant constant sections of the trivial bundle $`X_n\times VX_n`$ with flat connection

$$\underset{\begin{array}{c}j=1\\ ij\end{array}}{\overset{n}{}}\frac{\mathrm{\Omega }_{ij}}{z_iz_j}d(z_iz_j).$$

It follows from previous subsections that (18) is integrable if and only if the conditions (14)-(15) are satisfied. The monodromy representation of (18) can be extended to a representation of the braid group . Let $`V_1=V_2=\mathrm{}=V_n\mathrm{C}^2`$, then we obtain that the braid group $`B_n`$ acts on $`V^n(𝐂^2)^n`$. For quantum computing an appropriate choice of this action is necessary (see ).

The Knizhnik-Zamolodchikov equation was invented in conformal field theory. Its solutions describe $`n+1`$-point correlation functions on the Riemann sphere in the Wess-Zumino-Witten model. This model is uniquely determined by choice of a simple Lie algebra $`𝔤`$ for every $`k>0`$, $`k𝐙`$. Any model of the conformal field theory has a certain set of primary fields $`\{\varphi _j\}`$. They are in one-to-one correspondence with irreducible representations $`\{\rho _j\}`$ of the Kac-Moody algebra $`\widehat{𝔤}`$. Let $`𝔤=𝔰𝔩_2(𝐂)`$. For each nonnegative half integer $`j`$ there exists a unique irreducible $`𝔤`$-module $`V_j`$, called that of spin $`j`$, with highest weight $`j\alpha `$ and dim$`V_j=2j+1`$.

Theorem 4 , . The $`n`$-point function $`\mathrm{\Psi }=\varphi _1(z_1)\mathrm{}\varphi _n(z_n)`$, where $`\varphi _j`$ are primary fields, satisfies the system of differential equations (18) and the monodromy representation of this system is unitarizable.

For example consider the case $`n=2`$ with

$$\mathrm{\Omega }=\mathrm{\Omega }_{12}=\mathrm{\Omega }_{21}=\left(\begin{array}{cccc}1/2& 0& 0& 0\\ 0& 1/2& 1& 0\\ 0& 1& 1/2& 0\\ 0& 0& 0& 1/2\end{array}\right),$$

then (18) have the form

(19) $`{\displaystyle \frac{F}{z_1}}={\displaystyle \frac{1}{\lambda }}{\displaystyle \frac{\mathrm{\Omega }_{12}}{z_1z_2}}F,`$
(20) $`{\displaystyle \frac{F}{z_2}}={\displaystyle \frac{1}{\lambda }}{\displaystyle \frac{\mathrm{\Omega }_{21}}{z_2z_1}}F,`$
where $`F`$ is a function on $`X_2=\text{C}^2\{z_1=z_2\}`$ with values in $`V^2`$, $`dim_\text{C}V=2`$. The solutions to the above system are given by $`F(z)=e^{\frac{1}{\lambda }\mathrm{ln}(z_1z_2)\mathrm{\Omega }}C`$, where $`C`$ is a constant 4-vector. The image of the generator $`\sigma _1B_2`$ by the monodromy representation of the system (19)-(20) is $`e^{\frac{\pi i}{\lambda }\mathrm{\Omega }}`$ , which can be considered as a nontrivial 2-qubit gate.

In  is given an example of such 4-point function, for which the corresponding monodromy representation is $`\sigma _x`$. In  is considered the exact solution and integrability of the reduced BCS model of superconductivity considered from the conformal field theory point of view.

Georgian Academy of Sciences

Institute of Cybernetics

e-mail: giorgadze@rmi.acnet.ge