text:
# Fourier Techniques for Very Long Astrophysical Time Series Analysis

## 1 Introduction

The analysis of time series data is an important tool in many areas of astrophysics, including research involving white dwarfs, black holes, and neutron stars. In the study of neutron stars, time series analysis has had particular importance for pulsar research due to the high coherence of pulsar periodicities. While many techniques can be used to investigate the properties of these periodic signals, the computational efficiency of the Fast Fourier Transform (FFT) makes Fourier analysis the preferred approach for many applications (see e.g. Burns & Clark 1969). The literature contains literally hundreds of descriptions of various aspects of Fourier analysis, most of which deal with signal detection using the power spectrum.

Because of this concentration on power spectra, many researchers discard a wealth of information provided by the Fourier phases. Techniques which use this phase information exist that can provide insight into many useful signal properties. While many of these techniques have been known for some time (see e.g. Middleditch 1976), few have appeared in textbooks or refereed journals, and fewer still have been presented with any sort of derivation or insight into their assumptions and/or limitations.

A second, more practical problem with most astronomical Fourier analysis is its concentration on short time series. We define “short” to mean either that the full time series of binned or sampled data can fit into the core memory of your computer ($`N10^7`$ points), or for data consisting of events (such as photon arrival times in X-ray astronomy), that the time resolution ($`dt`$) of each event makes up a non-negligible fraction of the total time duration ($`T`$) of the data ($`T/dt10^7`$). Billion point FFTs have been used successfully in the past (e.g. Anderson 1992; Mattox et al. 1996; Middleditch et al. 2000), but each required the use of state-of-the-art supercomputing facilities. Today, such analyses are possible using clusters of workstations or even individual desktop machines. Many projects, such as pulsar searches of globular clusters, astero- or helioseismological observations, and gravitational wave experiments, require extremely large Fourier transforms in order to make the highest sensitivity (i.e. *coherent*) searches for pulsations, and to extract the maximum amount of information from signals found in these searches.

It is our goal in this paper to present some useful Fourier analysis techniques, that for various reasons are used only rarely when working with long time series. Most of our examples pertain to pulsar searches of very long time series, but the methods can be used in the Fourier analysis of virtually any coherent periodicity. This paper will briefly discuss the properties of the DFT, its response to periodic signals and noise, and methods for its computation. We will discuss methods for interpolating Fourier amplitudes, estimating a signal’s duration and centroid in time, accurately determining its frequency and frequency derivative, correcting for changing pulsation frequency during an observation, and estimating the phase- or amplitude modulation of a signal. Such techniques allow the detection of signals whose frequencies or amplitudes change significantly during an observation — for instance, due to orbital motion or pulsar spin-down.

## 2 Discrete Fourier Transform

In order to present more advanced Fourier techniques later, we first review some fundamentals of the DFT and its most common implementation, the FFT. Since thorough discussions of the Fourier transform in both its continuous and discrete variants exist in the literature (e.g. Bracewell 1999), we will mention only a few topics of particular relevance to astrophysical data analysis, closely following Middleditch (1976).

### 2.1 Introduction to the DFT

The $`k^{\mathrm{th}}`$ element of the Discrete Fourier Transform (DFT), of a uniformly-spaced time series $`n_j`$ ($`j=0,1,\mathrm{},N1`$) is defined as

$$A_k=\underset{j=0}{\overset{N1}{}}n_j\mathrm{e}^{2\pi \mathrm{i}jk/N},$$
(1)
where $`\mathrm{i}=\sqrt{1}`$ and $`k`$ is the Fourier frequency or wavenumber ($`k=0,1,\mathrm{},N1`$). For a time spacing $`dt`$ between successive data elements, the frequency of the $`k^{\mathrm{th}}`$ element is $`f_k=k/(Ndt)=k/T`$, where $`T`$ is the total time duration of the sequence being transformed. This frequency spacing of $`1/T`$, for evenly sampled, un-padded and non-windowed data, is often called the Independent Fourier Spacing (IFS)111The IFS is important when trying to determine the overall significance of a candidate from a search. The number of IFSs searched corresponds to the number of independent trials searched and should therefore be included in calculations that try to determine if a candidate was produced by noise. See Vaughan et al. (1994). , and defines the finest frequency resolution available while maintaining completely independent Fourier amplitudes. Fourier frequency $`N/2T`$ is known as the Nyquist frequency.

If we view the DFT summation in the complex plane, we see that it is a simple vector addition with each element rotated by $`2\pi k/N`$ from the previous element. If the $`n_j`$ have a constant value, the sum will form $`k`$ regular $`N/k`$-sided polygons with each polygon returning near to the origin, and with the last one returning exactly to the origin. Therefore, the DFT of a constant data string will be identically zero for all $`k>0`$, and equal to the sum of the data elements for $`k=0`$ (the “DC” frequency element).

For most astrophysical observations, the data points are real-valued. This property adds an important feature to Fourier transforms of these time series — they are symmetric about the Nyquist frequency, such that $`A_{Nk}=A_k^{}`$, where $`A_k^{}`$ represents the complex conjugate of $`A_k`$. This symmetry allows us to calculate the full DFT of a time series by computing amplitudes at only half of the normal Fourier frequencies, thereby speeding up computation of the DFT by nearly a factor of two (e.g., Press et al. 1992).

When deriving properties and techniques based on DFTs, it is often both computationally and intuitively easier to work with a time-normalized data series, where $`T=1`$ and $`u`$ represents the fraction of the observation complete at any given instant (such that $`0u1`$). In this case, instead of working with frequencies $`f`$, or integral wavenumbers $`k`$, we define $`r`$, which represents any real wavenumber (including non-integers). If the number of samples, $`N`$, from our data source, $`n(u)`$, is large, we can compute a continuous FT

$$A_r=N_0^1n(u)\mathrm{e}^{2\pi \mathrm{i}ru}du,$$
(2)
which is almost identical to our DFT when $`r=k`$, and produces very high accuracy estimates of the Fourier amplitudes at any frequency such that $`0rN/2`$. We will use this approximation in many of our derivations.

### 2.2 Computation of Very Long FFTs

The Fast Fourier Transform is a family of well-known computer algorithms which quickly calculate a DFT in only $`O(N\mathrm{log}_2(N))`$ operations as opposed to the $`O(N^2)`$ operations of a brute-force DFT. FFTs, their computation, and their origins have been described in numerous articles and books over the last few decades (see Bracewell 1999, for an introduction). Therefore, we will describe only a few special versions which have become generally known only recently and which are useful in the analysis of extremely long time series.

High energy pulsar searches using photon counting systems (infrared, optical, x-ray, and gamma-ray detectors) or pointed radio telescope searches (for example, searching globular clusters) often utilize very high sampling rates (i.e. $`1050`$ kHz) and very long integration times of hours, days, or even weeks. These observations result in time series with hundreds of millions or even billions of data points. The subsequent FFTs are impossible to perform using conventional FFT algorithms unless the full time series fits into the core memory of the computer being used. Utilizing special “out-of-core” FFT methods, we can transform such huge arrays using distributed memory parallel systems or with a single workstation and manageable numbers of passes through the data stored on external media. Most of these methods are based on the “four-step” FFT and/or special external media array permutation methods (Fraser 1976; Bailey 1990).

#### 2.2.1 The “Four-Step” FFT

Closely following the derivation described in Sweet & Wilson (1995), we can think of our one-dimensional time series as a “C-like” or row-ordered, two-dimensional matrix of size $`N=N_rN_c`$, where $`N_c`$ is the number of columns (i.e. the length of a row) and $`N_r`$ the number of rows (i.e. the length of a column). Using this data decomposition, the FFT is computed by *1.* FFTing the columns of the matrix, *2.* multiplying the data by complex “twiddle factors”, *3.* FFTing the rows of the matrix, and *4.* matrix transposing the result. If we define indices $`x=0,1,\mathrm{},N_c1`$, $`y=0,1,\mathrm{},N_r1`$, $`l=0,1,\mathrm{},N_r1`$, and $`m=0,1,\mathrm{},N_c1`$, we can write our signal and its DFT amplitudes as

$$n(x,y)=n_j,j=N_cy+x$$
(3)
and

$$A(l,m)=A_k,k=N_rm+l.$$
(4)
Substituting into the definition of the DFT (eqn. 1) and simplifying we get

$$A(l,m)=\underset{x=0}{\overset{N_c1}{}}\left[\mathrm{e}^{2\pi \mathrm{i}xl/N}\underset{y=0}{\overset{N_r1}{}}n(x,y)\mathrm{e}^{2\pi \mathrm{i}yl/N_r}\right]\mathrm{e}^{2\pi \mathrm{i}xm/N_c}.$$
(5)
Notice that the bracketed terms are really $`N_c`$ DFTs of length $`N_r`$ — FFTs of the matrix columns multiplied by the “twiddle factor” $`\mathrm{e}^{2\pi \mathrm{i}xl/N}`$. We denote these column FFTs as

$$A_c(x,l)=\mathrm{e}^{2\pi \mathrm{i}xl/N}\underset{y=0}{\overset{N_r1}{}}n(x,y)\mathrm{e}^{2\pi \mathrm{i}yl/N_r}.$$
(6)
We now see that the outer summation is $`N_r`$ DFTs of length $`N_c`$ — FFTs of the matrix rows composed of the $`A_c(x,l)`$ terms. We denote the result of these transforms as

$$A_r(m,l)=\underset{x=0}{\overset{N_c1}{}}A_c(x,l)\mathrm{e}^{2\pi \mathrm{i}xm/N_c}.$$
(7)
To recover the full FFT in its normal order, we simply transpose $`A_r(m,l)`$

$$A(l,m)=A_r^T(m,l).$$
(8)
The “four-step” algorithm only needs small portions of the data in memory at any one time. Unfortunately, the short FFTs are strided in memory222“Strided” means that the data is stored in a sequence of non-contiguous memory locations spaced by a constant amount of memory known as the “stride”. by instead of being stored contiguously. This results in a significant inefficiency in today’s cache-based processors and requires inter-node communications when distributing the short FFTs over many processors on parallel computer systems. These shortcomings can be overcome with the “six-step” algorithm (Bailey 1990).

#### 2.2.2 The “Six-Step” FFT

If the initial data set is transposed from a $`N_c\times N_r`$ matrix into a $`N_r\times N_c`$ matrix, the strided column FFTs of length $`N_r`$ become contiguous row FFTs. This memory locality facilitates the use of processor cache systems, greatly increasing memory response times, and allows independent calculation of the row FFTs by different processors in parallel systems. A second transpose operation after the row FFTs have been calculated counteracts the effects of the first transpose operation and makes the next set of row FFTs contiguous in memory as well. The FFT is finished with a final transpose operation which brings the data into normal order.

These transpose operations are relatively efficient on distributed memory parallel systems with fast inter-node communications, and add little time to the overall FFT. In fact, for applications with very large $`N`$, the time required to move the data from external media to computer memory and vice versa tends to dominate the FFT time. The “six-step” algorithm has become the “standard” FFT algorithm for distributed memory systems where the serial nature and large number of short FFTs exploit parallel computation strengths.

Unfortunately, if the data does not all fit into the core memory of a workstation or parallel machine, the transposes become extremely slow operations since data must be read to and written from much slower external media. Fraser (1976) devised optimized methods for dealing with such data permutations on external media including the “two-pass” FFT algorithm.

#### 2.2.3 The “Two-Pass” Out-of-Core FFT

Fraser (1976) and Bailey (1990) both describe how a very large data set may be Fourier transformed with only two read-write passes through externally stored data if a scratch area on the external media the same size as the input data set exists. The method uses the “four-step” algorithm with out-of-core transposes. These transpose algorithms allow blocked external media access and perform most of the transposition work in core memory. While external media access speeds and transfer rates are orders of magnitude slower than core-memory systems, the “two-pass” algorithm allows huge arrays to be transformed in manageable times.

## 3 Fourier Transforms of Real Data

### 3.1 Fourier Response to Noise

If our data $`n_j`$ are composed of some constant value $`c_j`$ plus random real-valued noise fluctuations $`d_j`$ with zero mean, the transform terms become

$$(c_j+d_j)\mathrm{e}^{2\pi \mathrm{i}jk/N}=c_j\mathrm{e}^{2\pi \mathrm{i}jk/N}+d_j\mathrm{e}^{2\pi \mathrm{i}jk/N}.$$
(9)
The linear nature of the Fourier transform allows us to treat the DFT of the $`d_j`$ independently from the constant length steps, $`c_j`$. Since the complex phase factor for a given $`j`$ and $`k`$ is fixed, the direction of each element in the sum is nearly fixed. However, since the sign of the $`d_j`$ may be either positive or negative, the vector direction of the $`j^{\mathrm{th}}`$ element may be reversed. Thus, the DFT of the $`d_j`$ creates a kind of random walk in the complex plane.

