trunk/SHAPES/FST_semi_memo.c

/***************************************************************************
  **************************************************************************
  
                           S2kit 1.0

          A lite version of Spherical Harmonic Transform Kit

   Peter Kostelec, Dan Rockmore
   {geelong,rockmore}@cs.dartmouth.edu
  
   Contact: Peter Kostelec
            geelong@cs.dartmouth.edu
  
   Copyright 2004 Peter Kostelec, Dan Rockmore

   This file is part of S2kit.

   S2kit 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.

   S2kit 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 S2kit; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

   See the accompanying LICENSE file for details.
  
  ************************************************************************
  ************************************************************************/


/********************************************************************

  FST_semi_memo.c - routines to perform convolutions on the
  2-sphere using a combination of semi-naive and naive algorithms.

  ASSUMES THAT ALL PRECOMPUTED-DATA IS IN MEMORY, AND NOT TO BE
  READ FROM THE DISK.

  The primary functions in this package are

  1) FST_semi_memo() - computes the spherical harmonic expansion.
  2) InvFST_semi_memo() - computes the inverse spherical harmonic transform.
  3) FZT_semi_memo() - computes the zonal harmonic transform.
  4) TransMult() - Multiplies harmonic coefficients using Driscoll-Healy
                    result.  Dual of convolution in "time" domain.
  5) Conv2Sphere_semi_memo() - Convolves two functins defined on the 2-sphere,
                          using seminaive transform

  and one utility function:

  1) seanindex(): Given bandwidth bw, seanindex(m,l,bw) will give the
     position of the coefficient f-hat(m,l) in the one-row array

  For descriptions on calling these functions, see the documentation
  preceding each function.

*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "fftw3.h"

#include "makeweights.h"
#include "cospmls.h"
#include "primitive.h"
#include "naive_synthesis.h"
#include "seminaive.h"
#include "FST_semi_memo.h"


/************************************************************************/


/*****************************************************************

  Given bandwidth bw, seanindex(m,l,bw) will give the position of the
  coefficient f-hat(m,l) in the one-row array that Sean stores the spherical
  coefficients. This is needed to help preserve the symmetry that the
  coefficients have: (l = degree, m = order, and abs(m) <= l)

  f-hat(l,-m) = (-1)^m * conjugate( f-hat(l,m) )

  Thanks for your help Mark!

  ******************************************************************/


int seanindex(int m,
              int l,
              int bw)
{     
  int bigL;

  bigL = bw - 1;

  if( m >= 0 )
    return( m * ( bigL + 1 ) - ( ( m * (m - 1) ) /2 ) + ( l - m ) );
  else
    return( ( ( bigL * ( bigL + 3 ) ) /2 ) + 1 +
            ( ( bigL + m ) * ( bigL + m + 1 ) / 2 ) + ( l - abs( m ) ) );
}


/************************************************************************/
/* performs a spherical harmonic transform using the semi-naive
   and naive algorithms

   bw -> bandwidth of problem
   size -> size = 2*bw -> dimension of input array (recall that
           sampling is done at twice the bandwidth)

   The inputs rdata and idata are expected to be pointers to
   size x size arrays. The array rdata contains the real parts
   of the function samples, and idata contains the imaginary
   parts.

   rcoeffs and icoeffs are expected to be pointers to bw x bw arrays,
   and will contain the harmonic coefficients in a "linearized" form.
   The array rcoeffs contains the real parts of the coefficients,
   and icoeffs contains the imaginary parts.

   spharmonic_pml_table should be a (double **) pointer to
   the result of a call to Spharmonic_Pml_Table.  Because this
   table is re-used in the inverse transform, and because for
   timing purposes the computation of the table is not included,
   it is passed in as an argument.  Also, at some point this
   code may be used as par of a series of convolutions, so
   reducing repetitive computation is prioritized.

   spharmonic_pml_table will be an array of (double *) pointers
   the array being of length TableSize(m,bw)

   workspace needs to be a double pointer to an array of size
   (8 * bw^2) + (7 * bw).

   cutoff -> what order to switch from semi-naive to naive
             algorithm.

   dataformat =0 -> samples are complex, =1 -> samples real

 
   Output Ordering of coeffs f(m,l) is
   f(0,0) f(0,1) f(0,2) ... f(0,bw-1)
          f(1,1) f(1,2) ... f(1,bw-1)
          etc.
                 f(bw-2,bw-2), f(bw-2,bw-1)
                               f(bw-1,bw-1)
                               f(-(bw-1),bw-1)
                 f(-(bw-2),bw-2) f(-(bw-2),bw-1)
          etc.
                  f(-2,2) ... f(-2,bw-1)
          f(-1,1) f(-1,2) ... f(-1,bw-1)

    
   This only requires an array of size (bw*bw).  If zero-padding
   is used to make the indexing nice, then you need a an
   (2bw-1) * bw array - but that is not done here.
   Because of the amount of space necessary for doing
   large transforms, it is important not to use any
   more than necessary.

*/

