AverageSpectrum.c

Go to the documentation of this file.

/*
*  Copyright (C) 2007 Bernd Machenschalk, Jolien Creighton, Kipp Cannon
*
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*  the Free Software Foundation; either version 2 of the License, or
*  (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with with program; see the file COPYING. If not, write to the
*  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
*  MA  02111-1307  USA
*/

#include <math.h>
#include <string.h>
#include <lal/LALComplex.h>
#include <lal/FrequencySeries.h>
#include <lal/LALStdlib.h>
#include <lal/LALConstants.h>
#include <lal/AVFactories.h>
#include <lal/Sequence.h>
#include <lal/TimeFreqFFT.h>
#include <lal/Units.h>
#include <lal/Window.h>
#include <lal/Date.h>

#include <lal/LALRCSID.h>
NRCSID (AVERAGESPECTRUMC,"$Id: AverageSpectrum.c,v 1.26 2008/09/30 00:52:56 kipp Exp $");




/**
 *
 * Compute a "modified periodogram," i.e., the power spectrum of a windowed
 * time series.
 *
 */
int XLALREAL4ModifiedPeriodogram(
    REAL4FrequencySeries        *periodogram,
    const REAL4TimeSeries       *tseries,
    const REAL4Window           *window,
    const REAL4FFTPlan          *plan
    )
{
  static const char *func = "XLALREAL4ModifiedPeriodogram";
  REAL4Sequence *work;
  REAL4 normfac;
  UINT4 k;
  int result;

  if ( ! periodogram || ! tseries || ! plan )
      XLAL_ERROR( func, XLAL_EFAULT );
  if ( ! periodogram->data || ! tseries->data )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( tseries->deltaT <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( periodogram->data->length != tseries->data->length/2 + 1 )
      XLAL_ERROR( func, XLAL_EBADLEN );
  if ( window )
  {
    if ( ! window->data )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->sumofsquares <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->data->length != tseries->data->length )
      XLAL_ERROR( func, XLAL_EBADLEN );
  }

  /* if the window has been specified, apply it to data */
  if ( window )
  {
    UINT4 j;

    /* make a working copy */
    work = XLALCutREAL4Sequence( tseries->data, 0, tseries->data->length );
    if ( ! work )
      XLAL_ERROR( func, XLAL_EFUNC );

    /* apply windowing to data */
    for ( j = 0; j < tseries->data->length; ++j )
      work->data[j] *= window->data->data[j];
  }
  else
    /* point to original data */
    work = tseries->data;

  /* compute the power spectrum of the (windowed) timeseries */
  /* CHECKME: are DC and Nyquist right? */
  result = XLALREAL4PowerSpectrum( periodogram->data, work, plan );
  /* destroy the workspace if it was created */
  if ( window )
    XLALDestroyREAL4Sequence( work );
  /* check for errors from the PowerSpectrum call */
  if ( result == XLAL_FAILURE )
    XLAL_ERROR( func, XLAL_EFUNC );

  /* normalize power spectrum to give correct units */
  /* CHECKME: is this the right factor? */
  if ( window )
    normfac = tseries->deltaT / window->sumofsquares;
  else
    normfac = tseries->deltaT / tseries->data->length;
  for ( k = 0; k < periodogram->data->length; ++k )
    periodogram->data->data[k] *= normfac;

  /* now set rest of metadata */
  periodogram->epoch  = tseries->epoch;
  periodogram->f0     = tseries->f0; /* FIXME: is this right? */
  periodogram->deltaF = 1.0 / ( tseries->data->length * tseries->deltaT );

  /* compute units */
  if ( ! XLALUnitSquare( &periodogram->sampleUnits, &tseries->sampleUnits ) )
    XLAL_ERROR( func, XLAL_EFUNC );
  if ( ! XLALUnitMultiply( &periodogram->sampleUnits,
                           &periodogram->sampleUnits, &lalSecondUnit ) )
    XLAL_ERROR( func, XLAL_EFUNC );

  return 0;
}

/**
 *
 * Compute a "modified periodogram," i.e., the power spectrum of a windowed
 * time series.
 *
 */