The statistical properties of this random walk for DFTs of pure noise have been well studied (see, e.g. Blackman & Tukey 1959), and result in power spectra distributed according to an exponential distribution (a $`\chi ^2`$ distribution with 2 degrees of freedom) with average and standard deviation equal to $`Nd_j^2`$. If we normalize the powers by dividing by $`Nd_j^2`$, the probability for a power $`P=\left|A_k\right|^2`$ in a single bin to equal or exceed a power $`P^{}`$ by chance is333This is different than the probability for an actual signal to produce a power $`P>P^{}`$ in the presence of noise (see §3.3). This difference is important in setting upper limits on the amplitudes of periodic signals as discussed in Vaughan et al. (1994).

$$Prob(PP^{})=\mathrm{e}^P^{}.$$
(10)
Similarly, if we sum $`m`$ properly normalized powers, the probability for the summed power $`P_m`$ to exceed a power $`P^{}`$ is given by

$$Prob(P_mP^{})=\underset{j=0}{\overset{m1}{}}\frac{(P^{})^j}{j!}\mathrm{e}^P^{},$$
(11)
which is the probability for a $`\chi ^2`$ distribution of $`2m`$ degrees of freedom to exceed $`2P^{}`$. Such an *incoherent* (since no phase information is used) summation of powers is often useful when searching for signals suspected of having power in many harmonics (see §3.2.2).

Proper normalization of the powers is essential for an accurate estimate of a signal’s statistical significance or lack thereof. We often cannot normalize our power spectrum by simply dividing by $`Nd_j^2`$, since frequency-dependent noise may be present throughout our power spectrum — perhaps as a result of instrumental, atmospheric, or astrophysical processes. Typically, these processes produce noise which increases in strength towards the low-frequency part of the spectrum and is correspondingly called *red noise*.

Techniques to flatten or remove this “coloured” noise component from the power spectrum are described by Israel & Stella (1996), and usually involve dividing short portions of the power spectrum by the locally determined average power level, $`P_{local}`$, such that

$$P_{k,norm}=\frac{\left|A_k\right|^2}{Nd_j^2}\frac{\left|A_k\right|^2}{P_{local}}=\frac{P_k}{P_{local}}.$$
(12)
As long as the number of averaged powers is small enough such that the power spectrum is roughly constant over the range in question, a *white noise* like power spectrum is produced with average and standard deviation of approximately one and an exponential distribution (eqn. 10).

Since strong narrow-band signals near some frequency of interest will skew a local power average upwards (and correspondingly decrease the calculated significance of a signal detection), it is important to exclude such powers from the calculation of $`P_{local}`$. A simple and effective way to accomplish this is by normalizing with a corrected local median power level instead of the local power average. An exponential distribution with unity mean and standard deviation has a median of $`\mathrm{ln}(2)`$. Therefore, dividing a section of raw powers by $`1/\mathrm{ln}(2)`$ times the local median value is theoretically equivalent to normalizing with the local mean, but has the advantage of being insensitive to high-power outliers in the spectrum.

More advanced algorithms for the removal of “coloured” noise and power normalization do exist. A simple example involves fitting polynomial models to portions of the power spectrum and then dividing by them out. These methods work well for Fourier frequencies near zero where the assumption of roughly equivalent power levels for the local powers may be unwarranted.

For the special case where the noise in our data is purely Poissonian (i.e. for binned photons in an optical or x-ray observation), we have $`d_j^2=n_j`$. In this case, our properly normalized power for the $`k^{\mathrm{th}}`$ DFT element is

$$P_{k,norm}=\frac{\left|A_k\right|^2}{Nd_j^2}=\frac{\left|A_k\right|^2}{n_{ph}},$$
(13)
where $`n_{ph}=Nn_j`$ is the sum of the $`n_j`$ (or the total number of photons for a photon-counting system), which is also equal to the “DC” frequency value of the FT (see §2.1). However, the same processes that caused the “coloured” noise discussed above can significantly alter this situation and require a power normalization based on local powers (see eqn. 12).

### 3.2 Fourier Response to Periodic Signals

One of the more useful properties of the FT for astronomical purposes is its response to periodic signals. Since all real periodic signals can be expanded into a series of sinusoids it is important to understand the FT response to a simple sine wave.

#### 3.2.1 Sinusoidal Signals

If we now let our $`n_j`$ represent a sampled cosinusoid of amplitude $`a`$, phase $`\varphi `$, and frequency $`f_r=r/T`$ (where wavenumber $`r`$ is an integer and $`f_r`$ an “integral frequency”), we can write

$`n_j`$ $`=`$ $`a\mathrm{cos}(2\pi f_rjdt+\varphi )`$ (14a)
$`=`$ $`a\mathrm{cos}(2\pi jr/N+\varphi )`$ (14b)
$`=`$ $`{\displaystyle \frac{a}{2}}\left(\mathrm{e}^{2\pi \mathrm{i}jr/N+\mathrm{i}\varphi }+\mathrm{e}^{2\pi \mathrm{i}jr/N\mathrm{i}\varphi }\right).`$ (14c)
From this expression, we see that the $`k^{\mathrm{th}}`$ element of the DFT is given by