void FST_semi_memo(double *rdata, double *idata,
                   double *rcoeffs, double *icoeffs, 
                   int bw,
                   double **seminaive_naive_table,
                   double *workspace,
                   int dataformat,
                   int cutoff,
                   fftw_plan *dctPlan,
                   fftw_plan *fftPlan,
                   double *weights )
{
  int size, m, i, j;
  int dummy0, dummy1 ;
  double *rres, *ires;
  double *rdataptr, *idataptr;
  double *fltres, *scratchpad;
  double *eval_pts;
  double pow_one;
  double tmpA, tmpSize ;

  size = 2*bw ;
  tmpSize = 1./ ((double ) size);
  tmpA = sqrt( 2. * M_PI ) ;

  rres = workspace;  /* needs (size * size) = (4 * bw^2) */
  ires = rres + (size * size); /* needs (size * size) = (4 * bw^2) */ 
  fltres = ires + (size * size) ; /* needs bw  */
  eval_pts = fltres + bw; /* needs (2*bw)  */
  scratchpad = eval_pts + (2*bw); /* needs (4 * bw)  */
 
  /* total workspace is (8 * bw^2) + (7 * bw) */

  /* do the FFTs along phi */
  fftw_execute_split_dft( *fftPlan,
                          rdata, idata,
                          rres, ires ) ; 
  /*
    normalize

    note that I'm getting the sqrt(2pi) in there at
    this point ... to account for the fact that the spherical
    harmonics are of norm 1: I need to account for
    the fact that the associated Legendres are
    of norm 1
  */
  tmpSize *= tmpA ;
  for( j = 0 ; j < size*size ; j ++ )
    {
      rres[j] *= tmpSize  ;
      ires[j] *= tmpSize  ;
    }
  
  /* point to start of output data buffers */
  rdataptr = rcoeffs;
  idataptr = icoeffs;
  
  for (m=0; m<bw; m++) {
    /*
      fprintf(stderr,"m = %d\n",m);
    */
    
    /*** test to see if before cutoff or after ***/
    if (m < cutoff){
      
      /* do the real part */
      SemiNaiveReduced(rres+(m*size), 
                       bw, 
                       m, 
                       fltres,
                       scratchpad,
                       seminaive_naive_table[m],
                       weights,
                       dctPlan);
      
      /* now load real part of coefficients into output space */  
      memcpy(rdataptr, fltres, sizeof(double) * (bw - m));
      
      rdataptr += bw-m;
      
      /* do imaginary part */
      SemiNaiveReduced(ires+(m*size),
                       bw, 
                       m, 
                       fltres,
                       scratchpad,
                       seminaive_naive_table[m],
                       weights,
                       dctPlan);

      /* now load imaginary part of coefficients into output space */  
      memcpy(idataptr, fltres, sizeof(double) * (bw - m));
      
      idataptr += bw-m;
      
    }
    else{
      /* do real part */      
      Naive_AnalysisX(rres+(m*size),
                      bw,
                      m,
                      weights,
                      fltres,
                      seminaive_naive_table[m],
                      scratchpad );
      memcpy(rdataptr, fltres, sizeof(double) * (bw - m));
      rdataptr += bw-m;
      
      /* do imaginary part */
      Naive_AnalysisX(ires+(m*size),
                      bw,
                      m,
                      weights,
                      fltres,
                      seminaive_naive_table[m],
                      scratchpad );
      memcpy(idataptr, fltres, sizeof(double) * (bw - m));
      idataptr += bw-m;
    }
  }
  
  /*** now do upper coefficients ****/
  
  /* now if the data is real, we don't have to compute the
     coefficients whose order is less than 0, i.e. since
     the data is real, we know that
     f-hat(l,-m) = (-1)^m * conjugate(f-hat(l,m)),
     so use that to get the rest of the coefficients
     
     dataformat =0 -> samples are complex, =1 -> samples real
     
  */
  
  if( dataformat == 0 ){
    
    /* note that m is greater than bw here, but this is for
       purposes of indexing the input data arrays.  
       The "true" value of m as a parameter for Pml is
       size - m  */
    
    for (m=bw+1; m<size; m++) {
      /*
        fprintf(stderr,"m = %d\n",-(size-m));
      */
      
      if ( (size-m) < cutoff )
        {
          /* do real part */
          SemiNaiveReduced(rres+(m*size), 
                           bw, 
                           size-m, 
                           fltres,
                           scratchpad,
                           seminaive_naive_table[size-m],
                           weights,
                           dctPlan );
      
          /* now load real part of coefficients into output space */  
          if ((m % 2) != 0) {
            for (i=0; i<m-bw; i++)
              rdataptr[i] = -fltres[i];
          }
          else {
            memcpy(rdataptr, fltres, sizeof(double) * (m - bw));
          }
          rdataptr += m-bw;
      
          /* do imaginary part */
          SemiNaiveReduced(ires+(m*size), 
                           bw, 
                           size-m, 
                           fltres,
                           scratchpad,
                           seminaive_naive_table[size-m],
                           weights,
                           dctPlan);
      
          /* now load imag part of coefficients into output space */  
          if ((m % 2) != 0) {
            for (i=0; i<m-bw; i++)
              idataptr[i] = -fltres[i];
          }
          else {
            memcpy(idataptr, fltres, sizeof(double) * (m - bw));
          }
          idataptr += m-bw;
        }
      else
        {
          Naive_AnalysisX(rres+(m*size),
                          bw,
                          size-m,
                          weights,
                          fltres,
                          seminaive_naive_table[size-m],
                          scratchpad);

          /* now load real part of coefficients into output space */  
          if ((m % 2) != 0) {
            for (i=0; i<m-bw; i++)
              rdataptr[i] = -fltres[i];
          }
          else {
            memcpy(rdataptr, fltres, sizeof(double) * (m - bw));
          }
          rdataptr += m-bw;

          /* do imaginary part */
          Naive_AnalysisX(ires+(m*size),
                          bw,
                          size-m,
                          weights,
                          fltres,
                          seminaive_naive_table[size-m],
                          scratchpad);

          /* now load imag part of coefficients into output space */  
          if ((m % 2) != 0) {
            for (i=0; i<m-bw; i++)
              idataptr[i] = -fltres[i];
          }
          else {
            memcpy(idataptr, fltres, sizeof(double) * (m - bw));
          }
          idataptr += m-bw;

        }

    }
  }
  else        /**** if the data is real ****/
    {            
      pow_one = 1.0;
      for(i = 1; i < bw; i++){
        pow_one *= -1.0;
        for( j = i; j < bw; j++){
          /*
          SEANINDEXP(dummy0,i,j,bw);
          SEANINDEXN(dummy1,-i,j,bw);
          */

          dummy0 = seanindex(i, j, bw);
          dummy1 = seanindex(-i, j, bw);

          rcoeffs[dummy1] =
            pow_one * rcoeffs[dummy0];
          icoeffs[dummy1] =
            -pow_one * icoeffs[dummy0];
        }
      }
    }
  
}