int XLALREAL8ModifiedPeriodogram(
    REAL8FrequencySeries        *periodogram,
    const REAL8TimeSeries       *tseries,
    const REAL8Window           *window,
    const REAL8FFTPlan          *plan
    )
{
  static const char *func = "XLALREAL8ModifiedPeriodogram";
  REAL8Sequence *work;
  REAL8 normfac;
  UINT4 k;
  int result;

  if ( ! periodogram || ! tseries || ! plan )
      XLAL_ERROR( func, XLAL_EFAULT );
  if ( ! periodogram->data || ! tseries->data )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( tseries->deltaT <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( periodogram->data->length != tseries->data->length/2 + 1 )
      XLAL_ERROR( func, XLAL_EBADLEN );
  if ( window )
  {
    if ( ! window->data )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->sumofsquares <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->data->length != tseries->data->length )
      XLAL_ERROR( func, XLAL_EBADLEN );
  }

  /* if the window has been specified, apply it to data */
  if ( window )
  {
    UINT4 j;

    /* make a working copy */
    work = XLALCutREAL8Sequence( tseries->data, 0, tseries->data->length );
    if ( ! work )
      XLAL_ERROR( func, XLAL_EFUNC );

    /* apply windowing to data */
    for ( j = 0; j < tseries->data->length; ++j )
      work->data[j] *= window->data->data[j];
  }
  else
    /* point to original data */
    work = tseries->data;

  /* compute the power spectrum of the (windowed) timeseries */
  /* CHECKME: are DC and Nyquist right? */
  result = XLALREAL8PowerSpectrum( periodogram->data, work, plan );
  /* destroy the workspace if it was created */
  if ( window )
    XLALDestroyREAL8Sequence( work );
  /* check for errors from the PowerSpectrum call */
  if ( result == XLAL_FAILURE )
    XLAL_ERROR( func, XLAL_EFUNC );

  /* normalize power spectrum to give correct units */
  /* CHECKME: is this the right factor? */
  if ( window )
    normfac = tseries->deltaT / window->sumofsquares;
  else
    normfac = tseries->deltaT / tseries->data->length;
  for ( k = 0; k < periodogram->data->length; ++k )
    periodogram->data->data[k] *= normfac;

  /* now set rest of metadata */
  periodogram->epoch  = tseries->epoch;
  periodogram->f0     = tseries->f0; /* FIXME: is this right? */
  periodogram->deltaF = 1.0 / ( tseries->data->length * tseries->deltaT );

  /* compute units */
  if ( ! XLALUnitSquare( &periodogram->sampleUnits, &tseries->sampleUnits ) )
    XLAL_ERROR( func, XLAL_EFUNC );
  if ( ! XLALUnitMultiply( &periodogram->sampleUnits,
                           &periodogram->sampleUnits, &lalSecondUnit ) )
    XLAL_ERROR( func, XLAL_EFUNC );

  return 0;
}


/**
 *
 * Use Welch's method to compute the average power spectrum of a time series.
 *
 * See: Peter D. Welch "The Use of Fast Fourier Transform for the Estimation
 * of Power Spectra: A Method Based on Time Averaging Over Short, Modified
 * Periodograms"  IEEE Transactions on Audio and Electroacoustics,
 * Vol. AU-15, No. 2, June 1967.
 *
 */