$`A_k`$ $`=`$ $`{\displaystyle \frac{a}{2}}{\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{e}^{2\pi \mathrm{i}jk/N}\left(\mathrm{e}^{2\pi \mathrm{i}jr/N+\mathrm{i}\varphi }+\mathrm{e}^{2\pi \mathrm{i}jr/N\mathrm{i}\varphi }\right)`$ (15a)
$`=`$ $`{\displaystyle \frac{a}{2}}{\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{e}^{2\pi \mathrm{i}j(kr)/N+\mathrm{i}\varphi }+\mathrm{e}^{2\pi \mathrm{i}j(k+r)/N\mathrm{i}\varphi },`$ (15b)
and represents the summation of two vectors in the complex plane. For $`kr`$, the first term traces out $`\left|kr\right|`$ complete “rotations” (pseudo-polygons which start and end at the origin) in the complex plane (since $`kr`$ is an integer), giving a net contribution of zero to the $`k^{\mathrm{th}}`$ DFT element. The second term traces out $`k+r`$ complete rotations and once again contributes nothing to the $`k^{\mathrm{th}}`$ DFT element (since $`k`$ and $`r`$ are both positive and therefore $`k+r0`$.

When $`k=r`$, however, the consecutive terms in the summation add coherently (i.e. in phase and therefore without rotation), since the rotation caused by the DFT exponential exactly cancels that from the signal exponential (the Fourier transform “derotates” the signal). As a result, each element of the sum is a step in the complex plane of magnitude $`a/2`$ in a direction parallel to the one set by the arbitrary initial phase of the signal, $`\varphi `$. For a cosinusoidal signal with integral frequency $`f_r`$, the DFT will be uniformly zero, except in the $`r^{\mathrm{th}}`$ frequency bin, where the response is

$$A_r=\frac{Na}{2}\mathrm{e}^{\mathrm{i}\varphi }.$$
(16)
The Fourier response is more complicated for sinusoids with non-integral frequencies (i.e. wavenumber $`r`$ is a non-negative real number). The $`k^{\mathrm{th}}`$ DFT element is still given by eqn. 15b, but not all of the signal ends up in a single DFT bin $`A_k`$. When $`k=[r]`$ (where $`[r]`$ is the nearest integer to $`r`$), the first term in eqn. 15b traces out a fraction ($`kr`$) of a complete rotation in the complex plane, while the second term traces out $`k+[r]`$ complete rotations, plus an additional fractional rotation.

When $`N`$ is large, these complete and fractional “rotations” can be treated as circles and arcs respectively. Therefore, the first term of eqn. 15b results in a semi-circular arc of length $`Na/2`$ along the arc, while the second term produces a semi-circular arc of length $`a/2(N\mathrm{mod}k+[r])`$ along the arc. The DFT response is simply a vector drawn from the origin to the end of the arc (see Fig. 1). Since virtually all astrophysical applications involve $`r1`$, where the first term dominates the response, we will ignore the second term in the rest of our analysis.

The final response is a chord subtending $`2\pi (kr)`$ radians of a circle of radius $`Na/4\pi (kr)`$. The equation for the length of a chord is

$$C=2\frac{s}{\mathrm{\Theta }}\mathrm{sin}\left(\frac{\mathrm{\Theta }}{2}\right),$$
(17)
where $`s`$ is the arc length and $`\mathrm{\Theta }`$ is the angle subtended by the chord. The curvature of the arc away from the signal’s starting phase $`\varphi `$ results in a phase change of $`e^{\mathrm{i}\pi (kr)}`$. Therefore, the DFT response and power are

$`A_k`$ $`=`$ $`\mathrm{e}^{\mathrm{i}\varphi }\mathrm{e}^{\mathrm{i}\pi (kr)}\mathrm{\hspace{0.17em}2}{\displaystyle \frac{Na/2}{2\pi (kr)}}\mathrm{sin}\left[{\displaystyle \frac{2\pi (kr)}{2}}\right]`$ (18a)
$`=`$ $`{\displaystyle \frac{Na}{2}}\mathrm{e}^{\mathrm{i}\varphi }\mathrm{e}^{\mathrm{i}\pi (kr)}{\displaystyle \frac{\mathrm{sin}[\pi (kr)]}{\pi (kr)}}`$ (18b)
$`=`$ $`A_o\mathrm{e}^{\mathrm{i}\pi (kr)}\mathrm{sinc}\left[\pi (kr)\right]`$ (18c)
$`P_k`$ $`=`$ $`\left|A_k\right|^2=P_o\mathrm{sinc}^2\left[\pi (kr)\right],`$ (18d)
where $`A_o=Na/2\mathrm{e}^{\mathrm{i}\varphi }`$ is the DFT response for an integral frequency signal (eqn. 16), $`P_o`$ is the corresponding Fourier power, and the $`\mathrm{sinc}`$ function is defined as $`\mathrm{sinc}(x)=\mathrm{sin}(x)/x`$. This result is easily confirmed by a direct integration of eqn. 2 where $`n(u)`$ is equal to eqn. 14c with $`j/Nu`$.

The $`\mathrm{sinc}`$ factor in eqn. 18c produces a loss of sensitivity for the standard FFT to most real-world signals (where $`r`$ is not an integer). This effect, often called “scalloping” (e.g. Middleditch et al. 1993), is shown in Fig. 2, and causes a worst-case (when $`\left|kr\right|=1/2`$) amplitude reduction of $`\left|A_k\right|=2\left|A_o\right|/\pi `$ — a nearly $`60\%`$ loss of signal power. On average, scalloping results in a $`23\%`$ loss of signal power (van der Klis 1989). It is important to remember, though, that this loss in sensitivity is due to the finite frequency resolution of the FFT algorithm rather than an intrinsic feature of the data itself. In §4.1 we discuss various methods to reduce or even eliminate this loss of sensitivity.

#### 3.2.2 Non-Sinusoidal Signals

Many real-world periodic signals are not sinusoidal. Fortunately, we can expand all real-valued pulsations as a series of $`m`$ sinusoidal components

$`n_j`$ $`=`$ $`{\displaystyle \underset{h=1}{\overset{m}{}}}a_h\mathrm{cos}(2\pi jhr/N+\varphi _h)`$ (19a)
$`=`$ $`{\displaystyle \underset{h=1}{\overset{m}{}}}{\displaystyle \frac{a_h}{2}}\left(\mathrm{e}^{2\pi \mathrm{i}jhr/N+\mathrm{i}\varphi _h}+\mathrm{e}^{2\pi \mathrm{i}jhr/N\mathrm{i}\varphi _h}\right)`$ (19b)
where $`h=1,2,\mathrm{},m`$ specifies the harmonic number (with $`h=1`$ known as the “fundamental”), and $`a_h`$ and $`\varphi _h`$ represent the amplitude and phase of each component respectively. Due to the linear nature of the FT, we can treat the harmonics as independent sinusoidal signals. Each of these sinusoids produces a Fourier response equivalent to eqn. 18c, except that $`A_o`$ becomes $`A_h=Na_h/2\mathrm{e}^{\mathrm{i}\varphi _h}`$.

For nearly sinusoidal pulsations only the first few terms of eqn. 19a contain significant amplitudes, $`a_h`$. This results in a similarly small number of significant peaks in the corresponding power spectrum of the data. Low duty-cycle pulsations (i.e. those with a pulse that is short compared to the pulse period), such as most radio pulsars, on the other hand, have dozens of significant terms in their expansions and therefore harmonics in their power spectra.

A useful pulsation model, particularly for radio and x-ray pulsars, can be constructed based on a modified von Mises distribution (MVMD)

$$f(\kappa ,t)=a\frac{\mathrm{e}^{\kappa \mathrm{cos}(2\pi f_rt+\varphi )}\mathrm{e}^\kappa }{\mathrm{I}_0(\kappa )\mathrm{e}^\kappa },$$
(20)
where $`0tT`$ is the instantaneous time, $`\mathrm{I}_0`$ is the modified Bessel function of zeroth order, and the shape parameter, $`\kappa `$, determines the width of the function (e.g. Mardia & Zemroch 1975). In the limit as $`\kappa 0`$, the MVMD becomes a sinusoid, while as $`\kappa \mathrm{}`$, it becomes a Gaussian (see Fig. 3). The integral of the MVMD over a single pulse period is simply $`a`$, all of which is pulsed (i.e. the pulsed fraction is one). The full-width at half-maximum (FWHM) as a fraction of a pulse, is

$$\mathrm{FWHM}_{\mathrm{MVMD}}=\pi ^1\mathrm{arccos}\left\{\mathrm{ln}\left[\mathrm{cosh}(\kappa )\right]\right\},$$
(21)
and the maximum value is

$$\mathrm{max}_{\mathrm{MVMD}}=\frac{2a\mathrm{cosh}(\kappa )}{\mathrm{I}_0(\kappa )\mathrm{e}^\kappa }.$$
(22)
The FT of the MVMD can be computed in a particularly convenient form for harmonic analysis. According to Abramowitz & Stegun (1972, eqn. 9.6.34), we can expand the exponential in the MVMD as

$$\mathrm{e}^{\kappa \mathrm{cos}(x)}=\mathrm{I}_0(\kappa )+2\underset{h=1}{\overset{\mathrm{}}{}}\mathrm{I}_h(\kappa )\mathrm{cos}(hx),$$
(23)
where $`\mathrm{I}_h`$ is the modified Bessel function of order $`h`$. When combined with the rest of the MVMD definition we have

$$f(\kappa ,t)=a+\frac{2a_{h=1}^{\mathrm{}}\mathrm{I}_h(\kappa )\mathrm{cos}(2\pi hf_rt+h\varphi )}{\mathrm{I}_0(\kappa )\mathrm{e}^\kappa }.$$
(24)
This expression is simply a “DC” term (since the integral over a pulse equals $`a`$) plus a series of independent cosinusoidal harmonics. After Fourier transforming (i.e. substituting into eqn. 2 with $`f_rtru`$), we are left with a series of harmonics of amplitude $`AN\mathrm{I}_h(\kappa )/\left[\mathrm{I}_0(\kappa )\mathrm{e}^\kappa \right]`$ and phase $`h\varphi `$ at Fourier frequency $`hr`$. It is important to note that as $`\kappa \mathrm{}`$ and the pulse becomes narrower, the Fourier amplitudes of the low order harmonics are twice that of a sinusoid with the same pulsed fraction (see the dashed line in Fig. 4). This fact, along with the large number of harmonics that low duty-cycle pulsations generate, can significantly increase search sensitivities to such pulsations (see Fig. 5).

Fig. 4 shows the approximate number of significant harmonics (meaning that a harmonic’s amplitude is greater than one-half the amplitude of the fundamental) generated by an MVMD pulsation, as well as a histogram of the duty-cycles of over 600 radio pulsars (the majority of which are from Taylor et al. 1995). Most radio pulsars have duty-cycles $`5\%`$ corresponding to $`10`$ significant harmonics — assuming a sufficient data sample rate.

### 3.3 Periodic Signals With Noise

When a periodic signal is present in a noisy time series, a sum of $`m`$ powers $`P_m`$, containing some amount of signal power $`P_s`$, is no longer described by a $`\chi ^2`$ distribution with $`2m`$ degrees of freedom (see §3.1). Groth (1975) calculated the expectation value and variance of $`P_m`$ as $`P_m=m+P_s`$ and $`P_m^2P_m^2=m+2P_s`$ respectfully. He also derived the exact probability density function for $`P_m`$ which can be integrated to give the probability that $`P_m`$ is greater than or equal to some power $`P^{}`$,

$$Prob((P_m;P_s)P^{})=\mathrm{e}^{\left(P^{}+P_s\right)}\underset{k=0}{\overset{\mathrm{}}{}}\underset{j=0}{\overset{k+m1}{}}\frac{(P^{})^jP_s^k}{j!k!}.$$
(25)
When $`P_s=0`$ this equation reduces to Eqn. 11.

The fact that the probability density function for a signal plus noise is different from a $`\chi ^2`$ distribution with $`2m`$ degrees of freedom is very important when trying to determine the sensitivity of a search for pulsations or an upper limit to the amplitude of a periodic signal present in a time series. Vaughan et al. (1994) describe a procedure444Note that Vaughan et al. (1994) use a power normalization that is a factor of two higher than ours. for correctly determining search sensitivities and upper limits using the equations of Groth (1975).

### 3.4 Photon Counting Data

Since many of today’s astronomical time series come from photon counting experiments, it is important to raise some of the issues particular to Fourier analysis of such data. If we can assume purely Poissonian statistics, a power spectrum of pure noise is flat, and can be normalized simply by dividing by the total number of photons in the data (the zeroth or “DC” frequency bin from the DFT — see §3.1). In addition to this difference in power spectrum normalization, the other points worth noting come from the fact that photon counting data is based on the measurement of events rather than the instantaneous sampling of a continuous process.

One important issue that is beyond the scope of this paper is dead-time correction. Dead-time effects modify a detector’s sensitivity to photons for some time after the detection of an earlier photon. These effects can cause complicated non-linear and frequency-dependent systematics during Fourier analysis. We refer the reader to Zhang et al. (1995) and references therein, for a thorough discussion of this topic.

#### 3.4.1 Binned vs. Sampled Data

Many high-energy telescopes and detectors produce time series of binned photons rather than the sampled data produced by radio telescopes. Since binning essentially averages a periodic signal’s instantaneous rate over the binning time ($`dt`$), it modifies the Fourier response to the signal. Binning removes phase information from the data and causes the Fourier response to sinusoidal pulsations to become frequency-dependent — resulting in decreased sensitivity at high frequencies (e.g. Middleditch 1976; Leahy et al. 1983).

The frequency-dependent loss in Fourier amplitude due to binning is $`\mathrm{sinc}(\pi f_rdt)`$ or $`\mathrm{sinc}(\pi r/N)`$. The binned data Fourier response to a sinusoid is therefore eqn. 18c times this factor. This decrease in sensitivity corresponds to a loss in signal power of about $`19.8\%`$ at half the Nyquist frequency and $`59.5\%`$ at the Nyquist frequency itself.

For Poissonian noise (i.e. from a photon-counting experiment that does not introduce count-rate dependent systematics), which is independent of the sampling interval, the Fourier response is flat over all frequencies. This is in contrast to a sinusoidal signal passing through the same system which suffers the frequency-dependent attenuation described above(Middleditch 1976). Such behavior is important when trying to estimate limits or amplitudes for pulsations in a time series (Vaughan et al. 1994).

#### 3.4.2 Low Count Rate Data

The Fourier analysis of gamma-ray or x-ray observations often places us in a very unique regime — very long integration times ($`10^4`$ s) with very low numbers of counts ($`10^3`$ photons). In addition, due to visibility constraints based on the orbits of the telescopes, large fractions of the time between the first and last photons may be devoid of counts.

Fourier analysis of such data can overwhelm present computational resources. For example, a $`10^6`$ s observation (about 11.6 days) with photon time-of-arrival (TOA) resolution of $`10^4`$ s would require a 10 Gpt FFT for a full-resolution analysis. Such FFTs, while possible, are extremely difficult to compute unless very special and dedicated hardware resources are available. If this data contains only a small number of photons, however, we can exactly compute the DFT over any frequency range and to any frequency resolution using a brute-force implementation of the FT.

If we treat each TOA as a sample of amplitude one, an exact DFT amplitude at arbitrary Fourier frequency $`r`$ becomes

$$A_r=\underset{j=1}{\overset{n_{ph}}{}}\mathrm{e}^{2\pi \mathrm{i}r(t_jt_o)/T},$$
(26)
where $`n_{ph}`$ is the number of photons, $`t_j`$ is the TOA of the $`j^{th}`$ photon, $`t_o`$ is the time of the start of the observation, and $`T`$ is the total duration of the observation. Very quick harmonic-summing searches of an observation are possible using this technique, with the added benefit that “scalloping” (see §3.2.1) is non-existent.

Since eqn. 26 only involves a summation over the number of photons, it can be computed quickly if $`n_{ph}`$ is relatively small. Great increases in computation speed can be had if we search a regular grid of Fourier frequencies. Trigonometric recurrences such as

$`\mathrm{cos}(\theta +\delta )`$ $`=`$ $`\mathrm{cos}(\theta )\left[\alpha \mathrm{cos}(\theta )+\beta \mathrm{sin}(\theta )\right]`$ (27a)
$`\mathrm{sin}(\theta +\delta )`$ $`=`$ $`\mathrm{sin}(\theta )\left[\alpha \mathrm{sin}(\theta )\beta \mathrm{cos}(\theta )\right],`$ (27b)
where $`\alpha =2\mathrm{sin}^2\left(\delta /2\right)`$ and $`\beta =\mathrm{sin}(\delta )`$, allow extremely efficient calculation of the complex exponential for each TOA (Press et al. 1992). This technique allows one to calculate billions of Fourier frequencies from a few hundred photons using only modest computational resources.

## 4 Improving the DFT Response to Arbitrary Signals

### 4.1 Fourier Interpolation

The potential for up to a 30% decrease in signal-to-noise ($`S/N\sqrt{P}`$) due to an essentially arbitrary difference between the signal frequency and the integer frequency of the nearest Fourier bin is clearly a drawback in the use of the DFT (see §3.2.1). However, if we could calculate the Fourier Transform (FT) at a higher frequency resolution than the $`1/T`$ spacing that results from the FFT, we could significantly reduce or eliminate scalloping and effectively flatten the Fourier response as a function of frequency.

One possibility for increasing the frequency resolution is to simply calculate the DFT by a brute force summation at frequencies between the integer frequencies. Such a technique is possible in special situations (see §3.4.2), but for most applications, the computational costs would be unacceptably high. Another well-known possibility is to “pad” the end of the time series with a large number of points with values equivalent to the mean of the data555Padding with the data mean is preferable to zero-padding since zero-padding introduces low frequency power into the Fourier response.. The padding adds no power to the data but it does increase the Fourier resolution since $`T`$ has been artificially increased by the padding. While this technique is simple and effective for short time series, the difficulties involved in performing very long FFTs (§2.2) makes this technique difficult when dealing with long time series.

Yet another way to calculate a higher-resolution Fourier response is to use the complex amplitudes produced by the standard FFT to interpolate responses at non-integer frequencies — a process known as “fine-binning” or “Fourier interpolation” (e.g. Middleditch et al. 1993). Similar techniques allow the full recovery of a signal’s theoretical coherent response provided that the signal’s behavior during the observation is either known or can be guessed.

The purpose of Fourier interpolation is to calculate a complex Fourier amplitude at an arbitrary frequency $`f_r=r/T`$, where $`r`$ is any real number, such that the result is sufficiently close to the exact calculation,

$$A_r=\underset{j=0}{\overset{N1}{}}n_j\mathrm{e}^{2\pi \mathrm{i}jr/N}.$$
(28)
We can rewrite this expression as

$$A_r=\underset{k=0}{\overset{N1}{}}A_k\mathrm{e}^{\mathrm{i}\pi (rk)}\mathrm{sinc}\left[\pi (rk)\right].$$
(29)
where the $`A_k`$ are the complex FFT amplitudes at the integer frequencies $`l`$ (see Appendix A and §4.2.2 for a derivation and discussion of this result).

The $`\mathrm{sinc}`$ function in eqn. 29, provides the key to computing an accurate interpolated amplitude using a relatively small number of operations. Since $`\mathrm{sinc}[\pi (rk)]0`$ as $`\pi (rk)\pm \mathrm{}`$, the expansion of $`A_r`$ in terms of the $`A_k`$ is dominated by the local Fourier amplitudes (i.e. where $`kr`$). We can therefore approximate $`A_r`$ as

$$A_r\underset{k=[r]m/2}{\overset{[r]+m/2}{}}A_k\mathrm{e}^{\mathrm{i}\pi (rk)}\mathrm{sinc}\left[\pi (rk)\right],$$
(30)
where $`[r]`$ is the nearest integer to $`r`$, and $`m`$ is the number of neighboring FFT bins used in the interpolation. Note that the interpolation is simply a correlation of the local FFT spectrum around the desired frequency element with a “template” response — in this case the theoretical response of a DFT to a sinusoid as described by eqn. 18c.

The upper panel in Fig. 6 shows both the raw FFT power spectrum (denoted by grey dots) and the interpolated power spectrum (the line connecting the dots) for a radio observation of the short-period binary pulsar PSR J1807$``$2459. The spectra cover a narrow frequency range near the pulsar’s rotational frequency and were calculated using $`m=32`$ and a frequency step size of $`\mathrm{\Delta }r=1/16`$ (compared to the raw FFT frequency step size of $`\mathrm{\Delta }r=1`$). Note that the interpolated spectrum reveals much more information regarding the shape of the frequency profile — including the true amplitude and location of the maximum in Fourier power. We will see in §5 how this can be used to deduce further information regarding the signal properties.

A computationally less expensive version of Fourier interpolation is “interbinning,” where we approximate the FT response at half-integer frequencies using only the nearest 2 integer frequency bins. By using the Fourier interpolation equation (eqn. 30) with $`m=2`$, ignoring an overall phase shift, and boosting the response such that the best-case response (at half-integer frequencies) is equivalent to the full response, we get

$$A_{k+\frac{1}{2}}\frac{\pi }{4}(A_kA_{k+1}).$$
(31)
This particular formulation of interbinning was reported by van der Klis (1989), and its response is shown in Fig. 2. Middleditch et al. (1993) have contributed a correction to this formula for use when the data are padded at the end.

Interbinning is extremely useful since such a computationally inexpensive calculation reduces the maximum loss of signal-to-noise from $`12/\pi `$ or $`36\%`$ at a frequency offset of ½ bin to $`7.4\%`$ at an offset of $`\pm \left(1\pi /4\right)`$ bins. This large but cheaply-obtained reduction in scalloping can be extremely beneficial when searching large numbers of FFT bins and interbins.

It is important to note that interbins as defined above have three different properties than integer FFT bins. First, they have different noise properties, which makes calculation of the significance of interbin powers much more difficult. Second, each interbin is correlated with the integer bins it was created from, meaning that interbins are not independent Fourier trials (see §2.1 for a discussion of the IFS). And finally, interbins do not recover the correct phase of a sinusoid at the interbin frequency. In general, since interbins are most commonly used during searches to simply identify signals in the power spectrum that would otherwise have been lost due to scalloping, these weaknesses do not degrade the usefulness of their calculation. When a signal is identified, a full-scale interpolation of the Fourier amplitudes around the signal using eqn. 30 allows accurate estimates to be made of the signal’s significance and other properties (see §5).

### 4.2 General DFT Response Correction

Fourier interpolation serves as a specific example of a much more general technique — the ability to completely recover the fully coherent response for virtually any signal. For Fourier interpolation, we can exactly calculate the response of any Fourier frequency based purely on the properties of the FT. To correct the Fourier transform’s response to a particular signal, we must know not only the properties of the FT, but the properties of the signal we are looking for as well. For the cases we will discuss, this ability comes in one of two forms: matched filtering in the Fourier domain using only the “local” Fourier amplitudes near the Fourier frequencies of interest (which we call the “correlation technique”), or the straightening of the curved Fourier vector addition in the complex plane (which we will call “vector bending”).

#### 4.2.1 Correcting for Constant Frequency Derivative

In order to illustrate these two methods, we demonstrate how to correct for a signal whose response is reduced due to a constant frequency derivative $`\dot{f}`$ (or in Fourier frequency bins, $`\dot{r}=\dot{f}T^2`$). The DFT operates, as we noted in §3.2.1, by “derotating” the vector addition of the data in the complex plane by changing the phases of each of the vector elements — causing a straight line to form for a sinusoidal signal with integral frequency. In the presence of a frequency derivative however, the signal frequency may change by one or more frequency bins over the course of the observation. The complex phase corrections provided by the DFT will fail to completely derotate the data, and pulsation power will be “smeared” across several nearby frequency bins — causing a decrease in the measured DFT response (e.g. Johnston & Kulkarni 1991). Fig. 1 illustrates this effect in the complex plane.

As with a frequency error, an uncorrected frequency derivative causes the vector addition to form an arc, although in this case quasi-parabolic rather than circular. The decreased DFT response equals the distance from the origin to the end of the arc. This distance is significantly shorter than that of a coherently detected signal which equals the the length along the arc.

Signals with constant or nearly constant $`\dot{f}`$ are quite common in pulsar astronomy — especially when dealing with time series of very long duration. In such long observations, even the very small spin-downs typical of pulsars can cause a signal to drift across numerous Fourier bins. The Doppler effects of binary pulsar orbits cause similar frequency drifting when the observation time is much shorter than the orbital period.

The “standard” method to correct for a constant frequency derivative is to “stretch” the original time series to compensate for the known or assumed $`\dot{f}`$. This process involves re-sampling the data ensemble $`n_j`$ using a transformation similar to

$$t^{}=t+\frac{2\dot{f}}{f_o}t^2,$$
(32)
where $`t`$ is the time used when sampling the original data, $`\dot{f}`$ is the frequency derivative, and $`f_o`$ is the initial frequency of the signal. Additional details and variations on the theme can be found in Middleditch & Kristian (1984); Anderson et al. (1990); Johnston & Kulkarni (1991) and Wood et al. (1991).

By stretching the data using the appropriate transform and then FFTing the corrected time series, we can recover the fully coherent response. Such techniques have been used with significant success in searches of relatively short time-series (e.g. Camilo et al. 2000). However, this technique runs into significant difficulties when trying large numbers of transformations using long time series where computation of the FFT is non-trivial. Both techniques that we will mention allow full corrections to be made to a signal without requiring multiple FFTs of the full original data set.

#### 4.2.2 Correlation Technique

The correlation technique is the more powerful of the two methods, and uses matched filtering in the Fourier domain to “sweep-up” signal spread over a number of frequency bins into a single bin. In astrophysical applications, we usually have some sort of “pure” signal (like a harmonic from a millisecond pulsar), whose frequency changes as a function of time due to some other process (such as orbital motion or pulsar spin-down). In the Fourier domain, these processes cause the perfect $`\mathrm{sinc}`$-like response of a harmonic to be spread over numerous local Fourier bins — in effect, the $`\mathrm{sinc}`$ response is convolved with a finite impulse response (FIR) filter (where finite in this case refers to a small portion, say $`m`$ bins, of the frequency range analyzed rather than a short period of time). If we can predict the complex form (and phase) of that FIR filter, we can recover the coherent response (i.e. the perfect $`\mathrm{sinc}`$-function) by correlating the appropriate Fourier bins with a “frequency-reversed” and complex-conjugated template that matches the filter.

In mathematical terms, consider a signal with a normalized Fourier response of $`A_{kr_o}`$, where $`kr_o`$ is simply the frequency offset of bin $`k`$ from some reference frequency $`r_o`$, which goes to zero as $`\left|kr_o\right|`$ approaches some number of bins $`m/2`$. For Fourier interpolation as described in §4.1, this response is equal to eqn. 18c without the $`A_o`$ factor (i.e. normalized to an amplitude of one for a coherent response). The complex-valued Fourier response of such a signal at frequency $`r_o`$ can be calculated with the sum

$$A_{r_o}\underset{k=[r_o]m/2}{\overset{[r_o]+m/2}{}}A_kA_{r_ok}^{}.$$
(33)
If $`r_o`$ is initially unknown (i.e. we are searching for a signal with the response shape as defined by the template but at an unknown frequency) we simply compute this summation at a range of frequencies $`r`$.

Calculating eqn. 33 over a range of evenly spaced frequencies is equivalent to correlating the raw FFT amplitudes with the template and is therefore most efficiently computed using short FFTs and the convolution theorem. With FFTs of length $`M`$, such that $`mMN/2`$, we can search a very-long FFT of length $`N/2`$ for any signal whose $`A_{kr_o}`$ we can compute, using overlap-and-save or overlap-and-add techniques (see e.g. Press et al. 1992). Such calculations have advantages over standard time-domain “stretching” techniques in that that they are memory local and can be easily parallelized — important properties when dealing with very-long time-series and modern distributed memory computer architectures.

Moving to our example of a signal with a constant frequency derivative, a single harmonic of the signal has the form

$`n(u)`$ $`=`$ $`a\mathrm{cos}\left[2\pi \left(r_ou+{\displaystyle \frac{\dot{r}}{2}}u^2\right)+\varphi \right]`$ (34a)
$`=`$ $`{\displaystyle \frac{a}{2}}\left[\mathrm{e}^{2\pi \mathrm{i}\left(r_ou+\frac{\dot{r}}{2}u^2\right)}\mathrm{e}^{\mathrm{i}\varphi }+\mathrm{e}^{2\pi \mathrm{i}\left(r_ou+\frac{\dot{r}}{2}u^2\right)}\mathrm{e}^{\mathrm{i}\varphi }\right],`$ (34b)
where we use the same notation as §2.1. Neglecting the second term as in §3.2.1 and Fourier transforming at some “center” or average frequency $`r_c^{}=r+\dot{r}/2`$, we get

$$A_{r_c^{}}=\frac{aN}{2}\mathrm{e}^{\mathrm{i}\varphi }_0^1\mathrm{e}^{\mathrm{i}\pi \left(\dot{r}u^2+2q_ru\right)}du,$$
(35)
where $`q_r=r_cr_c^{}`$, and the real “center” frequency of the signal is $`r_c=r_o+\dot{r}/2`$. This integral can be evaluated in closed form

$$_0^1\mathrm{e}^{\mathrm{i}\pi \left(\dot{r}u^2+2q_ru\right)}du=\frac{1}{\sqrt{2\dot{r}}}\mathrm{e}^{\mathrm{i}\pi \frac{q_r^2}{\dot{r}}}\left\{S\left(Z_r\right)S\left(Y_r\right)+\mathrm{i}\left[C\left(Y_r\right)C\left(Z_r\right)\right]\right\}$$
(36)
where $`Y_r=\sqrt{2/\dot{r}}q_r,`$ $`Z_r=\sqrt{2/\dot{r}}\left(q_r+\dot{r}\right),`$ and $`C(x)`$ and $`S(x)`$ are the Fresnel integrals

$$C(x)=_0^x\mathrm{cos}\left(\frac{\pi }{2}t^2\right)dt,S(x)=_0^x\mathrm{sin}\left(\frac{\pi }{2}t^2\right)dt.$$
(37)
The Fourier transform response then becomes

$$A_{r_c^{}}=\frac{aN}{2\sqrt{2\dot{r}}}\mathrm{e}^{\mathrm{i}\left(\varphi \pi \frac{q_r^2}{\dot{r}}\right)}\left\{S\left(Z_r\right)S\left(Y_r\right)+\mathrm{i}\left[C\left(Y_r\right)C\left(Z_r\right)\right]\right\}.$$
(38)
Using the correlation technique, the coherent response can be recovered by convolving local Fourier amplitudes with the “frequency-reversed” and complex-conjugated template as defined by eqn. 38. This response, at average Fourier frequency $`r_c`$ and Fourier frequency derivative $`\dot{r}`$ can therefore be written as

$$A_{r_c,\dot{r}}=\underset{k=\left[r\right]m/2}{\overset{k=\left[r\right]+m/2}{}}A_k\frac{1}{\sqrt{2\dot{r}}}\left(\mathrm{e}^{\mathrm{i}\pi \frac{q_k^2}{\dot{r}}}\left\{S\left(Z_k\right)S\left(Y_k\right)\mathrm{i}\left[C\left(Y_k\right)C\left(Z_k\right)\right]\right\}\right).$$
(39)
Eqn. 39 takes into account the fact that the signal has been “spread” relatively evenly into the $`\dot{r}`$ closest frequency bins to $`r_c`$, while an additional small amount of signal has “leaked” into bins further away — much like the non-zero wings of the $`\mathrm{sinc}`$-like response to a constant frequency signal. As a rule of thumb, the correct Fourier amplitude will be well-approximated if $`m`$ is chosen such that $`\dot{r}<m2\dot{r}`$.

Large scale searches of pulsations with constant frequency derivatives have been conducted using the correlation technique. A successful example is a search for pulsars in globular cluster NGC 6544 using data taken with the Parkes radio telescope (Ransom et al. 2001). The search was conducted using an FFT of 13865600 points over Fourier $`\dot{r}`$ values from $`100`$ to $`100`$ with a step-size of $`\mathrm{\Delta }\dot{r}=2`$ and included the calculation of amplitudes at half-bin frequency intervals. The search resulted in the detection of the 3.06 ms PSR J1807-2459 in a low-mass binary with orbital period 1.7 h, the second shortest radio pulsar orbital period known. A detailed view of the pulsar’s fundamental harmonic is shown in Fig. 6. The plot was calculated using the correlation technique with spacings of $`\mathrm{\Delta }r=0.0625`$ and $`\mathrm{\Delta }\dot{r}=0.25`$. The generation of such a piece of the $`f`$$`\dot{f}`$ plane takes only a fraction of a second on a rather modest workstation.

#### 4.2.3 Vector Bending

Vector bending is one of the simplest and most straightforward methods to correct a Fourier response that has been smeared over several local frequency bins. As we described in §2.1, the DFT can be thought of as the vector addition of $`N`$ complex numbers. This addition produces a straight line in the complex plane for a coherently detected sinusoid. For a sinusoid with a non-integral or time-varying pulsation frequency, the standard DFT addition produces a curved shape (see Fig. 1). Since the amplitude of the Fourier response is the distance between the origin and the end-point of the vector addition, any curvature in the vector addition implies non-optimal signal detection.

The precise shape of the response curve in the complex plane depends on the mismatch of the signal’s (possibly time-dependent) pulsation frequency and the frequency used in forming the DFT addition (i.e. the closest FFT bin). Regardless of the shape, though, for short enough segments of the curve, the segments differ little from straight lines. We can therefore approximate the shape of the curve as a sum of $`G`$ linear segments, each of which contains $`N/G`$ points from the full-resolution vector addition. In terms of the $`r^{\mathrm{th}}`$ DFT amplitude, we can write this as

$$A_r=\underset{g=0}{\overset{G1}{}}B_{r,g}=\underset{g=0}{\overset{G1}{}}\underset{h=0}{\overset{\frac{N}{G}1}{}}n_j\mathrm{e}^{2\pi \mathrm{i}jr/N},$$
(40)
where $`j=gN/G+h`$. This is equivalent to calculating and then summing $`G`$ independent DFTs (the $`B_{r,g}`$), each of which suffers virtually no loss in sensitivity when the curvature over a segment is small.

Using these vector addition segments or sub-vectors, we can correct for the loss of sensitivity due to curvature by simply straightening the vector addition. If we can predict the true pulsation frequency of a signal as a function of time, we can predict how much curvature will accumulate in each sub-vector and then remove it by rotating the segment appropriately.

For our example of a signal with constant frequency derivative, the instantaneous phase (ignoring the intrinsic phase of the signal) is equal to

$`\mathrm{\Phi }_{true}(u)`$ $`=`$ $`2\pi {\displaystyle r(u)du}`$ (41a)
$`=`$ $`2\pi {\displaystyle \left(r_o+\dot{r}u\right)du}`$ (41b)
$`=`$ $`2\pi r_ou+\pi \dot{r}u^2,`$ (41c)
where $`r_o=f_oT`$ is the initial pulsation frequency of the signal, and $`\dot{r}=\dot{f}T^2`$ is the frequency derivative. The process of taking a DFT removes an instantaneous phase equivalent to $`\mathrm{\Phi }_{DFT}(u)=2\pi ru`$ from the signal (see eqn. 2). So the instantaneous phase error is equal to

$`\mathrm{\Phi }_{error}(u)`$ $`=`$ $`\mathrm{\Phi }_{true}(u)\mathrm{\Phi }_{DFT}(u)`$ (42a)
$`=`$ $`2\pi (r_or)u+\pi \dot{r}u^2.`$ (42b)
Therefore, to correct a particular signal using vector bending, we first calculate the $`B_{r,g}`$ using eqn. 40 for a particular Fourier frequency $`r`$ (such as the frequency of a known pulsar). Now we attempt to un-bend the full Fourier response by summing the $`B_{r,g}`$ after correcting for the phase errors $`\mathrm{\Phi }_{error}(u_g)`$ as defined by eqn. 42a. The corrected response is equal to

$$A_{r,\dot{r}}=\underset{g=0}{\overset{G1}{}}B_{r,g}\mathrm{e}^{\mathrm{i}\mathrm{\Phi }_{error}(u_g)}.$$
(43)
A choice of $`G10^3`$ will essentially eliminate the loss of response for reasonable frequency offsets and frequency derivatives (i.e. less than a few 10s of Fourier bins).

While impractical for large-scale searches due to the fact that the $`B_{r,g}`$ must be recomputed every few $`r`$, vector bending offers significant computational advantages in certain situations. In particular, X-ray observations often consist of short ($`<1`$ hour) “on-source” segments separated by hours, days, or even weeks of “off-source” time (see also §3.4). An FFT of the entire time series might be prohibitively expensive. However, if we can determine an “initial guess” frequency (e.g. by FFTing one segment of the observation or from an ephemeris), we can quickly calculate the $`B_{k,g}`$ at this frequency from the “on-source” intervals alone. We can reconstruct the $`f`$$`\dot{f}`$ plane around our frequency of interest without ever creating the full “filled-in” or padded time series, let alone calculating a potentially huge FFT. Such techniques have allowed us to perform frequency analysis of very-long stretches of data from ROSAT observations of PSR B0540-59 (Eikenberry et al. 1998). Similarly, Fig. 1 shows the use of the method for a 2.4 day observation of the Crab pulsar.

## 5 Signal Property Estimation

Besides correcting for losses of sensitivity, Fourier interpolation and the other response correcting techniques mentioned in §4.2 allow us to determine other useful properties of a detected signal. Using detailed amplitude and phase information from a signal’s Fourier harmonics, we can estimate properties such as the statistical significance of the signal, the location and duration of the signal in the time series (the “centroid” and “purity” respectively), the precise pulsation frequency, and estimates of the measurement errors for Fourier power and phase.

The first step when estimating signal properties in the Fourier domain is to isolate the true peak of the Fourier response in power. This is easily accomplished by using the matched filtering techniques to generate an oversampled grid of amplitudes near and around the signal candidate (see Figs. 7 and 6). Simple optimization algorithms such as the downhill simplex method can then be used to refine the peak location (e.g. Press et al. 1992). Once the peak has been located, estimates of the first and second derivatives of power and phase with respect to Fourier frequency, obtained using Fourier interpolation, can be used to calculate various useful signal properties (Middleditch et al. 1993).

### 5.1 Power, Phase and Signal Amplitude

When the peak of the Fourier response has been located as a function of Fourier frequency and the other search parameters, the measured power is defined as

$$P_{meas}=\frac{\left|A_{r,\mathrm{}}\right|^2}{P_{norm}},$$
(44)
where $`P_{norm}`$ is the expected noise power and is usually described by one of $`Nd_j^2`$, $`P_{local}`$, or $`n_{ph}`$ as discussed in §3.1. Groth (1975, see §3.3) showed that since the measured power is a random variable due to the presence of noise, its variance is $`2P_{signal}+1`$, where $`P_{signal}`$ is the power caused by the signal. Since we do not *a priori* know the true signal power, a good estimate for the variance of the measured power is simply $`2P_{meas}1`$ since $`P_{meas}=P_{signal}+1`$.

Using $`P_{meas}`$, as well as the knowledge that a sinusoid of amplitude $`a`$ in a noisy time series produces a power with an expectation value of $`P_{meas}=a^2N^2/\left(4P_{norm}\right)+1`$ (see §3.2.1), we can estimate the signal amplitude as

$$a=\frac{2}{N}\sqrt{P_{norm}\left(P_{meas}1\right)}.$$
(45)
For binned data containing a signal with Fourier frequency $`r`$, the measured power should be multiplied by $`1/\mathrm{sinc}^2(\pi r/N)`$ to correct for the loss in sensitivity due to binning (see §3.4.1). Vaughan et al. (1994) provide detailed instructions on how to estimate upper limits on pulsation amplitudes as well as estimates of the overall sensitivity of a search.

The statistical significance of a signal is also determined by $`P_{meas}`$. The probability that noise caused a particular power to equal or exceed $`P_{meas}`$ is given by $`\mathrm{e}^{P_{meas}}`$ (eqn. 10 with $`P^{}=P_{meas}`$). But for a search over $`N_{IFS}`$ independent Fourier powers, the probability that at least one of the noise powers exceeds $`P_{meas}`$ is given by

$$Prob(P_{noise}P_{meas})=1\left(1\mathrm{e}^{P_{meas}}\right)^{N_{IFS}}.$$
(46)
Vaughan et al. (1994) show how to use this information to set detection thresholds that minimize the number of spurious candidate signals and give high confidence that signals with powers above the thresholds are real.

Using the real and imaginary parts of the peak Fourier response, we can also calculate the phase of the sinusoidal signal as

$$\varphi _{meas}=\mathrm{arctan}\left[\frac{\mathrm{Im}\left(A_{r,\mathrm{}}\right)}{\mathrm{Re}\left(A_{r,\mathrm{}}\right)}\right]$$
(47)
radians. Using similar arguments as for the measured power, the variance of the measured phase is approximately $`1/(2P_{meas}1)`$ radians.

### 5.2 Signal Location and Duration in Time

Astronomical observations of pulsations effectively consist of a *window* of on-source time where pulsations are present, and the rest of the time when they are not. For most of this paper we have assumed that a signal is present throughout the observation as evidenced by the limits of integration for eqn. 2, which in time-normalized units go from $`0u1`$, or equivalently, from $`0tT`$. In effect, pulsations such as that defined in eqn. 14a are multiplied by a square window function defined as 1 during the observation and 0 at all other times. This window function is simply a property of the DFT and is due to the finite duration of our observation.

It is possible, though, and for various reasons often likely, that a signal we are observing turns on and off or varies in intensity during an observation. The behavior of the signal itself effectively defines a new window function, $`W(u)`$. By measuring the moments of this window function with respect to time, we can determine approximately where in our data a signal is located and for how long.

The approximate location of a signal in a time-normalized time series is described by its centroid, $`\widehat{C}=u=t/T`$ which is proportional to the first moment of the window function with respect to time. More specifically,

$$\widehat{C}=\frac{_{\mathrm{}}^{\mathrm{}}uW(u)du}{_{\mathrm{}}^{\mathrm{}}W(u)du}.$$
(48)
Middleditch & Cordova (1982) wrote this in terms of the measured Fourier response as

$$\widehat{C}=\frac{1}{2\pi }\frac{\varphi (r_o)}{r}\pm \frac{1}{\sqrt{24P(r_o)}},$$
(49)
where $`P(r_o)`$ and $`\varphi (r_o)`$ is the phase measured at the peak of the Fourier response (i.e. $`\varphi _{meas}`$) and $`r`$ is the Fourier frequency (see Appendix D for a derivation). Signals present throughout an observation have $`\widehat{C}=1/2`$ while those present in only the first or second halves of the observation have $`\widehat{C}=1/4`$ or $`\widehat{C}=3/4`$ respectively.

The second moment of the window function with respect to time is related to the moment of inertia of a function and can therefore be used to estimate the root-mean-squared (rms) dispersion of the pulsations in time about the centroid. Using this information, Middleditch & Cordova (1982) defined a parameter called the “purity” (and symbolized by $`\alpha `$) as

$$\alpha =\frac{1}{\pi }\sqrt{\frac{3}{2P(r_o)}\frac{^2P(r_o)}{r^2}}\pm \frac{1}{\alpha \sqrt{10P(r_o)}},$$
(50)
where $`P(r_o)`$ is the measured power at the the peak of the Fourier response (i.e. $`P_{meas}`$, see Appendix E for a derivation). The scaling in eqn. 50 is chosen such that the rms dispersion of the signal about the centroid for a window function $`W(u)`$, is equivalent to that of a rectangular window function of duration $`\alpha `$ (in units of the time series length) centered on the centroid. A signal present throughout the data would have $`\alpha =1`$, while one present in only half the data (in a continuous section) would have $`\alpha =1/2`$. Signals present only at the start and end of an observation but absent in the middle have $`\alpha >1`$. Purity can also help to identify sidelobes caused by a periodic modulation of a signal as these Fourier amplitudes have $`\alpha =\sqrt{3}`$.

Since the location and duration of a signal in an time series affects the Fourier response, it is important to understand how eqn. 18c changes if a signal is present during only part of an observation. In Appendix F, we show that when close to the peak of a signal’s Fourier response,

$$A_rA_o\mathrm{e}^{2\pi \mathrm{i}\widehat{C}(rr_o)}\mathrm{sinc}\left[\pi \alpha (rr_o)\right]$$
(51)
where $`A_o=Na/2\mathrm{e}^{\mathrm{i}\varphi _o}`$, $`\varphi _o`$ is the intrinsic phase of the signal, and $`r_o`$ is the true signal frequency (in FFT bins). This equation demonstrates that for centroids different from $`1/2`$, the phase shift between consecutive FFT bins differs from the $`\pi `$ radians shown in eqn. 18c. Similarly, for purity values different from 1, neighboring FFT bins show more or less correlation with each other (i.e. the central peak of the sinc function changes its width).

### 5.3 The Pulsation Frequency and Frequency Derivative

The true pulsation frequency of the signal is located at the point where $`P/r=0`$. Furthermore, given the response in eqn. 51, we can show (Appendix B) that the uncertainty in this measurement (in Fourier bins) is given by

$$\sigma _r=\frac{3}{\pi \alpha \sqrt{6P_{meas}}}.$$
(52)
This uncertainty is considerably smaller than the often-quoted “frequency error” for the FFT of one bin width, which is simply the frequency resolution returned by the FFT algorithm.

If the correlation method is used to isolate a peak in the $`f`$$`\dot{f}`$ plane as shown in Figs 7 and 6, we can calculate the uncertainty in the measured $`\dot{r}`$ value by using similar arguments and methods as for the frequency uncertainty (see Appendix C for a derivation). The uncertainty in the $`\dot{r}`$ (in Fourier bins) is approximately

$$\sigma _{\dot{r}}=\frac{3\sqrt{10}}{\pi \alpha ^2\sqrt{P_{meas}}}.$$
(53)
## 6 Conclusions

In this paper we have described techniques that allow sophisticated and fully coherent Fourier analysis of very long time series. Most of these techniques use the wealth of information provided by the Fourier phases — information discarded during “standard” analyses based on raw power spectra.

Significant gains in sensitivity and efficiency are possible when using Fourier phase information during the search for periodic signals (using the Fourier domain matched filtering techniques described in §4.2) and when characterizing signals that are known to be present in the data (using the parameters described in §5). The methods of Fourier domain matched filtering allow efficient, memory local, and inherently parallel analysis of extremely long time series with only modest computational resources. Billion point out-of-core FFTs followed by fully coherent matched filtering pulsation searches are possible on “standard” workstations. More traditional time domain based techniques (such as acceleration searches performed by stretching or compressing the time series followed by large in-core FFTs) on similarly sized time series require specialized high-performance computing resources, assuming they can be performed at all.

As astronomical instruments become more sophisticated and specialized, time series of ever increasing duration and time resolution will appear. The Fourier domain techniques described in this paper should prove to be essential tools in their analysis.

*Acknowledgments* We would like to thank G. Fazio for supporting our research and encouraging us in this work. Additional thanks go to J. Grindlay, B. Schmidt, F. Seward, R. Narayan, V. Kaspi, and L. Greenhill for their encouragement and comments. Many of the computations for this paper were performed on equipment purchased with NSF grant PHY 9507695. S.M.R. acknowledges the support of a Tomlinson Fellowship awarded by McGill University. S.S.E. is supported in part by an NSF CAREER Grant.

## Appendix A Derivation of Fourier Interpolation

Following the derivation found in Middleditch et al. (1993), we begin with the definition of the $`k^{\mathrm{th}}`$ DFT element

$$A_k=\underset{j=0}{\overset{N1}{}}n_j\mathrm{e}^{2\pi \mathrm{i}jk/N},$$
(A1)
which we then rewrite by substituting the inverse DFT for the $`n_j`$

$`A_k`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{N1}{}}}\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{l=0}{\overset{N1}{}}}A_l\mathrm{e}^{2\pi \mathrm{i}jl/N}\right)\mathrm{e}^{2\pi \mathrm{i}jk/N}`$ (A2a)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{l=0}{\overset{N1}{}}}A_l{\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{e}^{2\pi \mathrm{i}j(kl)/N}.`$ (A2b)
The last summation can be computed exactly using the identity

