Scipy FFT
=========

.. contents::


In this tutorial, we dig a bit deeper in scipy's 1D FFT capabilities, exploring
rotational broadening in stellar spectra along the way. For a full list of
available functions, see http://docs.scipy.org/doc/scipy/reference/fftpack.html.

The flow of the programming is as follows (see, e.g., Carroll 1933, Dravins 1990,
Reiners 2001, Simon-Diaz 2006):

1. First, we construct a rotational broadening kernel assuming a linear limb-darkened star
2. We generate a synthetic spectral line, and convolve it with the kernel using ``fftconvolve``.
3. We add some noise on the spectrum.

Next, we assume we don't know the projected equatorial velocity, and derive it
from the generated spectrum, using the Fourier method. We compare an oversampled
computed Fourier transform (``fft``) with the analytical one. So the steps are:

4. Compute the ``fft`` of the generated spectrum and compare with the analytical one (of the broadening kernel only)
5. Compare the output vsini with input value.

Getting started
---------------

We need a lot of functions. First, we import the Fourier functions:

.. sourcecode:: python
    
    from scipy.fftpack import fft,fftfreq    
    from scipy.signal import fftconvolve

For the analytical computation of the Fourier transformation, Bessel functions
are required:

.. sourcecode:: python

    from scipy.special import j1

And of course the usual imports. An explicit import of ``pi``, ``sin``, ``cos``
and ``sqrt`` is done to avoid extensive use of ``np.`` which could make the code
less readable.

.. sourcecode:: python

    import pylab as plt    
    import numpy as np
    from numpy import pi,sin,cos,sqrt

To finalize the setup, we need some constants. We define them here explicitly,
but there exists a module `scipy.constants <http://docs.scipy.org/doc/scipy/reference/constants.html>`_, with an
extensive list of physical constants. However, it does not contain specifical astronomical constants such as the solar mass (see Prsa & Harmanec 2011).

.. sourcecode:: python
    
    midwave = 4480. # central wavelength of line
    cc = 2.99792458e+18 # velocity of light in angstrom/s    
    epsilon = 0.6 # linear limbdarkening parameter
    vsini = 75.*1e13 # input vsini in A/s (75 km/s)    
    dlam = 0.01 # spectrum resolution (A)
    delta = midwave*vsini/cc 

We will add some noise to a generate spectrum. To be able to reproduce the results,
we make sure that we always add the same noise, regardless of how many times
we run the script:

.. sourcecode:: python

    np.random.seed(1111)

The analytical broadening kernel
--------------------------------

The rotational broadening function *G* is only defined on a certain wavelength interval, but
we'll be lazy and solve that issue later. For now, we make the wavelength interval
broad enough, i.e. +/- 5 angstrom, using the spectral resolution given above. 

.. sourcecode:: python

    lambdas = np.linspace(-5,5,10./dlam+1) # wavelength interval
    y = 1-(lambdas/delta)**2 # transformation of wavelengths
    G = (2*(1-epsilon)*sqrt(y)+pi*epsilon/2.*y)/(pi*delta*(1-epsilon/3))  # the kernel
    print(G)

The ``sqrt`` is of course only defined for positive numbers. For negative numbers,
numpy silently fills in ``nan`` (not a number).


.. raw:: html

   <p class="flip1">Click to Show/Hide: NOTES on ``np.nan``</p> <div class="panel1">