int XLALREAL4AverageSpectrumWelch(
    REAL4FrequencySeries        *spectrum,
    const REAL4TimeSeries       *tseries,
    UINT4                        seglen,
    UINT4                        stride,
    const REAL4Window           *window,
    const REAL4FFTPlan          *plan
    )
{
  static const char *func = "XLALREAL4AverageSpectrumWelch";
  REAL4FrequencySeries *work; /* workspace */
  REAL4Sequence sequence; /* working copy of input time series data */
  REAL4TimeSeries tseriescopy; /* working copy of input time series */
  UINT4 numseg;
  UINT4 seg;
  UINT4 k;

  if ( ! spectrum || ! tseries || ! plan )
      XLAL_ERROR( func, XLAL_EFAULT );
  if ( ! spectrum->data || ! tseries->data )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( tseries->deltaT <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );

  /* construct local copy of time series */
  sequence = *tseries->data;
  tseriescopy = *tseries;
  tseriescopy.data = &sequence;
  tseriescopy.data->length = seglen;

  numseg = 1 + (tseries->data->length - seglen)/stride;

  /* consistency check for lengths: make sure that the segments cover the
   * data record completely */
  if ( (numseg - 1)*stride + seglen != tseries->data->length )
    XLAL_ERROR( func, XLAL_EBADLEN );
  if ( spectrum->data->length != seglen/2 + 1 )
    XLAL_ERROR( func, XLAL_EBADLEN );

  /* make sure window, if present, is appropriate */
  if ( window )
  {
    if ( ! window->data )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->sumofsquares <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->data->length != seglen )
      XLAL_ERROR( func, XLAL_EBADLEN );
  }

  /* clear spectrum data */
  memset( spectrum->data->data, 0,
      spectrum->data->length * sizeof( *spectrum->data->data ) );

  /* create frequency series data workspace */
  work = XLALCutREAL4FrequencySeries( spectrum, 0, spectrum->data->length );
  if( ! work )
    XLAL_ERROR( func, XLAL_EFUNC );

  for ( seg = 0; seg < numseg; seg++, tseriescopy.data->data += stride )
  {
    /* compute the modified periodogram; clean up and exit on failure */
    if ( XLALREAL4ModifiedPeriodogram( work, &tseriescopy, window, plan ) == XLAL_FAILURE )
    {
      XLALDestroyREAL4FrequencySeries( work );
      XLAL_ERROR( func, XLAL_EFUNC );
    }

    /* add the periodogram to the running sum */
    for ( k = 0; k < spectrum->data->length; ++k )
      spectrum->data->data[k] += work->data->data[k];
  }

  /* set metadata */
  spectrum->epoch       = work->epoch;
  spectrum->f0          = work->f0;
  spectrum->deltaF      = work->deltaF;
  spectrum->sampleUnits = work->sampleUnits;

  /* divide spectrum data by the number of segments in average */
  for ( k = 0; k < spectrum->data->length; ++k )
    spectrum->data->data[k] /= numseg;

  /* clean up */
  XLALDestroyREAL4FrequencySeries( work );

  return 0;
}

/**
 *
 * Use Welch's method to compute the average power spectrum of a time series.
 *
 * See: Peter D. Welch "The Use of Fast Fourier Transform for the Estimation
 * of Power Spectra: A Method Based on Time Averaging Over Short, Modified
 * Periodograms"  IEEE Transactions on Audio and Electroacoustics,
 * Vol. AU-15, No. 2, June 1967.
 *
 */
int XLALREAL8AverageSpectrumWelch(
    REAL8FrequencySeries        *spectrum,
    const REAL8TimeSeries       *tseries,
    UINT4                        seglen,
    UINT4                        stride,
    const REAL8Window           *window,
    const REAL8FFTPlan          *plan
    )
{
  static const char *func = "XLALREAL8AverageSpectrumWelch";
  REAL8FrequencySeries *work; /* workspace */
  REAL8Sequence sequence; /* working copy of input time series data */
  REAL8TimeSeries tseriescopy; /* working copy of input time series */
  UINT4 numseg;
  UINT4 seg;
  UINT4 k;

  if ( ! spectrum || ! tseries || ! plan )
      XLAL_ERROR( func, XLAL_EFAULT );
  if ( ! spectrum->data || ! tseries->data )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( tseries->deltaT <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );

  /* construct local copy of time series */
  sequence = *tseries->data;
  tseriescopy = *tseries;
  tseriescopy.data = &sequence;
  tseriescopy.data->length = seglen;

  numseg = 1 + (tseries->data->length - seglen)/stride;

  /* consistency check for lengths: make sure that the segments cover the
   * data record completely */
  if ( (numseg - 1)*stride + seglen != tseries->data->length )
    XLAL_ERROR( func, XLAL_EBADLEN );
  if ( spectrum->data->length != seglen/2 + 1 )
    XLAL_ERROR( func, XLAL_EBADLEN );

  /* make sure window, if present, is appropriate */
  if ( window )
  {
    if ( ! window->data )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->sumofsquares <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->data->length != seglen )
      XLAL_ERROR( func, XLAL_EBADLEN );
  }

  /* clear spectrum data */
  memset( spectrum->data->data, 0,
      spectrum->data->length * sizeof( *spectrum->data->data ) );

  /* create frequency series data workspace */
  work = XLALCutREAL8FrequencySeries( spectrum, 0, spectrum->data->length );
  if( ! work )
    XLAL_ERROR( func, XLAL_EFUNC );

  for ( seg = 0; seg < numseg; seg++, tseriescopy.data->data += stride )
  {
    /* compute the modified periodogram; clean up and exit on failure */
    if ( XLALREAL8ModifiedPeriodogram( work, &tseriescopy, window, plan ) == XLAL_FAILURE )
    {
      XLALDestroyREAL8FrequencySeries( work );
      XLAL_ERROR( func, XLAL_EFUNC );
    }

    /* add the periodogram to the running sum */
    for ( k = 0; k < spectrum->data->length; ++k )
      spectrum->data->data[k] += work->data->data[k];
  }

  /* set metadata */
  spectrum->epoch       = work->epoch;
  spectrum->f0          = work->f0;
  spectrum->deltaF      = work->deltaF;
  spectrum->sampleUnits = work->sampleUnits;

  /* divide spectrum data by the number of segments in average */
  for ( k = 0; k < spectrum->data->length; ++k )
    spectrum->data->data[k] /= numseg;

  /* clean up */
  XLALDestroyREAL8FrequencySeries( work );

  return 0;
}