$$\underset{j=0}{\overset{N1}{}}\mathrm{e}^{\mathrm{i}\alpha j}=\mathrm{e}^{\frac{\mathrm{i}\alpha }{2}(N1)}\frac{\mathrm{sin}\left(\frac{N\alpha }{2}\right)}{\mathrm{sin}\left(\frac{\alpha }{2}\right)},$$
(A3)
such that when $`N1`$ we have

$`{\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{e}^{2\pi \mathrm{i}j(kl)/N}`$ $`=`$ $`\mathrm{e}^{\mathrm{i}\pi (kl)(1\frac{1}{N})}{\displaystyle \frac{\mathrm{sin}\left[\pi (kl)\right]}{\mathrm{sin}\left[\pi (kl)/N\right]}}`$ (A4a)
$``$ $`\mathrm{e}^{\mathrm{i}\pi (kl)}{\displaystyle \frac{\mathrm{sin}\left[\pi (kl)\right]}{\pi (kl)/N}}`$ (A4b)
$``$ $`N\mathrm{e}^{\mathrm{i}\pi (kl)}\mathrm{sinc}\left[\pi (kl)\right].`$ (A4c)
Substituting this expression into eqn. A2b and changing the integer frequency $`k`$ into a continuous real-valued frequency $`r`$, we arrive at eqn. 29

