Wavelet analysis: a new significance test for signals dominated by intrinsic red-noise variability
aa r X i v : . [ a s t r o - ph . H E ] J un SUBMITTED TO IEEE TRANSACTIONS ON SIGNAL PROCESSING 1
Wavelet analysis: a new significance test for signalsdominated by intrinsic red-noise variability
Paweł Lachowicz,
Member, IEEE
Abstract —We develop a new statistical test for the waveletpower spectrum. We design it with purpose of testing signalswhich intrinsic variability displays in a Fourier domain a red-noise component described by a single, broken or doubly-broken power-law model. We formulate our methodology asstraightforwardly applicable to astronomical X-ray light curvesand aimed at judging the significance level for detected quasi-periodic oscillations (QPOs). Our test is based on a comparisonof wavelet coefficients derived for the source signal with theseobtained from the averaged wavelet decomposition of simulatedsignal which preserves the same broad-band model of variabilityas displayed by X-ray source. We perform a test for statisticallysignificant QPO detection in XTE J1550–564 microquasar andactive galaxy of RE J1034+396 confirming these results in thewavelet domain with our method. In addition, we argue on theusefulness of our new algorithm for general class of signalsdisplaying /f α -type variability. Index Terms —wavelet transforms, signal analysis
I. I
NTRODUCTION
Fourier transform has become a very popular tool of time-series analysis within last fifty years in almost all areas ofresearch. Applied to test signal frequency content, it wasquickly recognized as a helpful method for investigationof periodic variability [1]. It gained a particular interestin a wide field of astronomical research regarding differentclasses of objects with variable emission. In particular, time-series provided by recent high-energy detectors onboard X-ray orbiting NASA/ESA satellites [2]–[3] have been foundto trace remarkably rapid variability coming from neutronstar and black-hole systems [4]–[6]. Since the majority of X-ray emission takes its origin in accretion process of matterand gas onto central object [7]–[8], therefore the observedsignal variability is a direct manifestation of physical processestaking place in a very strong gravitational field [9].The accretion onto black-holes is one of the most energeticprocess in the Universe [10]. Energy spectra for the corre-sponding objects uncover both thermal and non-thermal natureof emission (e.g. [11]). The analysis of Fourier power spectra(also referred to as Power Spectral Density; PSD) of X-raylight curves x ( t ) has revealed that the overall shapes are veryoften well described, in general, by a power-law model, i.e. P ( f ) ∝ /f α where α ∈ ℜ denotes a slope ([5], [12]–[16]). Incase of α > one claims about signals dominated by intrinsicvariability of red-noise type. The red-noise character of sourcevariability is generally rarely met in nature [17]. Nevertheless,it seems that black-hole systems emitting in X-rays enlargethe sample significantly. P. Lachowicz is with Centre for Wavelets, Approximation and Informa-tion Processing at Temasek Laboratories, National University of Singapore,5A Engineering Drive 1,
To date the majority of efforts devoted to our curiosity ofunderstanding the X-ray variability of astronomical objectswere based on the systematic study of PSD shapes and lookingfor occasionally observed quasi-periodic oscillations (QPOs;[5], [18]–[19]). Since Fourier PSD contains no informationabout time evolution of detected periodicities, the need ofuse of the time-frequency techniques emerged. The applicationof Short-Term Fourier Transform (STFT; e.g. [20]) provideda good solution in this domain and has been found helpfulin localizing QPOs (e.g. [21]–[22]). However, regardless thefrequency of dectected QPOs, STFT keeps the same time-frequency resolution as dictated by a length of a slidingwindow. In most cases it prevents a detection of oscillationslasting shorter than the window span. On the contrary, awavelet analysis occurred to be more successful where thetime-frequency resolution was scale (frequency) dependant.Within recent fifteen years it has attracted attention of manyresearchers [23]–[24].The wavelet power spectrum (scalogram) for discrete evenlysampled X-ray signal x n ( n = 1 , ..., N ; T = N ∆ t ) given ascount-rate [cts s − ] can be defined as the normalized squareof the modulus of the wavelet transform: W = ξ | w n ( a m ) | (1)where ξ denotes a normalization factor and w a discrete formof the continuous wavelet transform in a function of twoparameters: scale a and localized time index n [24]–[26]: w n ( a m ) = (cid:18) πa m ∆ t (cid:19) / N X j =1 ˆ x j ˆ ψ ⋆ (2 πa m f j ) e i πf j n ∆ t (2)where the discrete Fourier transform of signal x n is given by: ˆ x j = N X n =1 ( x n − ¯ x ) e − i (2 πjn/N ) (3)and j denotes a frequency index, ν j = ( j/ ( N ∆ t ) for j ≤ N/ − j/ ( N ∆ t ) for j > N/ (4)and ¯ x a mean value of x n [26]. Following [27] one canassume normalization factor of w n ( a m ) to be ξ = 2∆ t ¯ x − which provides that integrated wavelet spectrum over time,also referred to as a global wavelet power spectrum [26], G ( a m ) = ξN N X n =1 | w n ( a m ) | , (5)will hold the units of (rms/mean) Hz − , i.e. it gives thevariance, relative to the mean, within a given frequency range SUBMITTED TO IEEE TRANSACTIONS ON SIGNAL PROCESSING of integration [27]–[28]. For the computation of (1) the use ofa dyadic grid of scales, i.e. a m = a m ∆ m , m = 0 , ..., M, (6)can be implemented [26], where the smallest scale a = 2∆ t corresponds to the reversed Nyquist frequency, and M = ∆ m − log ( N ∆ ta − ) . (7)A grid resolution is given by ∆ m . Since the search forperiodicity and quasi-periodicity in X-ray light curves seemsto be the main target, it is reasonable to assume the analyzingfunction to be Morlet wavelet, ψ ( t ) = π − / e i (2 πf t ) e − t / , (8)which oscillates due to a term ∝ e i t , where πf parameteris set to 6.In order to popularize wavelet analysis as a new tool inresearch over time-series, an additional requirement in theform of significance test was strongly desired. Torrence &Compo [26] (TC98) were first who addressed that issueextensively providing the public with an excellent guide andsoftware ready-to-use. In general, a map of wavelet powerspectrum shows a distribution of spectral components in thetime–scale (or time–frequency) plane. Single oscillation, ofsufficiently high amplitude to be detected, marks itself in themap as a wavelet peak. To every peak a certain statisticalsignificance can be assessed. We say that a wavelet peak issignificant at assumed per cent of confidence when it is abovea certain background spectrum. The latter can be defined bythe mean Fourier power spectrum of analyzed time-series.TC98 in their work adopted the background power spec-trum, P j , after [29] (see also Eq.(16) in [26]) which is validfor application as long as the analyzed signal is consideredto be a realization of univariate lag-1 autoregressive (AR1)process [30]–[32], x l = αx l − + z l , where α is lag-1autocorrelation, z l denotes a random variable drawn from theGaussian distribution (with zero mean and variance σ ) and x = 0 . However, the PSD shape provided by [29] is not agood representation of broad-band PSDs in case of consideredhere X-ray sources. This is because of at least two reasons:(i) not every X-ray light curve can be described by AR1process [33] (although [34]), (ii) PSDs are dominated by thered-noise spectral components which shapes are mainly fittedwith single, broken, doubly broken or even bending power-lawmodels.TC98 showed that for a time-series modeled as the AR1process its local wavelet spectrum follows the mean Fourierspectrum given by their Eq.(16) [26]. By the “local waveletspectrum” they define a vertical slice (along the frequencyaxis) through a wavelet map. Empirical justification is doneon the way of Monte Carlo simulations which were performedconsidering only one local wavelet power spectrum per sim-ulation and, in addition, selected for a particular moment oftime ( k = 256 of 512 points). The corresponding distributionof the local wavelet spectrum (at each localized time index n and scale a m ), ought to have χ distribution, | w k ( a m ) | σ ⇒ P j χ , (9) where χ denotes the χ distribution with two degrees offreedom and an arrow “ ⇒ ” means “is distributed as”. Theabove implication is performed based on the assumption thata random variable x l is normally distributed.A noteworthy work in this domain has been recently doneby Zhongfu Ge [35] who generalized TC98’s significance test.Ge left the assumption of signal to be the realization of AR1process and proposed Gaussian White Noise process to serveas the null hypothesis.What TC98 and Ge did not show is that the local waveletpower spectra, P j , at nearby times n are correlated due tovariable length of the wavelet function, and the correlation isstronger at low frequencies as shown by [36]. This introducesdepartures off χ distribution in any finite data set. Thus amore correct procedure is required to take this effect intoaccount.Maraun, Kurths & Holschneider [37] proposed another testfor wavelet significance, namely, areawise test which aimsat comparing the wavelet peak size to the expected time-frequency uncertainty. They noticed that for any process withunknown distribution a number of wavelet peaks forming inthe wavelet map different shapes (circles, ovals, patches, etc.)become difficult in interpretation. The extracted picture ismessy. Some peaks are due to red-noise process of the sourceintrinsic variability and the others are due to white-noise. Asshown by them, the areawise test may reduce a number ofsignificant peaks even by 90 per cent.To judge the reality of each wavelet peak it seems bestto set the significance level highly enough in order to rejectspurious results. However, not always this is a right way, forinstance, in case of X-ray sources of rather low count-rate(low S/N ratio) and strong mixing between red- and white-noise components in high-frequency region. On the other hand,this also does not mean that the analysis does not containany scientific results. For example, the same situation tookplace for COBE detections of microwave background peaks[38]. Most of them were spurious but the fluctuations werestatistically different from expectations of pure noise [39].None of this method can be applied as a fully trustfultest working well for astronomical X-ray signals as long astheir assumptions rely on the analyses that use improprietebackground spectra. In this paper we address that issue andpropose a new significance test.II. W AVELET SIGNIFICANCE FOR RED - NOISE SIGNALSOF POWER - LAW FORM
We design a method for determination of significance ofwavelet spectral features (in short, wavelet peaks) calculatedfor any X-ray light curve which intrinsic variability is a real-ization of a sum of red- and white-noise (Poisson) processes.The main idea standing for that is to compare the distributionof wavelet power between real data and simulated signal. Here,by the simulated light curve (also referred to as TK light curve)we will understand a time-series which displays the sameprofile of red-noise variability as that one found for X-raydata. Before providing the reader with ready-to-use algorithm(Sect. II-C), some essential information on the background
ACHOWICZ: A NEW SIGNIFICANCE TEST FOR WAVELET ANALYSIS 3 assumptions need to be clarified. We provide them within thefollowing two subsections (Sect. II-A and II-B).
A. Simulations of /f α process Let’s assume that x ( t ) describes a part of stochastic con-tinuous process X ( t ) generated by the physical system ofunknown nature which produces the variability. Based solelyon x ( t ) a reproduction of the overall information and prop-erties of the system is very difficult. Recorded data canfacilitate a determination of a model which would reproducethe system best. The observed light curve of X-ray emission isthe realization of the physical process(es) for which calculatedFourier power spectrum (PSD) reveals the variability of power-law /f α -form (see references in Section I). Therefore, weare able to approximate one of the system’s parameters (i.e.its variability) by a power-law model which, as we believe,constitutes the best PSD description of the underlying process.Aperiodic broad-band variability as displayed in X-ray lightcurves can be thought, in the first approximation, as the real-ization of a linear stochastic process, weakly non-stationary.Timmer & Konig [40] proposed an algorithm (hereafter alsoreferred to as TK) for generating linear aperiodic signals, l ( t ) ,that exhibit P ( f ) ∝ /f α power-law spectrum, l ( t ) ∝ X f p P ( f ) cos [2 πf − φ ( f )] , (10)where the phase is randomized φ ( f ) ∈ [0 , π ] . Since theamplitudes are taken as p P ( f ) , [40] noticed that in orderto create power-law signal it is essential to randomize also theamplitude.A light curve obtained in this manner displays /f α -shape, however as argued by [33], it does not preservefundamental statistical properties of X-ray time-series foraccreting black-hole systems. The latter have been found tobe rather non-linear showing rms–flux relation at all time-scales [33], [41] (rms is defined as a square root of variancecalculated for a signal). According to [33], formally non-lineartype of variability for any linear input light curve can beobtained by taking simple exponential transformation of (10),i.e. x exp ( t ) = exp[ l ( t )] . To make the picture complete everysimulated light curve should include white noise. Following[25] it can be expressed as: x ( t ) = POISDEV [ x exp ( t )∆ t ] / ∆ t (11)where x ( t ) denotes simulated signal with Poisson noise, ∆ t is a light curve resolution (bin time) and POISDEV representsa subroutine generating a random number with Poisson distri-bution [42].
B. Selection of best /f α model for PSD of X-ray source A description of Fourier PSD of X-ray object under inves-tigation can be obtained in the process of numerical fittingof power-law model. As reviewed by [5], a small but sub-stantial subset of X-ray black-hole (and neutron star) systemsapart from exhibiting /f α -like broad-band variability, alsoshow quasi-periodic oscillations (QPOs). Hardly ever metas tiny spikes, they are much often resolved as broadened features characterized by the quality factor of Q = f / ∆ f where ∆ f denotes the peak’s full width at half maximum(FWHM) and f QPO central frequency. The broader peakthe higher ∆ f thus the lower Q . So far, the majority ofQPO peaks have been found to be surprisingly well describedby Lorentzians, e.g. [12], [16]. Lorentzians are defined as L ( f ) ∝ ∆ f / [ f − f ) + (∆ f / ] with a centroid frequency f and, if detected, they constitute an intergal component ofthe overall PSD model. In practice, fitting of power-law modelas a background model may or may not include fitting of L ( f ) component simultaneously. In the latter case, PSD data pointsdescribing Lorentzian can be simply omitted and these twoaforementioned approaches rely on the individual preferencesand fitting experience.In this work our wish will be to generate simulated X-raylight curve which will own the same power-law backgroundmodel of variability as fitted to the object’s PSD. This steprequires that the overall rms and the average count rate forsimulated time-series must match the corresponding samevalues for X-ray signal. For Fourier PSD normalized to unitof (rms/mean) Hz − [27], the root-mean-square (rms) can beobtained by integration of underlying PSD in a given band offrequencies [28] as follows:rms = σ = Z f f P ( f ) df ≈ j X j P j ∆ f (12)where /N ∆ t ≤ f < f ≤ / t . We argue that rms values,both obtained for simulated signal and X-ray light curve oughtto be calculated in the same ( f , f ) range of frequencies whatcan guarantee the best match. The integration (12) shouldomit PSD data points corresponding to QPO component(s),if present. C. A new algorithm for wavelet significance
We define a new algorithm leading to determination ofsignificance matrix for derived wavelet power spectrum by thefollowing steps: • Two-dimensional matrix of wavelet power spectrum (1)for a time-series x n is a set of ( M + 1) × N waveletcoefficients where M corresponds to (7) and N denotesa number of data points. Let us denote the results ofcomputation of (1) for X-ray source as: W src( M +1) × N . (13)Please notice that in practice instead of time–scale planeone often refers to wavelet spectrum displayed in time–frequency plane. If required, for Morlet wavelet, a trans-formation from scales to Fourier frequencies can beapplied according to f j = (1 . a m ) − relation [26]. • For a given discrete signal x n calculate its Fourier powerspectrum and find the best model which describes it(Section II-B). If the model for considered data set isknown in advance (e.g. from the existing literature), usethe best fitted parameters as the reference ones. • Generate one long simulated light curve (with the same ∆ t ), best at least two or three orders of magnitude longer SUBMITTED TO IEEE TRANSACTIONS ON SIGNAL PROCESSING than your signal under analysis. Use Timmer & K¨onig’smethod which allows to obtain simulated signals fromFourier power spectra of assumed shape (Section II-A).Take care to keep the rms variability at the same level asin the data (Section II-B). • From the simulated TK light curve, select randomly K signals of the same time duration as in original data set,i.e. N ∆ t , and compute for each of them the waveletpower spectrum according to (1). Shall us denote theresulting matrices as: W i ( M +1) × N where i = 1 , ..., K. (14)In calculation of (14) for normalization factor ξ (seeSect. I) assume every time a mean count rate ¯ x of thelocal TK signal i . • Find the maximum value between all W i matrices: ∀ W i ∃ x where x = max( W i ) . (15) • For every scale m of every W i calculate the histogram ofwavelet power distribution from 0 to x with a histogramresolution defined as ∆ h = x/ ( M + 1) . In consequence,the following matrices shall be created: H i ( M +1) × ( x/ ∆ h ) (16)consisting of M + 1 histograms. • Calculate a matrix containing averaged histograms forevery scale m : ¯ H ( M +1) × ( x/ ∆ h ) = 1 K K X i =1 H i . (17)This step is introduced in order to smooth out differentdistribution of wavelet power appearing in every H i foreach scale. Please note the higher K the more smoothedaveraged wavelet power histograms one obtains. • For every of M + 1 histograms stored in ¯ H , determinethe p -th quantile x p . Namely, let ¯ H ( m ) be a continuousrandom variable, then for < p < , the p -th quantile isthe x p such that Pr ( ¯ H ( m ) ≤ x p ) = p [44]. For instance,in deriving of the significance of wavelet peaks at the 95per cent confidence level one aims in determination of0.95-th quantile x p in a function of scale m (frequency).In practice, a numerical integration of normalized ¯ H ( m ) must be performed until an integral value exceeds 0.95.This determines quantile x p with an error of ∆ h for eachhistogram: X p = x p x p ...x mp ...x Mp ( M +1) × . (18) • Expand matrix X p to have dimensions ( M + 1) × N bycalculating a linear algebraic product of the matrices X p A term normalized means that an integral of ¯ H ( m ) must equal 1. and I : W bkg( M +1) × N obs ( i, j ) = X k =1 X p ( i, k ) I ( k, j ) , (19)where I is a matrix of ones of × N dimension. Denotethe result as a background wavelet matrix W bkg . • Finally, a matrix of (1 − p )100% significance, S , willidentify significant wavelet peaks appearing in computedwavelet power spectrum W src when the following con-dition will be met: S ( M +1) × N obs = W src W bkg > . (20)In proposed procedure a distribution of wavelet power ata given wavelet scale (frequency) is directly compared to thedistribution of power obtained for simulated light curves. Sincethe Fourier PSD model and rms variability in simulated time-series are assumed to match the same quantities as for theX-ray source data, they serve as a good reference level totest any excess of wavelet power above it. In addition, bytesting a wavelet power scale by scale, one takes into accountan influence of correlation of wavelet power over a range offrequency and time due to varying length of a wavelet function(Sect. I, [36]).III. A PPLICATION TO X- RAY OBSERVATIONS
In this section we provide the reader with two examples ofapplication of our new algorithm to real observations of black-hole systems emitting in X-rays: microquasar XTE J1550–564 [45] and active galactic nuclei of RE J1034+396. Incase of the former object, the source is known because ofdisplaying a strong X-ray aperiodic variability with a variety oflow-frequency (0.08–18 Hz) and high-frequency QPOs (100–285 Hz) [46]–[48]. Here, we examine its activity peaked at thefrequency of 4 Hz and derive corresponding QPO significancein the wavelet domain. As for the latter source, we select theactive galaxy of RE J1034+396 which has been recently foundto uncover the very first ever significant signature (over σ )of quasi-periodicity among all active galaxies [49]. Inspiredby Fourier results of [49], we verify their findings by determi-nation of wavelet significance for detected ≃ . × − Hzvariability.
A. XTE J1550–564
In order to examine XTE J1550–564 light curve variabilitywhich has been dominated by intrinsic red-noise component in − – Hz frequency band, after [50] we select
RXTE /PCAobservation of 40099-01-24-01 (1998-09-29) from the publicHEASARC archive. We reduce the data with LHEASOFTpackage ver. 6.6.1 applying the standard PCA selection criteriafor above data set (see http://heasarc.gsfc.nasa.gov for details).As the final product we extract a 2–13 keV light curve with abin size of ∆ t = 2 − s for the first part of PCA observation(i.e. N = 80065 , T = N ∆ t = 2502 s).Fig. 1 presents Fourier power spectrum density calculatedfrom 78 averaged spectra based on 1024 point data segments.The overall PSD shape can be split into two main compo-nents: power-law underlying spectrum and broad line profile ACHOWICZ: A NEW SIGNIFICANCE TEST FOR WAVELET ANALYSIS 5 − . P o w e r [(r m s / m ea n ) / H z ] Frequency (Hz)
Fig. 1. Power spectrum density of XTE J1550–564 as observed by
RXTE /PCA on 1998-09-29. Black solid line denotes fitted model composed ofdouble broken power-law model (red solid line) and Lorentzian (dotted line)representing red-noise background spectrum and quasi-periodic oscillation(QPO), respectively. Power spectrum for simulated TK time-series has beenmarked in blue (circle markers). representing red-noise source variability and QPO component,respectively. In the χ fitting process we found that a modelcomposed of three-segment broken power-law (i.e. with twobreak frequencies) and Lorentzian profile line [51] fit the PSDbest yielding reduced- χ of 10.9 at 59 degrees of freedom. Thecorresponding power-law slope values and break frequencieshave been determined as follows: α = 0 . ± . , α =0 . ± . , α = 3 . ± . , f = 1 . ± . Hz and f = 8 . ± . , whereas for Lorentzian the central frequencyof f = 4 . ± . Hz and its FWHM ∆ f = 0 . ± . Hzwere fitted providing QPO quality factor of Q ≃ . . Wemarked both model components in Fig. 1 by red solid lineand dotted line, respectively, where the overall best model hasbeen indicated by black solid line.Based on the above power-law model we generated a point long simulated signal (64 segments of 1024 points) withTK method. In simulation we assumed the sampling time of ∆ t = 2 − s and ensured that resulting rms value (here, derivedin 0.5–2 Hz and 8–12 Hz frequency band; see Sect. II-B fordetails) was nearly the same as compared to corresponding rmsvalues for X-ray source. We calculated the averaged FourierPSD for TK time-series and plotted it in Fig. 1 by blue circlemarkers. One can clearly notice that TK PSD reproduces wellthe assumed double broken power-law model component.Following the recipe given in Section II-C we aim atdetermination of wavelet confidence at 0.9995 level. From theentire XTE J1550–564 light curve, we select a 3072 point long test signal (96 s) for which we compute wavelet significancematrix (Eq. 20). Upper plot of Fig. 2 presents a dependenceof 0.9995-th wavelet quantile (Eq. 18) in a function of Fourierfrequency in 1–10 Hz band, both for X-ray data segment(black solid line) and TK simulation (red dashed line). It isstraightforward to note that at the significance of × − the quasi-periodic variability around f ≃ Hz is high abovethe mean (simulated) red-noise level. Finally, the bottom mapof Fig. 2 presents the wavelet power spectrum (Eq. 1) for −2 −1 Frequency [Hz] x p ( . − t h quan t il e ) Fig. 2. (upper)
The dependence of 0.9995-th quantiles in a functionof frequency calculated for selected 96 s data signal of XTE J1550–564(solid black line) and obtained from wavelet analysis of simulated TK series(red dashed line). See Section III-A for details. (bottom)
Correspondingwavelet power spectrum displayed with derived significance contours (0.9995confidence level). For better clarity of map reading we plot 96 s long waveletspectrum between 37–94 s only. considered data segment with overlaid contours of derivedsignificance.
B. RE J1034+396
A Seyfert 1 galaxy of RE J1034+396 has been foundand confirmed as the first active galaxy showing highlysignificant quasi-periodic oscillation (QPO) at the frequency f ≃ . × − Hz [49]. We reproduce their PSD in Fig. 3 ascomputed for signal segment of
XMM-Newton
EPIC observa-tion of RE J1034+396 conducted on 2007-05-31. Solid blackline crossing PSD denotes best-fitted power-law spectrumdescribed as: lg P ( f ) = − .
35 lg f − . (21)(M. Gierlinski, private communication). Upper two dashedlines mark computed confidence levels for detected period-icity at 99.73% ( σ ) and 99.99% confidence levels. In orderto verify aforementioned significance of QPO detection, wedownloaded from XMM-Newton
Science Operations Centrethe same EPIC pn+MOS1/MOS2 observation data set andreduced it in the same manner as described in [49]. We limitoriginal X-ray light curve (0.3–10 keV) to 601 point long time-series (corresponding to segment 2 of [49]) where the samplingtime of ∆ t = 100 s has been used. SUBMITTED TO IEEE TRANSACTIONS ON SIGNAL PROCESSING . P o w e r [ (r m s / m ean ) H z − ] P ( f ) Poisson noise −5 −4 −3 f (Hz) Fig. 3. Fourier power spectral density of an active galaxy of RE J1034+396[49] displaying a significant quasi-periodicity at frequency f ≃ . × − Hz ( > . of confidence; upper dotted line). Solid line denotes best power-law model fitted to the data P ( f ) ∝ f − . , whereas bottom horizontaldotted lines marks expected white-noise level. Reproduced with kind permis-sion of Marek Gierlinski. We compute wavelet power spectrum (1) in − – − Hzband and make use of our new algorithm in order to determinetwo contour plots marking 99.73% and 99.99% confidencefor wavelet results. For the purpose of wavelet backgroundsimulation, we generate point TK signal based on assumedsingle power-law model (Eq. 21). Since the underlying PSDmodel is known in advance, we perform an integration of itbetween − Hz and − Hz to find the correspondingrms value in this frequency range (rms ≃ . %). Within thesimulation process of TK time-series we check to make surethat in the same frequency band the rms of TK signal staysin agreement with derived value. Based on that, we determinetwo significance matrices (Eq. 20) and plot them together withwavelet power spectrum. The final form of all calculationspresents Fig. 4 where an outter and inner contour denote99.73% and 99.99% confidence level, respectively.Remarkably, a detected quasi-periodic oscillation at f ≃ . × − Hz appears to be very coherent in time thoughshowing small but noticeable period drift. We checked that theentire periodicity remains significant at the confidence level of99% throughout the signal segment duration. In contrast, asmarked in Fig. 4, only a small fraction of QPO activity exceedstop level of . significance whereas about 70% of detectedQPO signal stays significant at σ level.IV. C ONCLUSIONS
In this paper we reviewed a recent status of knowledgedevoted to wavelet significance tests and proposed a newmethod designed especially for signals which intrinsic vari-ability displays in Fourier power spectral density a power-lawform. Since a large number of astronomical X-ray sources hasbeen found to be characterized by the similar PSD shapes, we
Fig. 4. RE J1034+396: Wavelet power spectrum calculated for the seconddata segment of
XMM-Newton
EPIC/pn+MOS observation, as defined by [49].Quasi-periodic modulation is clearly detected at f ≃ . × − Hz frequencyand last throughout a whole light curve duration. Outer and inner contourscorrespond to derived 99.73% ( σ ) and 99.99% confidence level, respectively,with the use of new algorithm as provided in Section II-C. focused our attention to that class of objects and their time-series analysis.We presented a flexible algorithm helpful in determinationof the significance of quasi-periodic oscillations in the waveletdomain. Our method has been based on a comparison ofwavelet coefficients calculated for X-ray time-series with thosecomputed for simulated signal, where the latter displays thesame PSD broad-band model of variability as X-ray object.We argue that the usefulness of our new algorithm can beeasily extended to non-astronomical time-series which PSDcan be described by power-law model as well. If required,the omission of exponential transformation of linear signalgenerated with the use of TK method and inclusion of Poissonnoise should be avoided (Section II-A).A CKNOWLEDGMENT
The author would like to thank Prof. Didier Barret forinspiration which led to the accomplishment of this paper andDr. Marek Gierli´nski for a kind providing with an originalversion of Fig. 3. R
EFERENCES[1] R. N. Bracewell, ”The Fourier Transform and Its Applications”, McGraw-Hill, 1965, 2nd ed. 1978, revised 1986[2] J. H. Swank, “The XTE Mission”,
Bulletin of the American AstronomicalSociety , vol. 26, pp. 1420, 1994[3] C. Gabriel, M. Guainazzi and L. Metcalfe, “XMM-Newton: Passing FiveYears of Successful Science Operations”, Astronomical Data AnalysisSoftware and Systems XIV, vol. 347, pp. 425, 2005[4] P. Uttley, I. M. McHardy and I. E. Papadakis, “Measuring the broad-bandpower spectra of active galactic nuclei with RXTE”,
Monthly Notices ofthe Royal Astronomical Society , vol. 332, issue 1, pp. 231-250, 2002[5] W. H. G. Lewin and M. van der Klis, “Compact stellar X-ray sources”,Cambridge Astrophysics Series, No. 39. Cambridge, UK: CambridgeUniversity Press, 2006[6] R. Sunyaev, and M. Revnivtsev, “Fourier power spectra at high frequen-cies: a way to distinguish a neutron star from a black hole”,
Astronomyand Astrophysics , vol. 358, pp. 617-623, 2000[7] N. I. Shakura and R.A. Sunyaev, “Black holes in binary systems.Observational appearance”,
Astronomy and Astrophysics , vol. 24, pp. 337,1973[8] C. M. Urry and P. Padovani, “Unified Schemes for Radio-Loud ActiveGalactic Nuclei”,
Publications of the Astronomical Society of the Pacific ,vol. 107, p. 803
ACHOWICZ: A NEW SIGNIFICANCE TEST FOR WAVELET ANALYSIS 7 [9] M. van der Klis, “Neutron Star QPOs as Probes of Strong Gravity andDense Matter”, X-ray Timing 2003: Rossie and Beyond. AIP ConferenceProceedings, vol. 714, pp. 371-378, 2004[10] J. Frank, A. King A. and D. J. Raine, “Accretion Power in Astrophysics”,Cambridge University Press, 2002[11] C. Done and M. Gierli´nski, “Observing the effects of the event horizonin black holes”,
Monthly Notice of the Royal Astronomical Society , vol.342, issue 4, pp. 1041-1055, 2003[12] K. Pottschmidt, J. Wilms, M. A. Nowak, G. G. Pooley, T. Gleissner,W. A. Heindl, D. M. Smith, R. Remillard and R. Staubert, “Long termvariability of Cygnus X-1. I. X-ray spectral-temporal correlations in thehard state”,
Astronomy and Astrophysics , vol. 407, pp. 1039-1058, 2003[13] A. Lawrence, M. G. Watson, K. A. Pounds,and M. Elvis, “Low-frequency divergent X-ray variability in the Seyfert galaxy NGC4051”,
Nature , vol. 325, p. 694, 1987[14] I. M. McHardy, I. E. Papadakis, P. Uttley, M. J. Page and K. O. Mason,“Combined long and short time-scale X-ray variability of NGC 4051 withRXTE and XMM-Newton”,
Monthly Notices of the Royal AstronomicalSociety , vol. 348, issue 3, pp. 783-801, 2004[15] T. Belloni and G. Hasinger, “Variability in the noise properties of CygnusX-1”,
Astronomy and Astrophysics , vol. 227, no. 2, L33-L36, 1990[16] M. A. Nowak, B. A. Vaughan, J. Wilms, J. B. Dove and M. C. Begelman,“Rossi X-Ray Timing Explorer Observation of Cygnus X-1. II. TimingAnalysis”,
The Astrophysical Journal , vol. 510, issue 2, pp. 874-891,1999[17] W. Press, ”Flicker noises in astronomy and elsewhere”, Comments onAstrophysics, vol. 7, p. 103, 1978[18] M. van der Klis, “The QPO phenomenon”,
Astronomische Nachrichten ,vol. 326, issue 9, p.798-803, 2005[19] M. van der Klis, “Comparing Black Hole and Neutron Star Variability”,
Astrophysics and Space Science , vol. 300, issue 1-3, pp. 149-157, 2005[20] L. Cohen, “Time-frequency analysis”, Englewood Cliffs, NJ: Prentice-Hall, 1995[21] J. Wilms, M. A. Nowak, K. Pottschmidt, W. A. Heindl, J. B. Dove andM. C. Begelman, “Discovery of recurring soft-to-hard state transitions inLMC X-3”,
Monthly Notices of the Royal Astronomical Society , vol. 320,issue 3, pp. 327-340, 2001[22] D. Barret, W. Klu´zniak, J. F. Olive, S. Paltani and G. K. Skinner, “Onthe high coherence of kHz quasi-periodic oscillations”,
Monthly Noticesof the Royal Astronomical Society , vol. 357, issue 4, pp. 1288-1294, 2005[23] Y. Meyer and S. Roques, “Progress in wavelet analysis and appli-cations”, 1993, Proceedings of the International Conference ”Waveletsand Applications”, Toulouse, France, June 1992, Gif-sur-Yvette: EditionsFrontieres, ed. by Meyer, Yves; Roques, Sylvie, 1993[24] P. S. Addison, “Illustrated Wavelet Transform Handbook”, IoP, 2002[25] P. Lachowicz and B. Czerny, “Wavelet analysis of millisecond variabilityof Cygnus X-1 during its failed state transition”,
Monthly Notices of theRoyal Astronomical Society , vol. 361, issue 2, pp. 645-658, 2005[26] C. Torrence and G. P. Compo, “A Practical Guide to Wavelet Analysis”,
Bulletin of the American Meteorological Society , vol. 79, issue 1, pp.61-78, 1998[27] S. Vaughan, R. Edelson, R. S. Warwick and P. Uttley, “On characterizingthe variability properties of X-ray light curves from active galaxies”,
Monthly Notices of the Royal Astronomical Society , vol. 345, issue 4,pp. 1271-1284, 2003[28] M. van der Klis, “Quasi-periodic oscillations and noise in low-mass X-ray binaries”,
Annual review of astronomy and astrophysics . Volume 27.Palo Alto, CA, Annual Reviews, Inc., pp. 517-553, 1989[29] D. L. Gilman, F. J. Fuglister and J. M. Mitchell Jr., “On the PowerSpectrum of Red Noise”,
Journal of Atmospheric Sciences , vol. 20, issue2, pp. 182-184, 1963[30] S. L. Marple, Jr., “Digital Spectral Analysis with Applications”, PrenticeHall, Englewood Cliffs, 1987, Chapter 8[31] C. Chatfield, “The Analysis of Time Series. An Introduction.”, Chapman& Hall, A CRC Press Company, Sixth Ed., 2004[32] S. M. Kay, “Modern Spectral Estimation: Theory and Application”.Englewood Cliffs, NJ: Prentice-Hall, 1988[33] P. Uttley, I. M. McHardy and S. Vaughan, “Non-linear X-ray variabilityin X-ray binaries and active galaxies”,
Monthly Notices of the RoyalAstronomical Society , vol. 359, p. 345, 2005[34] K. Pottschmidt, M. Koenig, J. Wilms and R. Staubert, “Analyzing short-term X-ray variability of Cygnus X-1 with Linear State Space Models”,
Astronomy and Astrophysics , vol. 334, p. 201, 1998[35] Z. Ge, “Significance test for wavelet power and the wavelet powerspectrum”,
Ann. Geophys. , vol. 25, p. 2259, 2007[36] D. Maraun and J. Kurths, “Cross Wavelet Analysis. Significance Testingand Pitfalls”,
Nonlin. Proc. Geoph. , vol. 11(4), pp. 505-514, 2004 [37] D. Maraun, J. Kurths and M. Holschneider, “Nonstationary Gaussianprocesses in wavelet domain: Synthesis, estimation, and significancetesting”,
Physical Review E , vol. 75, issue 1, id. 016707, 2007[38] G. F. Smoot, et al., “Structure in the COBE differential microwaveradiometer first-year maps”,
Astrophysical Journal , vol. 396, L1 , 1992[39] R. Fabbri and S. Torres, “Peak statistics on COBE maps”,
Astronomyand Astrophysics , vol. 307, pp. 703-707, 1996[40] J. Timmer and M. Koenig, “On generating power law noise”,
Astronomyand Astrophysics , vol. 300, p. 707, 1995[41] P. Uttley and I. M. McHardy, “The flux-dependent amplitude of broad-band noise variability in X-ray binaries and active galaxies”,
MonthlyNotices of the Royal Astronomical Society , vol. 323, L26, 2001[42] W. H. Press, S. A. Teukolsky, W. T. Vetterling W.T. and B. P. FlanneryB.P., “Numerical recipes in FORTRAN. The art of scientific computing”,Cambridge: University Press, 2nd ed., 1992[43] M. van der Klis, “Quantifying Rapid Variability in Accreting CompactObjects”, Statistical Challenges in Modern Astronomy II, 321, 1997[44] Ch. Clapham, The Concice Oxford Dictonary of Mathematics, OxfordUniversity Press, 1996[45] D. A. Smith, “XTE J1550-564”, IAU Circ., 7008, 1, 1998[46] J. A. Orosz, et al. “Dynamical Evidence for a Black Hole in theMicroquasar XTE J1550-564”,
The Astrophysical Journal , vol. 568, p.845, 2002[47] W. Cui, S. N. Zhang, W. Chen and E. H. Morgan, “Strong aperiodicX-ray variability and quasi-periodic oscillationin X-ray nova XTE J1550-564,
The Astrophysical Journal Letters , vol. 512, L43, 1999[48] R. A. Remillard, J. E. McClintock, G. J. Sobczak, C. D. Bailyn,J. A. Orosz, E. H. Morgan and A. M. Levine, “X-Ray Nova XTEJ1550-564: Discovery of a Quasi-periodic Oscillation near 185 Hz”,