/* 
 *
 * Median Method: use median average rather than mean.
 *
 */

/* compute the median bias */
REAL8 XLALMedianBias( UINT4 nn )
{
  const UINT4 nmax = 1000;
  REAL8 ans = 1;
  UINT4 n = (nn - 1)/2;
  UINT4 i;

  if ( nn >= nmax )
    return LAL_LN2;

  for ( i = 1; i <= n; ++i )
  {
    ans -= 1.0/(2*i);
    ans += 1.0/(2*i + 1);
  }

  return ans;
}

/* cleanup temporary workspace... ignore xlal errors */
static void median_cleanup_REAL4( REAL4FrequencySeries *work, UINT4 n )
{
  int saveErrno = xlalErrno;
  UINT4 i;
  for ( i = 0; i < n; ++i )
    if ( work[i].data )
      XLALDestroyREAL4Vector( work[i].data );
  XLALFree( work );
  xlalErrno = saveErrno;
  return;
}
static void median_cleanup_REAL8( REAL8FrequencySeries *work, UINT4 n )
{
  int saveErrno = xlalErrno;
  UINT4 i;
  for ( i = 0; i < n; ++i )
    if ( work[i].data )
      XLALDestroyREAL8Vector( work[i].data );
  XLALFree( work );
  xlalErrno = saveErrno;
  return;
}

/* comparison for floating point numbers */
static int compare_REAL4( const void *p1, const void *p2 )
{
  REAL4 x1 = *(const REAL4 *)p1;
  REAL4 x2 = *(const REAL4 *)p2;
  return (x1 > x2) - (x1 < x2);
}
static int compare_REAL8( const void *p1, const void *p2 )
{
  REAL8 x1 = *(const REAL8 *)p1;
  REAL8 x2 = *(const REAL8 *)p2;
  return (x1 > x2) - (x1 < x2);
}


/**
 *
 * Median Method: use median average rather than mean.  Note: this will
 * cause a bias if the segments overlap, i.e., if the stride is less than
 * the segment length -- even though the median bias for Gaussian noise
 * is accounted for -- because the segments are not independent and their
 * correlation is non-zero.
 *
 */