$$A_r=\underset{l=0}{\overset{N1}{}}A_l\mathrm{e}^{\mathrm{i}\pi (rl)}\mathrm{sinc}\left[\pi (rl)\right].$$
(A5)
## Appendix B Derivation of Frequency Uncertainty

For a given Fourier frequency offset $`\mathrm{\Delta }_r=rr_o`$, where $`r_o`$ is the Fourier frequency of the signal, the magnitude of the Fourier response and the power go as (see Appendix F)

$`\left|A(r)\right|`$ $`=`$ $`\left|A_o\right|\mathrm{sinc}\left(\pi \alpha \mathrm{\Delta }_r\right)`$ (B1a)
$`P(r)`$ $`=`$ $`P_{meas}\mathrm{sinc}^2\left(\pi \alpha \mathrm{\Delta }_r\right)`$ (B1b)
where $`\alpha `$ is the signal purity. We can expand the $`\mathrm{sinc}`$ function in order to approximate the expression for the power near the peak of the response as

$$P(r)=P_{meas}\left[\frac{\mathrm{sin}(\pi \alpha \mathrm{\Delta }_r)}{\pi \alpha \mathrm{\Delta }_r}\right]^2P_{meas}\left[1\frac{\left(\pi \alpha \mathrm{\Delta }_r\right)^2}{3}\right].$$
(B2)
Taking the derivative of power with respect to $`r`$ and solving for $`\mathrm{\Delta }_r`$ we obtain