Revision:	1294
Committed:	Thu Jun 24 13:44:53 2004 UTC (21 years, 4 months ago) by chrisfen
Content type:	text/plain
File size:	11574 byte(s)
Log Message:	Removed some unnecessary C functions from the SHT code
#	User	Rev	Content
1	chrisfen	1287	/***************************************************************************
2			**************************************************************************
3
4			S2kit 1.0
5
6			A lite version of Spherical Harmonic Transform Kit
7
8			Peter Kostelec, Dan Rockmore
9			{geelong,rockmore}@cs.dartmouth.edu
10
11			Contact: Peter Kostelec
12			geelong@cs.dartmouth.edu
13
14			Copyright 2004 Peter Kostelec, Dan Rockmore
15
16			This file is part of S2kit.
17
18			S2kit is free software; you can redistribute it and/or modify
19			it under the terms of the GNU General Public License as published by
20			the Free Software Foundation; either version 2 of the License, or
21			(at your option) any later version.
22
23			S2kit is distributed in the hope that it will be useful,
24			but WITHOUT ANY WARRANTY; without even the implied warranty of
25			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26			GNU General Public License for more details.
27
28			You should have received a copy of the GNU General Public License
29			along with S2kit; if not, write to the Free Software
30			Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
31
32			See the accompanying LICENSE file for details.
33
34			************************************************************************
35			************************************************************************/
36
37
38
39			/********************************************************************
40
41			FST_semi_memo.c - routines to perform convolutions on the
42			2-sphere using a combination of semi-naive and naive algorithms.
43
44			ASSUMES THAT ALL PRECOMPUTED-DATA IS IN MEMORY, AND NOT TO BE
45			READ FROM THE DISK.
46
47			The primary functions in this package are
48
49			1) FST_semi_memo() - computes the spherical harmonic expansion.
50			2) InvFST_semi_memo() - computes the inverse spherical harmonic transform.
51			3) FZT_semi_memo() - computes the zonal harmonic transform.
52			4) TransMult() - Multiplies harmonic coefficients using Driscoll-Healy
53			result. Dual of convolution in "time" domain.
54			5) Conv2Sphere_semi_memo() - Convolves two functins defined on the 2-sphere,
55			using seminaive transform
56
57			and one utility function:
58
59			1) seanindex(): Given bandwidth bw, seanindex(m,l,bw) will give the
60			position of the coefficient f-hat(m,l) in the one-row array
61
62			For descriptions on calling these functions, see the documentation
63			preceding each function.
64
65			*/
66
67			#include <math.h>
68			#include <stdio.h>
69			#include <stdlib.h>
70			#include <string.h>
71
72			#include "fftw3.h"
73
74			#include "makeweights.h"
75			#include "cospmls.h"
76			#include "primitive.h"
77			#include "naive_synthesis.h"
78			#include "seminaive.h"
79			#include "FST_semi_memo.h"
80
81
82
83			/************************************************************************/
84
85
86			/*****************************************************************
87
88			Given bandwidth bw, seanindex(m,l,bw) will give the position of the
89			coefficient f-hat(m,l) in the one-row array that Sean stores the spherical
90			coefficients. This is needed to help preserve the symmetry that the
91			coefficients have: (l = degree, m = order, and abs(m) <= l)
92
93			f-hat(l,-m) = (-1)^m * conjugate( f-hat(l,m) )
94
95			Thanks for your help Mark!
96
97			******************************************************************/
98
99
100			int seanindex(int m,
101			int l,
102			int bw)
103			{
104			int bigL;
105
106			bigL = bw - 1;
107
108			if( m >= 0 )
109			return( m * ( bigL + 1 ) - ( ( m * (m - 1) ) /2 ) + ( l - m ) );
110			else
111			return( ( ( bigL * ( bigL + 3 ) ) /2 ) + 1 +
112			( ( bigL + m ) * ( bigL + m + 1 ) / 2 ) + ( l - abs( m ) ) );
113			}
114
115
116			/************************************************************************/
117			/* performs a spherical harmonic transform using the semi-naive
118			and naive algorithms
119
120			bw -> bandwidth of problem
121			size -> size = 2*bw -> dimension of input array (recall that
122			sampling is done at twice the bandwidth)
123
124			The inputs rdata and idata are expected to be pointers to
125			size x size arrays. The array rdata contains the real parts
126			of the function samples, and idata contains the imaginary
127			parts.
128
129			rcoeffs and icoeffs are expected to be pointers to bw x bw arrays,
130			and will contain the harmonic coefficients in a "linearized" form.
131			The array rcoeffs contains the real parts of the coefficients,
132			and icoeffs contains the imaginary parts.
133
134			spharmonic_pml_table should be a (double **) pointer to
135			the result of a call to Spharmonic_Pml_Table. Because this
136			table is re-used in the inverse transform, and because for
137			timing purposes the computation of the table is not included,
138			it is passed in as an argument. Also, at some point this
139			code may be used as par of a series of convolutions, so
140			reducing repetitive computation is prioritized.
141
142			spharmonic_pml_table will be an array of (double *) pointers
143			the array being of length TableSize(m,bw)
144
145			workspace needs to be a double pointer to an array of size
146			(8 * bw^2) + (7 * bw).
147
148			cutoff -> what order to switch from semi-naive to naive
149			algorithm.
150
151			dataformat =0 -> samples are complex, =1 -> samples real
152
153
154			Output Ordering of coeffs f(m,l) is
155			f(0,0) f(0,1) f(0,2) ... f(0,bw-1)
156			f(1,1) f(1,2) ... f(1,bw-1)
157			etc.
158			f(bw-2,bw-2), f(bw-2,bw-1)
159			f(bw-1,bw-1)
160			f(-(bw-1),bw-1)
161			f(-(bw-2),bw-2) f(-(bw-2),bw-1)
162			etc.
163			f(-2,2) ... f(-2,bw-1)
164			f(-1,1) f(-1,2) ... f(-1,bw-1)
165
166
167			This only requires an array of size (bw*bw). If zero-padding
168			is used to make the indexing nice, then you need a an
169			(2bw-1) * bw array - but that is not done here.
170			Because of the amount of space necessary for doing
171			large transforms, it is important not to use any
172			more than necessary.
173
174			*/
175
176			void FST_semi_memo(double rdata, double idata,
177			double rcoeffs, double icoeffs,
178			int bw,
179			double **seminaive_naive_table,
180			double *workspace,
181			int dataformat,
182			int cutoff,
183			fftw_plan *dctPlan,
184			fftw_plan *fftPlan,
185			double *weights )
186			{
187			int size, m, i, j;
188			int dummy0, dummy1 ;
189			double rres, ires;
190			double rdataptr, idataptr;
191			double fltres, scratchpad;
192			double *eval_pts;
193			double pow_one;
194			double tmpA, tmpSize ;
195
196			size = 2*bw ;
197			tmpSize = 1./ ((double ) size);
198			tmpA = sqrt( 2. * M_PI ) ;
199
200			rres = workspace; /* needs (size * size) = (4 * bw^2) */
201			ires = rres + (size * size); /* needs (size * size) = (4 * bw^2) */
202			fltres = ires + (size * size) ; /* needs bw */
203			eval_pts = fltres + bw; /* needs (2bw) /
204			scratchpad = eval_pts + (2bw); / needs (4 * bw) */
205
206			/* total workspace is (8 * bw^2) + (7 * bw) */
207
208			/* do the FFTs along phi */
209			fftw_execute_split_dft( *fftPlan,
210			rdata, idata,
211			rres, ires ) ;
212			/*
213			normalize
214
215			note that I'm getting the sqrt(2pi) in there at
216			this point ... to account for the fact that the spherical
217			harmonics are of norm 1: I need to account for
218			the fact that the associated Legendres are
219			of norm 1
220			*/
221			tmpSize *= tmpA ;
222			for( j = 0 ; j < size*size ; j ++ )
223			{
224			rres[j] *= tmpSize ;
225			ires[j] *= tmpSize ;
226			}
227
228			/* point to start of output data buffers */
229			rdataptr = rcoeffs;
230			idataptr = icoeffs;
231
232			for (m=0; m<bw; m++) {
233			/*
234			fprintf(stderr,"m = %d\n",m);
235			*/
236
237			/* test to see if before cutoff or after */
238			if (m < cutoff){
239
240			/* do the real part */
241			SemiNaiveReduced(rres+(m*size),
242			bw,
243			m,
244			fltres,
245			scratchpad,
246			seminaive_naive_table[m],
247			weights,
248			dctPlan);
249
250			/* now load real part of coefficients into output space */
251			memcpy(rdataptr, fltres, sizeof(double) * (bw - m));
252
253			rdataptr += bw-m;
254
255			/* do imaginary part */
256			SemiNaiveReduced(ires+(m*size),
257			bw,
258			m,
259			fltres,
260			scratchpad,
261			seminaive_naive_table[m],
262			weights,
263			dctPlan);
264
265			/* now load imaginary part of coefficients into output space */
266			memcpy(idataptr, fltres, sizeof(double) * (bw - m));
267
268			idataptr += bw-m;
269
270			}
271			else{
272			/* do real part */
273			Naive_AnalysisX(rres+(m*size),
274			bw,
275			m,
276			weights,
277			fltres,
278			seminaive_naive_table[m],
279			scratchpad );
280			memcpy(rdataptr, fltres, sizeof(double) * (bw - m));
281			rdataptr += bw-m;
282
283			/* do imaginary part */
284			Naive_AnalysisX(ires+(m*size),
285			bw,
286			m,
287			weights,
288			fltres,
289			seminaive_naive_table[m],
290			scratchpad );
291			memcpy(idataptr, fltres, sizeof(double) * (bw - m));
292			idataptr += bw-m;
293			}
294			}
295
296			/* now do upper coefficients **/
297
298			/* now if the data is real, we don't have to compute the
299			coefficients whose order is less than 0, i.e. since
300			the data is real, we know that
301			f-hat(l,-m) = (-1)^m * conjugate(f-hat(l,m)),
302			so use that to get the rest of the coefficients
303
304			dataformat =0 -> samples are complex, =1 -> samples real
305
306			*/
307
308			if( dataformat == 0 ){
309
310			/* note that m is greater than bw here, but this is for
311			purposes of indexing the input data arrays.
312			The "true" value of m as a parameter for Pml is
313			size - m */
314
315			for (m=bw+1; m<size; m++) {
316			/*
317			fprintf(stderr,"m = %d\n",-(size-m));
318			*/
319
320			if ( (size-m) < cutoff )
321			{
322			/* do real part */
323			SemiNaiveReduced(rres+(m*size),
324			bw,
325			size-m,
326			fltres,
327			scratchpad,
328			seminaive_naive_table[size-m],
329			weights,
330			dctPlan );
331
332			/* now load real part of coefficients into output space */
333			if ((m % 2) != 0) {
334			for (i=0; i<m-bw; i++)
335			rdataptr[i] = -fltres[i];
336			}
337			else {
338			memcpy(rdataptr, fltres, sizeof(double) * (m - bw));
339			}
340			rdataptr += m-bw;
341
342			/* do imaginary part */
343			SemiNaiveReduced(ires+(m*size),
344			bw,
345			size-m,
346			fltres,
347			scratchpad,
348			seminaive_naive_table[size-m],
349			weights,
350			dctPlan);
351
352			/* now load imag part of coefficients into output space */
353			if ((m % 2) != 0) {
354			for (i=0; i<m-bw; i++)
355			idataptr[i] = -fltres[i];
356			}
357			else {
358			memcpy(idataptr, fltres, sizeof(double) * (m - bw));
359			}
360			idataptr += m-bw;
361			}
362			else
363			{
364			Naive_AnalysisX(rres+(m*size),
365			bw,
366			size-m,
367			weights,
368			fltres,
369			seminaive_naive_table[size-m],
370			scratchpad);
371
372			/* now load real part of coefficients into output space */
373			if ((m % 2) != 0) {
374			for (i=0; i<m-bw; i++)
375			rdataptr[i] = -fltres[i];
376			}
377			else {
378			memcpy(rdataptr, fltres, sizeof(double) * (m - bw));
379			}
380			rdataptr += m-bw;
381
382			/* do imaginary part */
383			Naive_AnalysisX(ires+(m*size),
384			bw,
385			size-m,
386			weights,
387			fltres,
388			seminaive_naive_table[size-m],
389			scratchpad);
390
391			/* now load imag part of coefficients into output space */
392			if ((m % 2) != 0) {
393			for (i=0; i<m-bw; i++)
394			idataptr[i] = -fltres[i];
395			}
396			else {
397			memcpy(idataptr, fltres, sizeof(double) * (m - bw));
398			}
399			idataptr += m-bw;
400
401			}
402
403			}
404			}
405			else /** if the data is real **/
406			{
407			pow_one = 1.0;
408			for(i = 1; i < bw; i++){
409			pow_one *= -1.0;
410			for( j = i; j < bw; j++){
411			/*
412			SEANINDEXP(dummy0,i,j,bw);
413			SEANINDEXN(dummy1,-i,j,bw);
414			*/
415
416			dummy0 = seanindex(i, j, bw);
417			dummy1 = seanindex(-i, j, bw);
418
419			rcoeffs[dummy1] =
420			pow_one * rcoeffs[dummy0];
421			icoeffs[dummy1] =
422			-pow_one * icoeffs[dummy0];
423			}
424			}
425			}
426
427			}
428