int XLALREAL4AverageSpectrumMedian(
    REAL4FrequencySeries        *spectrum,
    const REAL4TimeSeries       *tseries,
    UINT4                        seglen,
    UINT4                        stride,
    const REAL4Window           *window,
    const REAL4FFTPlan          *plan
    )
{
  static const char *func = "XLALREAL4AverageSpectrumMedian";
  REAL4FrequencySeries *work; /* array of frequency series */
  REAL4 *bin; /* array of bin values */
  REAL4 biasfac; /* median bias factor */
  REAL4 normfac; /* normalization factor */
  UINT4 reclen; /* length of entire data record */
  UINT4 numseg;
  UINT4 seg;
  UINT4 k;

  if ( ! spectrum || ! tseries || ! plan )
      XLAL_ERROR( func, XLAL_EFAULT );
  if ( ! spectrum->data || ! tseries->data )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( tseries->deltaT <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );

  reclen = tseries->data->length;
  numseg = 1 + (reclen - seglen)/stride;

  /* consistency check for lengths: make sure that the segments cover the
   * data record completely */
  if ( (numseg - 1)*stride + seglen != reclen )
    XLAL_ERROR( func, XLAL_EBADLEN );
  if ( spectrum->data->length != seglen/2 + 1 )
    XLAL_ERROR( func, XLAL_EBADLEN );

  /* make sure window, if present, is appropriate */
  if ( window )
  {
    if ( ! window->data )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->sumofsquares <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->data->length != seglen )
      XLAL_ERROR( func, XLAL_EBADLEN );
  }

  /* create frequency series data workspaces */
  work = XLALCalloc( numseg, sizeof( *work ) );
  if ( ! work )
    XLAL_ERROR( func, XLAL_ENOMEM );
  for ( seg = 0; seg < numseg; ++seg )
  {
    work[seg].data = XLALCreateREAL4Vector( spectrum->data->length );
    if ( ! work[seg].data )
    {
      median_cleanup_REAL4( work, numseg ); /* cleanup */
      XLAL_ERROR( func, XLAL_EFUNC );
    }
  }

  for ( seg = 0; seg < numseg; ++seg )
  {
    REAL4Vector savevec; /* save the time series data vector */
    int code;

    /* save the time series data vector */
    savevec = *tseries->data;

    /* set the data vector to be appropriate for the even segment */
    tseries->data->length  = seglen;
    tseries->data->data   += seg * stride;

    /* compute the modified periodogram for the even segment */
    code = XLALREAL4ModifiedPeriodogram( work + seg, tseries, window, plan );
    
    /* restore the time series data vector to its original state */
    *tseries->data = savevec;

    /* now check for failure of the XLAL routine */
    if ( code == XLAL_FAILURE )
    {
      median_cleanup_REAL4( work, numseg ); /* cleanup */
      XLAL_ERROR( func, XLAL_EFUNC );
    }
  }

  /* create array to hold a particular frequency bin data */
  bin = XLALMalloc( numseg * sizeof( *bin ) );
  if ( ! bin )
  {
    median_cleanup_REAL4( work, numseg ); /* cleanup */
    XLAL_ERROR( func, XLAL_ENOMEM );
  }

  /* compute median bias factor */
  biasfac = XLALMedianBias( numseg );

  /* normaliztion takes into account bias */
  normfac = 1.0 / biasfac;

  /* now loop over frequency bins and compute the median-mean */
  for ( k = 0; k < spectrum->data->length; ++k )
  {
    /* assign array of even segment values to bin array for this freq bin */
    for ( seg = 0; seg < numseg; ++seg )
      bin[seg] = work[seg].data->data[k];

    /* sort them and find median */
    qsort( bin, numseg, sizeof( *bin ), compare_REAL4 );
    if ( numseg % 2 ) /* odd number of evens */
      spectrum->data->data[k] = bin[numseg/2];
    else /* even number... take average */
      spectrum->data->data[k] = 0.5*(bin[numseg/2-1] + bin[numseg/2]);

    /* remove median bias */
    spectrum->data->data[k] *= normfac;
  }

  /* set metadata */
  spectrum->epoch       = work->epoch;
  spectrum->f0          = work->f0;
  spectrum->deltaF      = work->deltaF;
  spectrum->sampleUnits = work->sampleUnits;

  /* free the workspace data */
  XLALFree( bin );
  median_cleanup_REAL4( work, numseg );

  return 0;
}

/**
 *
 * Median Method: use median average rather than mean.  Note: this will
 * cause a bias if the segments overlap, i.e., if the stride is less than
 * the segment length -- even though the median bias for Gaussian noise
 * is accounted for -- because the segments are not independent and their
 * correlation is non-zero.
 *
 */