.. note::

    * in some versions of scipy and numpy, ``nan`` is not silently introduced,
      but a warning is printed to the screen. You can avoid these print
      statements via ``np.seterr(all='ignore')``.
    * In the example, we were lazy and let the ``nan`` be introduced in the
      arrays. In general, you should be **very careful** with this behaviour,
      since many statements will silently break if there are nan's around:

      .. sourcecode:: python
        
          x = np.random.normal(size=100)    
          x.std()
          x[10] = np.nan
          np.std(x) # std is not defined when there is a nan!
      
      Some functions have an alternative that can handle nans:

      .. sourcecode:: python

          np.nanmax(x)

      but not all of them:

      .. sourcecode:: python

          hasattr(np,'nanstd'), hasattr(np,'nansum')
      
    * You might think ``np.nan`` is a good way to **mask elements in an array**,
      e.g. when you want to compute the mean of an array. However, in the
      previous note we've seen that numpy functions will return ``np.nan`` when
      used on an array containing at least on ``np.nan``. Numpy has a seperate
      array class for arrays with masked values, appropriately named **masked arrays**:

      .. sourcecode:: python
        
         x = np.random.normal(size=(4,4))
         x
         mask = x<0
         y = np.ma.array(x,mask=mask)
         y

      Notice that numpy actually uses a specific fill value for the masked values,
      usually ``1e+20``, but you can control that value. Then finally, we can use
      the normal numpy functions again:

      .. sourcecode:: python
        
         x.std(),y.std()
         x.min(),y.min()

    * ``matplotlib`` can handle ``np.nan``:

      .. sourcecode:: python
      
         X,Y = np.mgrid[-5:5:500j,-5:5:500j]
         Z = X*Y
         Z[np.sqrt((X-1)**2+(Y-2)**2)<1.] = np.nan
         plt.imshow(Z);
         plt.colorbar();
         
    .. image:: matplotlibnan.png
      :scale: 50

.. raw:: html

   </div>


We remove all the wavelength 
and kernel elements wherever the kernel is not defined:

.. sourcecode:: python

    keep = -np.isnan(G) # returns boolean array with `False` where there are nans
    lambdas,G = lambdas[keep],G[keep] # crop the arrays not to contain nans
    lambdas.min(),lambdas.max()

The broadening kernel is only defined between -1.12 and 1.12 A. We could have
foreseen this, since

.. sourcecode:: python

    1.12/midwave*(3e5) # approximate velocity in km/s of 1.12 angstrom at 4480A

Finally make a plot of the rotational kernel.

.. sourcecode:: python

    plt.plot(lambdas,G)

.. image:: rotkernel.png
  :scale: 50

The synthetic line
------------------

Suppose a line has an intrinsic Gaussian profile in absence of rotation:

.. sourcecode:: python

    wavelengths = np.arange(midwave-5,midwave+5+dlam,dlam)
    spec_line = 1-0.5*np.exp( -(wavelengths-midwave)**2/(2*0.1**2)) # width of 0.1A

This needs to be convolved with the rotational broadening kernel. Scipy (and
numpy) have a ``convolve`` function that does not use the FFT, but here we choose
to use the FFT version. We construct the new wavelength array for the convolved
spectrum, and make sure the equivalent width has not changed during the
convolution process:

.. sourcecode:: python

    spec_conv = fftconvolve(1-spec_line,G,mode='full')
    N = len(spec_conv) # equals len(spec_line)+len(G)+1 when mode='full'
    wavelengths_conv = np.arange(-N/2,N/2,1)*dlam + midwave
    EW_before = np.trapz(1-spec_line,x=wavelengths) # trapezoidal integration
    EW_after = np.trapz(spec_conv,x=wavelengths_conv)
    spec_conv = 1-spec_conv/EW_after*EW_before

.. raw:: html

   <p class="flip2">Click to Show/Hide: NOTE on numerical integration</p> <div class="panel2">

.. note::

   Scipy contains a package dedicated to integration: ``scipy.integrate``. It
   also contains the ``trapz`` function, implements the Simpson rule, and
   provides ``quad`` and ``dblquad``, that can integrate functions numerically
   between limits (i.e., without you having to specify an ``x`` and ``y``
   array):

    .. sourcecode:: python
        
        >>> from scipy.integrate import quad
        >>> quad(lambda t: np.exp(-t) / t**5,1,np.inf)
        (0.0704542374617277, 2.8164841189102693e-09)
        

   The return value is a tuple, with the first element holding the estimated
   value of the integral and the second element holding an upper bound on the
   error.
    
   The ``integrate`` package also contains an ODE solver ``odeint``. See e.g.
   the `Zombie Apocalypse example <http://www.scipy.org/Cookbook/Zombie_Apocalypse_ODEINT>`_.

.. raw:: html

   </div>

Finally, we simulate an observed profile by adding random Gaussian noise:

.. sourcecode:: python

    spec_obs = spec_conv + np.random.normal(size=len(spec_conv),scale=0.02)

A plot of the unbroadened (black) and broadened line (red/blue) can be made with:

.. sourcecode:: python

    plt.figure(10);    
    plt.plot(wavelengths,spec_line,'k-');    
    plt.plot(wavelengths_conv,spec_conv,'r-');
    plt.plot(wavelengths_conv,spec_obs,'b-');    
    plt.xlim(4478,4482)

.. image:: speclines.png
  :scale: 50

The analytical Fourier transform
--------------------------------

Let's get back to the rotational kernel. The analytical Fourier transform is
easily computed. First, we need to define the domain in frequency space on which
to compute the transform, and then we evaluate it. Note the use of scipy's
Bessel function:

.. sourcecode:: python

    x = np.linspace(0,30,1000)[1:]  # exclude zero    
    g = 2. / (x*(1-epsilon/3.)) * ( (1-epsilon)*j1(x) +  epsilon* (sin(x)/x**2 - cos(x)/x))
    g = g**2 # convert to power

.. sourcecode:: python

    x /= (2*pi*delta)
    

Computing the Fourier transform
-------------------------------

When calculating the FFT with ``fft``, a complex array is returned. You can get
the real and imaginary part with ``y.real`` and ``y.imag``, and the norm and 
phase angle via ``np.abs(y)`` and ``np.angle(y)``. In this case,
we are only interested in the power. Also, we need to increase the frequency
resolution of the Fourier transform, so that we can nicely resolve the zeros.
This is done via zero padding, and there exists and optional keyword ``n`` to
the function ``fft`` to do just that. If you do not specify ``n``, ``n`` equals
the length of the input array, otherwise the input array is zero-padded until
it is of size ``n``. We choose to make our input array 100 times larger. Finally,
the Fourier method is sensitive to noise outside of the spectral line region.
Therefore, we cut the spectrum to have little continuum remaining.

.. sourcecode:: python

    keep = np.abs(midwave-wavelengths_conv)<1.2 # minimize contribution of continuum
    spec_to_transform = (1-spec_obs)[keep] # we need the continuum at zero
    new_n = 100*len(spec_to_transform) # new length for zeropadding    
    spec_fft = np.abs(fft(spec_to_transform,n=new_n))**2 # power of FFT

The domain of the Fourier transform is retrieved via ``fftfreq``, but we exclude
the negative frequencies:
    
.. sourcecode:: python
    
    x_fft = fftfreq(len(spec_fft),d=dlam)    
    keep = x_fft>=0 # only positive frequencies    
    x_fft, spec_fft = x_fft[keep], spec_fft[keep]

Measuring the vsini
-------------------

The measured vsini corresponds to the first zero in the Fourier transform. Due
to discretization and noise, zero will never be reached, so we need to look for
local minima: those are the locations where the first derivative switches sign
from negative to positive:

.. sourcecode:: python
    
    neg_to_pos = (np.diff(spec_fft[:-1])<=0) & (np.diff(spec_fft[1:])>=0)
    minima = x_fft[1:-1][neg_to_pos]

Let's compare that to the minima from the analytical transform:

.. sourcecode:: python
    
    neg_to_pos = (np.diff(g[:-1])<=0) & (np.diff(g[1:])>=0)
    minima_an = x[1:-1][neg_to_pos]

The frequency domain can be converted to velocity (km/s) as follows:

.. sourcecode:: python

    q1 = 0.610 + 0.062*epsilon + 0.027*epsilon**2 + 0.012*epsilon**3 + 0.004*epsilon**4        
    vsini_zeros = cc/midwave*q1/minima/1e13 
    vsini_zeros_an = cc/midwave*q1/minima_an/1e13

The first zero should be the input vsini:

.. sourcecode:: python

    print(vsini_zeros_an[:10])
    print(vsini_zeros[:10])

Finally, we can make a nice plot: We use the s/A scaling on the x-axis, but
set the tickmarks to correspond to km/s.

.. sourcecode:: python
    
    plt.figure();    
    plt.plot(x,g/g.max());    
    plt.plot(x_fft,spec_fft/spec_fft.max());    
    plt.gca().set_yscale('log');    
    plt.xlim(0.0,0.05);    
    plt.ylim(1e-8,1e1);    
    plt.grid();    
    plt.xticks(minima_an[:5]*midwave/q1/cc*1e13,['{0:.2f}'.format(i) for i in vsini_zeros_an[:5]])

.. image:: vsini.png
  :scale: 50