$$\mathrm{\Delta }_r=\frac{3}{2\pi ^2\alpha ^2P_{meas}}\frac{P}{r}.$$
(B3)
As expected, when the Fourier frequency equals the true frequency of the pulsations (i.e. $`\mathrm{\Delta }_r=0`$), the Fourier response peaks and $`P(r)/r=0`$.

In order to estimate the uncertainty in $`\mathrm{\Delta }_r`$, we apply standard propagation of errors to arrive at

$$\sigma _r=\frac{3}{2\pi ^2\alpha ^2P_{meas}}\sigma _{P^{}(r)},$$
(B4)
where we have replaced $`P(r)/r`$ with $`P^{}(r)`$ to simplify the notation. The derivative of the power at the true frequency can be approximated using finite differences as

$$P^{}(r)\frac{P\left(r_o+\mathrm{\Delta }_r\right)P\left(r_o\mathrm{\Delta }_r\right)}{2\mathrm{\Delta }_r}=\frac{P^+P^{}}{2\mathrm{\Delta }_r},$$
(B5)
where we have simply renamed $`P\left(r_o+\mathrm{\Delta }_r\right)`$ and $`P\left(r_o\mathrm{\Delta }_r\right)`$. The uncertainty in $`P^{}(r)`$ can also be approximated using finite differences and error propagation. Since $`P^+`$ and $`P^{}`$ are highly correlated when $`\mathrm{\Delta }_r1`$, their uncertainties are also correlated giving

$$\sigma _{P^{}(r)}\frac{1}{2\mathrm{\Delta }_r}\left(\sigma _{P^+}+\sigma _P^{}\right)\frac{\sigma _{P^+}}{\mathrm{\Delta }_r}.$$
(B6)
Now, we turn to the question of the uncertainty in $`P^+`$, closely following Middleditch (1976). The amplitude of the Fourier response at the true signal frequency can be written as

$$\sqrt{P_{meas}}=\underset{j=0}{\overset{N1}{}}y_j\mathrm{cos}(\varphi _j)$$
(B7)
where $`y_j=n_j/\sqrt{P_{norm}}`$ are the points in our time series as defined in eqn. 1, but scaled using the appropriate $`P_{norm}`$ such that the measured power, $`P_{meas}`$, is properly normalized (see §3.1). Similarly, the $`\varphi _j`$ represent the pulsation phases at times $`u=j/N`$, but rotated by the measured Fourier phase, $`\varphi _{meas}`$, such that the result of the vector addition lies along the real axis in the complex plane (i.e. the final complex phase is zero). In effect, this transform isolates components of the data that are parallel to the final Fourier response.

At a small frequency offset $`\mathrm{\Delta }_r`$ from the true frequency, we can expand the power in a similar fashion as

$$\sqrt{P^+}=\underset{j=0}{\overset{N1}{}}y_j\mathrm{cos}\left(\varphi _j+\delta _{\varphi _j}\right),$$
(B8)
where $`\delta _{\varphi _j}`$ are the “phase errors” introduced by the frequency offset. The phase errors add curvature to the vector addition and are defined as

$$\delta _{\varphi _j}=2\pi \alpha \mathrm{\Delta }_r\left(\frac{j}{N}\frac{1}{2}\right)=2\pi \alpha \mathrm{\Delta }_r\left(u\frac{1}{2}\right),$$
(B9)
where $`u`$ is the normalized time, $`u=t/T=j/N`$. The $`1/2`$ term in eqn. B9 removes the accumulated phase error over the course of the observation (i.e. $`_0^12\pi \alpha \mathrm{\Delta }_rudu=\pi \alpha \mathrm{\Delta }_r`$) and makes the vector summation of eqn. B7 finish on the real axis. Expanding the cosine in eqn. B8 gives

$$\sqrt{P^+}=\underset{j=0}{\overset{N1}{}}y_j\mathrm{cos}(\varphi _j)\mathrm{cos}\left(\delta _{\varphi _j}\right)\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}(\varphi _j)\mathrm{sin}\left(\delta _{\varphi _j}\right).$$
(B10)
Considering the uncertainties in the separate terms of eqn. B10, since the cosine term is an even function of $`\delta _{\varphi _j}`$, it is symmetric about $`r_o`$, and therefore does not contribute to the uncertainty in the $`P^{}(r)`$ measurement as defined by eqn. B5. For the sine term, given that the $`\mathrm{cos}(\varphi _j)`$ “derotates” the signal onto the real axis by definition, we see that

$$\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}(\varphi _j)=0.$$
(B11)
Furthermore, since the $`\varphi _j`$ and $`\delta _{\varphi _j}`$ are uncorrelated, the average value of this term will be zero,

$$\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}(\varphi _j)\mathrm{sin}\left(\delta _{\varphi _j}\right)=0,$$
(B12)
and has no systematic effect on $`P^+`$. However, we can calculate the fluctuations introduced by this term

$$\left(\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}\left(\varphi _j\right)\mathrm{sin}\left(\delta _{\varphi _j}\right)\right)^2=\left(\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}\left(\varphi _j\right)\right)^2\mathrm{sin}^2\left(\delta _{\varphi _j}\right).$$
(B13)
where the cross-terms average to zero since $`\varphi _j`$ and $`\delta _{\varphi _j}`$ are uncorrelated. Due to the normalization of the $`y_j`$, the sum component averages to $`1/2`$. The $`\delta _{\varphi _j}`$ component has an average of

$`\mathrm{sin}^2\left(\delta _{\varphi _j}\right)`$ $`=`$ $`{\displaystyle _0^1}\mathrm{sin}^2\left(\delta _{\varphi _j}\right)du`$ (B14a)
$``$ $`{\displaystyle _0^1}\delta _{\varphi _j}^2du`$ (B14b)
$``$ $`{\displaystyle _0^1}(2\pi \alpha \mathrm{\Delta }_r)^2\left(u^2u+{\displaystyle \frac{1}{4}}\right)du`$ (B14c)
$``$ $`{\displaystyle \frac{(\pi \alpha \mathrm{\Delta }_r)^2}{3}}.`$ (B14d)
Therefore eqn. B13 is equal to $`(\pi \alpha \mathrm{\Delta }_r)^2/6`$, and the variance of $`\sqrt{P^+}`$ will be