int XLALREAL8AverageSpectrumMedian(
    REAL8FrequencySeries        *spectrum,
    const REAL8TimeSeries       *tseries,
    UINT4                        seglen,
    UINT4                        stride,
    const REAL8Window           *window,
    const REAL8FFTPlan          *plan
    )
{
  static const char *func = "XLALREAL8AverageSpectrumMedian";
  REAL8FrequencySeries *work; /* array of frequency series */
  REAL8 *bin; /* array of bin values */
  REAL8 biasfac; /* median bias factor */
  REAL8 normfac; /* normalization factor */
  UINT4 reclen; /* length of entire data record */
  UINT4 numseg;
  UINT4 seg;
  UINT4 k;

  if ( ! spectrum || ! tseries || ! plan )
      XLAL_ERROR( func, XLAL_EFAULT );
  if ( ! spectrum->data || ! tseries->data )
      XLAL_ERROR( func, XLAL_EINVAL );
  if ( tseries->deltaT <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );

  reclen = tseries->data->length;
  numseg = 1 + (reclen - seglen)/stride;

  /* consistency check for lengths: make sure that the segments cover the
   * data record completely */
  if ( (numseg - 1)*stride + seglen != reclen )
    XLAL_ERROR( func, XLAL_EBADLEN );
  if ( spectrum->data->length != seglen/2 + 1 )
    XLAL_ERROR( func, XLAL_EBADLEN );

  /* make sure window, if present, is appropriate */
  if ( window )
  {
    if ( ! window->data )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->sumofsquares <= 0.0 )
      XLAL_ERROR( func, XLAL_EINVAL );
    if ( window->data->length != seglen )
      XLAL_ERROR( func, XLAL_EBADLEN );
  }

  /* create frequency series data workspaces */
  work = XLALCalloc( numseg, sizeof( *work ) );
  if ( ! work )
    XLAL_ERROR( func, XLAL_ENOMEM );
  for ( seg = 0; seg < numseg; ++seg )
  {
    work[seg].data = XLALCreateREAL8Vector( spectrum->data->length );
    if ( ! work[seg].data )
    {
      median_cleanup_REAL8( work, numseg ); /* cleanup */
      XLAL_ERROR( func, XLAL_EFUNC );
    }
  }

  for ( seg = 0; seg < numseg; ++seg )
  {
    REAL8Vector savevec; /* save the time series data vector */
    int code;

    /* save the time series data vector */
    savevec = *tseries->data;

    /* set the data vector to be appropriate for the even segment */
    tseries->data->length  = seglen;
    tseries->data->data   += seg * stride;

    /* compute the modified periodogram for the even segment */
    code = XLALREAL8ModifiedPeriodogram( work + seg, tseries, window, plan );
    
    /* restore the time series data vector to its original state */
    *tseries->data = savevec;

    /* now check for failure of the XLAL routine */
    if ( code == XLAL_FAILURE )
    {
      median_cleanup_REAL8( work, numseg ); /* cleanup */
      XLAL_ERROR( func, XLAL_EFUNC );
    }
  }

  /* create array to hold a particular frequency bin data */
  bin = XLALMalloc( numseg * sizeof( *bin ) );
  if ( ! bin )
  {
    median_cleanup_REAL8( work, numseg ); /* cleanup */
    XLAL_ERROR( func, XLAL_ENOMEM );
  }

  /* compute median bias factor */
  biasfac = XLALMedianBias( numseg );

  /* normaliztion takes into account bias */
  normfac = 1.0 / biasfac;

  /* now loop over frequency bins and compute the median-mean */
  for ( k = 0; k < spectrum->data->length; ++k )
  {
    /* assign array of even segment values to bin array for this freq bin */
    for ( seg = 0; seg < numseg; ++seg )
      bin[seg] = work[seg].data->data[k];

    /* sort them and find median */
    qsort( bin, numseg, sizeof( *bin ), compare_REAL8 );
    if ( numseg % 2 ) /* odd number of evens */
      spectrum->data->data[k] = bin[numseg/2];
    else /* even number... take average */
      spectrum->data->data[k] = 0.5*(bin[numseg/2-1] + bin[numseg/2]);

    /* remove median bias */
    spectrum->data->data[k] *= normfac;
  }

  /* set metadata */
  spectrum->epoch       = work->epoch;
  spectrum->f0          = work->f0;
  spectrum->deltaF      = work->deltaF;
  spectrum->sampleUnits = work->sampleUnits;

  /* free the workspace data */
  XLALFree( bin );
  median_cleanup_REAL8( work, numseg );

  return 0;
}


/* 
 *
 * Median-Mean Method
 *
 */


/* cleanup temporary workspace... ignore xlal errors */
static void median_mean_cleanup_REAL4( REAL4FrequencySeries *even, REAL4FrequencySeries *odd, UINT4 n )
{
  int saveErrno = xlalErrno;
  UINT4 i;
  for ( i = 0; i < n; ++i )
  {
    if ( even[i].data )
      XLALDestroyREAL4Vector( even[i].data );
    if ( odd[i].data )
      XLALDestroyREAL4Vector( odd[i].data );
  }
  XLALFree( even );
  XLALFree( odd );
  xlalErrno = saveErrno;
  return;
}
static void median_mean_cleanup_REAL8( REAL8FrequencySeries *even, REAL8FrequencySeries *odd, UINT4 n )
{
  int saveErrno = xlalErrno;
  UINT4 i;
  for ( i = 0; i < n; ++i )
  {
    if ( even[i].data )
      XLALDestroyREAL8Vector( even[i].data );
    if ( odd[i].data )
      XLALDestroyREAL8Vector( odd[i].data );
  }
  XLALFree( even );