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 (20 years ago) by chrisfen
Content type:	text/plain
File size:	11574 byte(s)
Log Message:	Removed some unnecessary C functions from the SHT code
#	Content
1	/***************************************************************************
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