$$\sigma _{\sqrt{P^+}}^2=\sqrt{P^+}\frac{(\pi \alpha \mathrm{\Delta }_r)^2}{6}.$$
(B15)
To get the standard deviation we take the square root of this expression. Propagating errors to get the uncertainty in $`P^+`$ adds a factor of $`2\sqrt{P^+}`$ to give

$$\sigma _{P(r_o+\mathrm{\Delta }_r)}=2\pi \alpha \mathrm{\Delta }_r\sqrt{\frac{P_{meas}}{6}},$$
(B16)
where we have used the approximation $`P^+P_{meas}`$. Finally, substituting into equations B6 and B4, we have

$$\sigma _r=\frac{3}{\pi \alpha \sqrt{6P_{meas}}}.$$
(B17)
A much simpler derivation of eqn. B17 is possible if we realize that properly normalized powers times two are distributed according to a $`\chi ^2`$ distribution with 2 degrees of freedom (see §3.1). For such a distribution, a $`1\sigma `$ error corresponds to a change in the measured $`\chi ^2`$ of $`1/2`$. In the case of powers, the $`\pm 1\sigma `$ errors can therefore be found by starting with the expansion of the power around $`P(r_o)=P_{meas}`$ as given in eqn. B2,

$$P_{meas}\left[1\frac{\left(\pi \alpha \mathrm{\Delta }_r\right)^2}{3}\right]=P_{meas}\frac{1}{2},$$
(B18)
and then solving for the $`\mathrm{\Delta }_r`$ that corresponds to $`P_{meas}1/2`$. This gives

$$\mathrm{\Delta }_r=\sigma _r=\frac{3}{\pi \alpha \sqrt{6P_{meas}}}.$$
(B19)
## Appendix C Derivation of $`\dot{f}`$ Uncertainty

If a peak in the $`f`$$`\dot{f}`$ plane has been isolated using techniques similar to those shown in §4.2.2, we can calculate the error in the measurement of the true Fourier frequency derivative $`\dot{r_o}=\dot{f_o}T^2`$ in a manner similar to that for the frequency uncertainty as described in Appendix B. Signals with non-zero frequency derivatives have Fourier peaks which are located off the $`\dot{r}=0`$ line in the $`f`$$`\dot{f}`$ plane, but the *shapes* of those peaks are independent of $`\dot{r}`$ and in fact depend only on the window function of the signal (see Appendix F). The shape of the response in power as a function of $`\mathrm{\Delta }_{\dot{r}}=\dot{r}\dot{r_o}`$ at the “correct” Fourier frequency $`r=r_o`$ is described by eqn. 38, and can be written as

$$A(r_o,\dot{r})=A_o\mathrm{e}^{\mathrm{i}\pi \frac{\mathrm{\Delta }_{\dot{r}}}{4}}\sqrt{\frac{2}{\alpha ^2\mathrm{\Delta }_{\dot{r}}}}\left[S\left(\sqrt{\frac{\alpha ^2\mathrm{\Delta }_{\dot{r}}}{2}}\right)\mathrm{i}C\left(\sqrt{\frac{\alpha ^2\mathrm{\Delta }_{\dot{r}}}{2}}\right)\right],$$
(C1)
when $`q_r`$ is defined as

$$q_r=r_cr_c^{}=\left(r_o+\frac{\dot{r_o}}{2}\right)\left(r_o+\frac{\dot{r_o}+\alpha ^2\mathrm{\Delta }_{\dot{r}}}{2}\right)=\frac{\alpha \mathrm{\Delta }_{\dot{r}}}{2}.$$
(C2)
This definition of $`q_r`$ keeps the magnitude of the Fourier response symmetric about $`q_r`$ no matter what the value of $`\mathrm{\Delta }_{\dot{r}}`$. The power as a function of $`\mathrm{\Delta }_{\dot{r}}`$ is therefore

$$P(r_o,\dot{r})=P_o\frac{2}{\alpha ^2\mathrm{\Delta }_{\dot{r}}}\left\{\left[S\left(\sqrt{\frac{\alpha ^2\mathrm{\Delta }_{\dot{r}}}{2}}\right)\right]^2+\left[C\left(\sqrt{\frac{\alpha ^2\mathrm{\Delta }_{\dot{r}}}{2}}\right)\right]^2\right\}.$$
(C3)
If we expand the Fresnel integrals about $`0`$ as

$$C(x)x\frac{\pi ^2}{40}x^5\mathrm{},\text{and}S(x)\frac{\pi }{6}x^3\mathrm{},$$
(C4)
and then substitute, the power becomes

$$P(r_o,\dot{r})P_o\left[1\frac{(\pi \alpha ^2\mathrm{\Delta }_{\dot{r}})^2}{180}\right].$$
(C5)
Taking the derivative of $`P(r_o,\dot{r})`$ with respect to $`\dot{r}`$ and solving for $`\mathrm{\Delta }_{\dot{r}}`$ gives

$$\mathrm{\Delta }_{\dot{r}}=\frac{90}{\pi ^2\alpha ^4P_o}\frac{P(r_o,\dot{r})}{\dot{r}}.$$
(C6)
From propagation of errors,

$$\sigma _{\dot{r}}=\frac{90}{\pi ^2\alpha ^4P_o}\sigma _{P^{}(r_o,\dot{r})},$$
(C7)
where $`P^{}(r_o,\dot{r})=P(r_o,\dot{r})/\dot{r}`$. Similarly, after using a finite difference estimate of the power derivative following Appendix B, we find

$$\sigma _{P^{}(r_o,\dot{r})}=\frac{\sigma _{P(r_o,\mathrm{\Delta }_{\dot{r}})}}{\mathrm{\Delta }_{\dot{r}}},$$
(C8)
where $`P(r_o,\mathrm{\Delta }_{\dot{r}})`$ represents the power as measured at $`\dot{r}=\dot{r_o}+\mathrm{\Delta }_{\dot{r}}`$.

Closely following Appendix B, we represent $`\sqrt{P(r_o,\mathrm{\Delta }_{\dot{r}})}`$ as a sum of the parallel components of the properly normalized time series points, $`y_j`$, as

$$\sqrt{P(r_o,\mathrm{\Delta }_{\dot{r}})}=\underset{j=0}{\overset{N1}{}}y_j\mathrm{cos}\left(\varphi _j+\delta _{\varphi _j}\right).$$
(C9)
The $`\delta _{\varphi _j}`$ are the “phase errors” introduced when $`\mathrm{\Delta }_{\dot{r}}0`$, and are defined as

$$\delta _{\varphi _j}=\pi \alpha ^2\mathrm{\Delta }_{\dot{r}}\left(u^2u+\frac{1}{6}\right),$$
(C10)
where the $`u`$ term comes from keeping the response symmetric about $`r_o`$ (i.e. $`q_r=\alpha ^2\mathrm{\Delta }_{\dot{r}}/2`$) and the $`1/6`$ removes the accumulated phase error over the course of the observation (i.e. $`_0^1\pi \alpha ^2\mathrm{\Delta }_{\dot{r}}\left(u^2u\right)du=\pi \alpha ^2\mathrm{\Delta }_{\dot{r}}/6`$) and makes the vector summation of eqn. C9 finish on the real axis.

Expanding the cosine term of eqn. C9 gives

$$\sqrt{P(r_o,\mathrm{\Delta }_{\dot{r}})}=\underset{j=0}{\overset{N1}{}}y_j\mathrm{cos}\left(\varphi _j\right)\mathrm{cos}\left(\delta _{\varphi _j}\right)\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}\left(\varphi _j\right)\mathrm{sin}\left(\delta _{\varphi _j}\right).$$
(C11)
The first term shortens both $`P(r_o,\mathrm{\Delta }_{\dot{r}})`$ and $`P(r_o,\mathrm{\Delta }_{\dot{r}})`$ by the same amount and does not affect the derivative of power. For the sine term, since $`\varphi _j`$ and $`\delta _{\varphi _j}`$ are uncorrelated, the average value is zero (see Appendix B), but its fluctuations are important. To calculate the fluctuations, we square the terms and get

$$\left(\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}\left(\varphi _j\right)\mathrm{sin}\left(\delta _{\varphi _j}\right)\right)^2=\left(\underset{j=0}{\overset{N1}{}}y_j\mathrm{sin}\left(\varphi _j\right)\right)^2\mathrm{sin}^2\left(\delta _{\varphi _j}\right).$$
(C12)
Due to the normalization of the $`y_j`$, $`\left(_{j=0}^{N1}y_j\mathrm{sin}\left(\varphi _j\right)\right)^2`$ averages to $`1/2`$ as before, and we can directly calculate $`\mathrm{sin}^2\left(\delta _{\varphi _j}\right)`$ as

$`\mathrm{sin}^2\left(\delta _{\varphi _j}\right)`$ $`=`$ $`{\displaystyle _0^1}\mathrm{sin}^2\left(\delta _{\varphi _j}\right)du`$ (C13a)
$``$ $`{\displaystyle _0^1}\delta _{\varphi _j}^2du`$ (C13b)
$`=`$ $`\pi ^2\alpha ^4\mathrm{\Delta }_{\dot{r}}^2{\displaystyle _0^1}\left(u^2u+{\displaystyle \frac{1}{6}}\right)^2du`$ (C13c)
$`=`$ $`{\displaystyle \frac{\pi ^2\alpha ^4\mathrm{\Delta }_{\dot{r}}^2}{180}}.`$ (C13d)
The fluctuations from the sine term are therefore $`(\pi \alpha ^2\mathrm{\Delta }_{\dot{r}})^2/360`$, and since squaring eqn. C9 doubles the errors, the standard deviation of $`P(r_o,\mathrm{\Delta }_{\dot{r}})`$ is

$$\sigma _{P(r_o,\mathrm{\Delta }_{\dot{r}})}=\sqrt{P(r_o,\mathrm{\Delta }_{\dot{r}})}\frac{2\pi \alpha ^2\mathrm{\Delta }_{\dot{r}}}{\sqrt{360}}\sqrt{P_{meas}}\frac{\pi \alpha ^2\mathrm{\Delta }_{\dot{r}}}{\sqrt{90}}.$$
(C14)
Substituting into equations C8 and then C7 as in Appendix B gives us the uncertainty in the frequency derivative

$$\sigma _{\dot{r}}=\frac{3\sqrt{10}}{\pi \alpha ^2\sqrt{P_{meas}}}.$$
(C15)
## Appendix D Derivation of Centroid

The centroid is a measure of the approximate location of a signal in a time series as estimated by the first moment of the signal with respect to time (see §5.2). We can think of a sinusoidal signal in our data as being always present but modulated in intensity by some window function $`W(u)`$, where $`u=t/T`$ is the normalized time and $`T`$ is the length of the observation. A “normal” observation of a pulsar of constant intensity would therefore have $`W(u)=1`$ when $`0<u<1`$ and $`W(u)=0`$ at all other times (i.e. a square window). Our signal is therefore described by

$$s(u)=a\mathrm{cos}\left(2\pi r_ou+\varphi _o\right)W(u),$$
(D1)
where $`a`$ is the amplitude, $`r_o=f_oT`$ is the Fourier frequency, and $`\varphi _o`$ is the phase of the sinusoid at time $`u=0`$.

Since the centroid of a function is proportional to the first moment of the function with respect to time, we can easily calculate the centroid using the Moment Theorem of Fourier transforms. Bracewell (1999) does this and defines the centroid as

$$u=\frac{A^{}(0)}{2\pi \mathrm{i}A(0)},$$
(D2)
where $`A(0)`$ and $`A^{}(0)`$ are the Fourier transform and its first derivative with respect to $`r`$ measured at $`r=0`$.

Eqn. D2 is not directly applicable for our sinusoidal signal since the information about the window function in eqn. D1 has been shifted to the frequency of the sinusoid in accordance with the Modulation Theorem of Fourier transforms. Accordingly, we can apply the Modulation Theorem to eqn. D2 which gives us

$$u_{s(u)}=\frac{A^{}(r_o)}{2\pi \mathrm{i}A(r_o)}.$$
(D3)
Finally, we can write $`A(r)`$ in phasor form as $`A(r)=a(r)\mathrm{e}^{\mathrm{i}\varphi (r)}`$ where $`a(r)`$ and $`\varphi (r)`$ represent the Fourier amplitude and phase as functions of the Fourier frequency $`r`$. The derivative with respect to Fourier frequency can be written

$$A^{}(r)=\frac{a(r)}{r}\mathrm{e}^{\mathrm{i}\varphi (r)}+a(r)\mathrm{i}\frac{\varphi (r)}{r}\mathrm{e}^{\mathrm{i}\varphi (r)}.$$
(D4)
At the frequency of our signal, the amplitude is $`A(r_o)=a(r_o)\mathrm{e}^{\mathrm{i}\varphi (r_o)}`$ and the Fourier response is at its peak, making $`a(r_o)/r=0`$. Therefore

$$A^{}(r_o)=A(r_o)\mathrm{i}\frac{\varphi (r_o)}{r}.$$
(D5)
Substituting eqn. D5 into eqn. D3 we arrive at

$$\widehat{C}=u_{s(u)}=\frac{1}{2\pi }\frac{\varphi (r_o)}{r},$$
(D6)
which is equivalent to eqn. 49.

We can also estimate the uncertainty on the measured value of the centroid. Following Appendix B from eqns. B8$``$B15, we note that the same noise fluctuations introduced when offsetting from the true pulsation frequency by a small amount $`\mathrm{\Delta }_r`$ that effect the Fourier amplitude will also effect the Fourier phases. When $`\mathrm{\Delta }_r`$ is small and the powers are properly normalized, the amplitude fluctuations of eqn. B13 with variance $`\pi ^2\mathrm{\Delta }_r^2/6`$ correspond to an uncertainty in the phase measurement of $`\sigma _{\varphi (r_o+\mathrm{\Delta }_r)}\pi \mathrm{\Delta }_r/\sqrt{6P(r_o)}`$ radians and therefore

$$\sigma _{\widehat{C}}=\frac{1}{2\pi }\frac{\sigma _{\varphi (r_o+\mathrm{\Delta }_r)}}{\mathrm{\Delta }_r}\frac{1}{\sqrt{24P(r_o)}}.$$
(D7)
## Appendix E Derivation of Purity

The “purity” of a signal (see §5.2) is a measure of the rms dispersion of the pulsations in time with respect to the centroid and is directly proportional to the variance in time of the window function of the sinusoid from eqn. D1. The time variance is defined as

$$\left(uu\right)^2=\frac{_{\mathrm{}}^{\mathrm{}}\left(uu\right)^2W(u)du}{_{\mathrm{}}^{\mathrm{}}W(u)du},$$
(E1)
but can be written as

$`\left(uu\right)^2`$ $`=`$ $`u^2u^2`$ (E2a)
$`=`$ $`{\displaystyle \frac{A^{\prime \prime }(0)}{4\pi ^2A(0)}}+{\displaystyle \frac{\left[A^{}(0)\right]^2}{4\pi ^2\left[A(0)\right]^2}},`$ (E2b)
using the Moment Theorem for Fourier transforms (see e.g. Bracewell 1999). Since our signal is sinusoidal (see Appendix D), application of the Modulation Theorem gives

$`\left(uu\right)^2_{s(u)}`$ $`=`$ $`{\displaystyle \frac{A^{\prime \prime }(r_o)}{4\pi ^2A(r_o)}}+{\displaystyle \frac{\left[A^{}(r_o)\right]^2}{4\pi ^2\left[A(r_o)\right]^2}}`$ (E3a)
$`=`$ $`{\displaystyle \frac{A^{\prime \prime }(r_o)}{4\pi ^2A(r_o)}}\widehat{C}^2.`$ (E3b)
Using $`A(r)=a(r)\mathrm{e}^{\mathrm{i}\varphi (r)}`$ and remembering that $`a(r_o)/r=0`$, the second derivative of the Fourier amplitude is

$$A^{\prime \prime }(r_o)=\mathrm{i}A(r_o)\frac{^2\varphi (r_o)}{r^2}A(r_o)\left(\frac{\varphi (r_o)}{r}\right)^2+\frac{A(r_o)}{a(r_o)}\frac{^2a(r_o)}{r^2}.$$
(E4)
From Appendix D, we know that $`\varphi (r_o)/r=2\pi \widehat{C}`$ which makes the first term of eqn. E4 equal to zero since $`^2\varphi (r_o)/r^2=0`$ and the second term equal to $`4\pi ^2A(r_o)\widehat{C}^2`$. The second derivative of power at the peak response can be written as

$$\frac{^2P(r_o)}{r^2}=\frac{^2}{r^2}\left(A^{}(r_o)A(r_o)\right)=2a(r_o)\frac{^2a(r_o)}{r^2},$$
(E5)
making the third term equal to $`\frac{A(r_o)}{2\left(a(r_o)\right)^2}\frac{^2P(r_o)}{r^2}`$. Substituting into eqn. E3a and simplifying gives

$$\left(uu\right)^2_{s(u)}=\frac{1}{8\pi ^2P(r_o)}\frac{^2P(r_o)}{r^2}.$$
(E6)
If we normalize the variance using the value obtained for a signal present throughout the observation (i.e. a square window where $`W(u)=1`$ from $`0u1`$ and zero elsewhere — which we will call a *pure* signal) where

$$\left(uu\right)^2_{pure}=u^2u^2=\frac{1}{3}\left(\frac{1}{2}\right)^2=\frac{1}{12},$$
(E7)
and then taking the square-root, we are left with

$$\alpha =\sqrt{\frac{\left(uu\right)^2_{s(u)}}{\left(uu\right)^2_{pure}}}=\frac{1}{\pi }\sqrt{\frac{3}{2P(r_o)}\frac{^2P(r_o)}{r^2}},$$
(E8)
which is equivalent to eqn. 50.

In order to estimate the uncertainty in the measurement of $`\alpha `$, we note that by squaring eqn. E8 and substituting the finite difference approximation of

$$\frac{^2P(r_o)}{r^2}\frac{P\left(r_o+\mathrm{\Delta }_r\right)+P\left(r_o\mathrm{\Delta }_r\right)2P\left(r_o\right)}{\mathrm{\Delta }_r^2}=\frac{P^++P^{}2P\left(r_o\right)}{\mathrm{\Delta }_r^2},$$
(E9)
where $`\mathrm{\Delta }_r`$ corresponds to a small frequency offset from the measured peak power $`P(r_o)`$, we get

$$\alpha ^2\frac{3}{2\pi ^2P(r_o)}\frac{P^++P^{}2P\left(r_o\right)}{\mathrm{\Delta }_r^2}.$$
(E10)
This means that

$$\sigma _\alpha \frac{3}{4\pi ^2\alpha \mathrm{\Delta }_r^2P(r_o)}\sigma _{\mathrm{\Delta }_P},$$
(E11)
where $`\mathrm{\Delta }_P=P^++P^{}2P\left(r_o\right)`$ and the extra factor of $`2\alpha `$ comes from converting $`\sigma _{\alpha ^2}`$ to $`\sigma _\alpha `$.

We can then expand the Fourier amplitude around the peak of the signal as in eqn. B10, where we see that the amplitude fluctuations from the sine term (which is anti-symmetric about $`r_o`$) will cancel after the addition of $`P(r_o+\mathrm{\Delta }_r)`$ and $`P(r_o\mathrm{\Delta }_r)`$ in the finite difference approximation shown above. Conversely, the fluctuations due to the cosine term (which is symmetric about $`r_o`$) will add. These fluctuations can be computed by taking the variance of the first non-constant term in the expansion of $`\mathrm{cos}\left(\delta _j\right)`$ which is $`\delta _j^2/2=2\pi ^2\mathrm{\Delta }_r^2(u1/2)^2`$. The computation gives

$`\left(\delta _j^2/2\delta _j^2/2\right)^2`$ $`=`$ $`\left(\delta _j^2/2\right)^2\delta _j^2/2^2`$ (E12a)
$`=`$ $`4\pi ^4\mathrm{\Delta }_r^4\left\{{\displaystyle _0^1}\left(u{\displaystyle \frac{1}{2}}\right)^4du\left[{\displaystyle _0^1}\left(u{\displaystyle \frac{1}{2}}\right)^2du\right]^2\right\}`$ (E12b)
$`=`$ $`\pi ^4\mathrm{\Delta }_r^4\left[{\displaystyle \frac{1}{20}}\left({\displaystyle \frac{1}{6}}\right)^2\right]`$ (E12c)
$`=`$ $`{\displaystyle \frac{\pi ^4\mathrm{\Delta }_r^4}{45}},`$ (E12d)
and since the variations in the $`\mathrm{cos}\left(\varphi _j\right)`$ term average to $`1/2`$, just like the $`\mathrm{sin}\left(\varphi _j\right)`$ term did in eqn. B13, the standard deviation of $`\sqrt{P(r_o+\mathrm{\Delta }_r)}`$ is equal to

$$\sigma _{\sqrt{P(r_o+\mathrm{\Delta }_r)}}=\sqrt{P(r_o+\mathrm{\Delta }_r)}\frac{\pi ^2\mathrm{\Delta }_r^2}{\sqrt{90}}\pi ^2\mathrm{\Delta }_r^2\sqrt{\frac{P(r_o)}{90}}.$$
(E13)
Doubling the error when converting from amplitudes to powers and then doubling it again due to the addition of the errors from the power offsets in eqn. E9 gives

$$\sigma _{\mathrm{\Delta }_P}=4\pi ^2\mathrm{\Delta }_r^2\sqrt{\frac{P(r_o)}{90}}.$$
(E14)
Finally, by substituting into eqn. E11, we get

$$\sigma _\alpha =\frac{3}{4\pi ^2\alpha \mathrm{\Delta }_r^2P(r_o)}4\pi ^2\mathrm{\Delta }_r^2\sqrt{\frac{P(r_o)}{90}}=\frac{1}{\alpha \sqrt{10P(r_o)}}.$$
(E15)
## Appendix F Centroid, Purity, and Fourier Response

In order to consider the effects of centroid and purity on the Fourier response to a sinusoidal signal as described by eqn. D1, we initially assume a Fourier response equal to eqn. 18c of

$$A_r=A_o\mathrm{e}^{\mathrm{i}\pi \left(rr_o\right)}\mathrm{sinc}\left[\pi (rr_o)\right]$$
(F1)
where $`A_o=Na/2\mathrm{e}^{\mathrm{i}\varphi _o}`$, $`\varphi _o`$ is the intrinsic phase of the signal, and $`r_o`$ is the true signal frequency (in FFT bins). This response is correct only for signals with a square window function (i.e. $`W(u)=1`$ from $`0u1`$ and zero elsewhere).

From eqn. F1 we see that a change in Fourier frequency of a single Fourier bin causes a change in the measured Fourier phase of $`\pi `$ radians. This phase change is also visible from the centroid equation for a pure signal with $`\widehat{C}=1/2`$, where $`\mathrm{d}\varphi (r)=\pi `$ for every $`\mathrm{d}r=1`$. Rewriting the centroid equation and integrating, we see that the Fourier phase near the peak response goes as

$$\varphi (r)=2\pi r\widehat{C}+c.$$
(F2)
When $`r=r_o`$, $`\varphi =\varphi _o`$, allowing us to solve for the constant of integration, $`c=\varphi _o+2\pi r_o\widehat{C}`$. Substituting into eqn. F2, we see that the phase of the Fourier response is equal to

$$\varphi (r_o)=\varphi _o2\pi \widehat{C}(rr_o).$$
(F3)
Therefore, for signals that have centroids different from $`1/2`$, the phase change across a single Fourier bin is different from the usual $`\pi `$ radians.

The purity, as described in §5 and Appendix E, is the effective duration of a square window which reproduces the measured rms dispersion of the signal in time about the centroid. Since the Fourier response to a square window goes as $`\mathrm{sinc}(\pi fT)`$, where $`fT=r`$, we can see that replacing the window of length $`T`$ with one of effective duration $`\alpha T`$ causes the Fourier response to go as $`\mathrm{sinc}(\pi \alpha fT)`$. This fact is also approximately true for more complicated window functions as long as $`\left|rr_o\right|1`$. The Fourier response to a windowed sinusoid is therefore

$$A_rA_o\mathrm{e}^{2\pi \mathrm{i}\widehat{C}(rr_o)}\mathrm{sinc}\left[\pi \alpha (rr_o)\right].$$
(F4)
Numerical simulations show that this approximation is valid for purity values $`\alpha 1.5`$.

These same “effective duration” arguments also apply to the shape of the response in the $`\dot{f}`$ direction of the $`f`$$`\dot{f}`$ plane. A change in the effective duration of a signal causes a change in the $`\dot{f}`$ response since $`\dot{f}T^2\dot{f}\alpha ^2T^2`$, or equivalently, $`\dot{r}\dot{r}\alpha ^2`$. The results of this change can be calculated by directly substituting $`\dot{r}\alpha ^2`$ for $`\dot{r}`$ in eqn. 38, as was done in Appendix C.