AcopFFT xref

View Javadoc

1   /*
2    * Copyright (c) 2006 Stiftung Deutsches Elektronen-Synchroton,
3    * Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY.
4    *
5    * THIS SOFTWARE IS PROVIDED UNDER THIS LICENSE ON AN "../AS IS" BASIS.
6    * WITHOUT WARRANTY OF ANY KIND, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
7    * TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR PARTICULAR PURPOSE AND
8    * NON-INFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
9    * FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10   * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
11   * THE USE OR OTHER DEALINGS IN THE SOFTWARE. SHOULD THE SOFTWARE PROVE DEFECTIVE
12   * IN ANY RESPECT, THE USER ASSUMES THE COST OF ANY NECESSARY SERVICING, REPAIR OR
13   * CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL PART OF THIS LICENSE.
14   * NO USE OF ANY SOFTWARE IS AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
15   * DESY HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS,
16   * OR MODIFICATIONS.
17   * THE FULL LICENSE SPECIFYING FOR THE SOFTWARE THE REDISTRIBUTION, MODIFICATION,
18   * USAGE AND OTHER RIGHTS AND OBLIGATIONS IS INCLUDED WITH THE DISTRIBUTION OF THIS
19   * PROJECT IN THE FILE LICENSE.HTML. IF THE LICENSE IS NOT INCLUDED YOU MAY FIND A COPY
20   * AT HTTP://WWW.DESY.DE/LEGAL/LICENSE.HTM
21   */
22  
23  package de.desy.acop.chart;
24  
25  public class AcopFFT
26  {
27    private int data_size, half_size, ary_size, reverseFFT;
28  	private float q_wert;
29  	private float[]	q_array, f_array, adr0_r, adr0_i, adr1_r, adr1_i;
30  
31    public AcopFFT(){};
32  
33    public int fFT(int n_points, double[] x_array, double[] outf_array, double[] maxPerCut,  int specktrum)
34    { int i,j, oddFlag;
35  	  int[] error = new int[1];
36         
37      if ( n_points < 2 || n_points > AcopConst.max_points ) return -1;
38        
39      if ( q_array == null || n_points > ary_size )
40      {
41          q_array = new float[n_points];
42          f_array = new float[n_points];
43          adr0_r = new float[n_points*2];
44          adr0_i = new float[n_points*2];
45          adr1_r = new float[n_points*2];
46          adr1_i = new float[n_points*2];
47          ary_size = n_points;
48      }
49  	  else for (i=0; i<n_points; i++)
50      {
51          adr0_r[i] = 0;
52          adr0_i[i] = 0;
53          adr1_r[i] = 0;
54          adr1_i[i] = 0;
55      }
56      for (j=0; j<n_points; j++) q_array[j] = (float)x_array[j];   /* input data points    */
57      half_size = n_points / 2;
58      oddFlag = n_points % 2;
59        
60  	  // Fast FT if possible
61      data_size = n_points ;
62      reverseFFT = specktrum == AcopConst.FFT_RI_INVERSE ? 1 : 0;
63      c_fourier ( data_size, 0) ;     // FFT first call = 0 (have to) */
64  
65  
66      // 0x1: spektrum + 0; 0x2: phase + 0; 0x3: spek + phase; else: real + imaginary part
67      if ( specktrum < AcopConst.FFT_SP || specktrum==AcopConst.FFT_REAL || specktrum==AcopConst.FFT_IMAGINARY )
68        for (i=half_size; i<n_points; i++)
69        {
70          f_array[i] = 0;
71          outf_array[i] = 0;
72        }
73        if ( oddFlag != 0 ) // last point
74        {
75          f_array[n_points-1]=0;
76          f_array[n_points-1]=0;
77        }
78  
79      for (j=0; j<=half_size; j++)	// first half => FFT Power spectrum
80  	  {								// second half => FFT phase (-pi -> pi)
81  	    if ( (specktrum == AcopConst.FFT_SPEKTRUM) || (specktrum == AcopConst.FFT_SP) )
82  	    {
83  	      f_array[j] = (float)Math.sqrt(adr0_r[j]*adr0_r[j]+adr0_i[j]*adr0_i[j]);
84  		    if ( specktrum == AcopConst.FFT_SP && j>0 && j+half_size<n_points ) f_array[j+half_size] = (float)Math.atan2(adr0_i[j],adr0_r[j]);
85  	    }
86  	    else if ( specktrum == AcopConst.FFT_PHASE ) f_array[j] = (float)Math.atan2(adr0_i[j],adr0_r[j]);
87        else if ( specktrum == AcopConst.FFT_RI )
88        {
89          f_array[j] = adr0_r[j];
90          if ( j>0 && j+half_size<n_points ) f_array[j+half_size] = adr0_i[j];
91        }
92        else if ( specktrum == AcopConst.FFT_REAL ) f_array[j] = adr0_r[j];
93        else if ( specktrum == AcopConst.FFT_IMAGINARY ) f_array[j] = adr0_i[j];
94        else
95        {
96          f_array[j] = adr0_r[j];
97          if ( j>0 && j+half_size<n_points ) f_array[j+half_size] = adr0_r[j+half_size];
98        }
99  	  }
100     for (j=0; j<n_points; j++)
101     {
102       if ( (specktrum < AcopConst.FFT_SP || specktrum==AcopConst.FFT_REAL || specktrum==AcopConst.FFT_IMAGINARY) && j > half_size )break;
103       outf_array[j] = (double)f_array[j];
104     } 
105        
106 	  maxPerCut[0] = max_find ( half_size, error, maxPerCut[1], maxPerCut[2]);   // find Q value
107 	  
108 	  return 0;
109   }
110     
111   private double max_find(int size, int[] error, double percent, double cut)
112   { int j, i_max ;
113 	  double r_max, r_cut, i_weight ;
114 
115 	  error[0] = 0 ;
116     r_max =  f_array[0] ;       /* find maximum position i_max  */
117     i_max = 0;
118 	  for (j=0; j<=size; j++)  if ( f_array[j] > r_max)
119 	  {
120 	  	r_max =  f_array[j] ;
121 		  i_max = j ;
122 	  }
123 	  r_cut = r_max * percent ;       /* above this wert, take into account   */
124 	  r_max = 0.0;                            // r_max = -1000.0 ;
125 	  i_weight = 0.0 ;
126 	  for (j=0; j<=size; j++)         /* make weight to all high points       */
127 	    if ( f_array[j] > r_cut) 
128 	  {
129   		r_max +=  f_array[j] ;
130 	  	i_weight +=  f_array[j] * j ;
131 	  }
132 	  if ( r_max < AcopConst.zero_check)
133 	  {
134 	    i_weight = i_max;
135 	    error[0] |= 0x1;
136     }
137     else i_weight /= r_max ;
138 	
139 	  if ( Math.abs(i_weight - i_max) > cut ) error[0] |= 0x2 ;
140 	  
141 	  return (i_weight) ;                                             /* return Q*freqency    */
142   }
143 
144   private void c_fourier(int now_data, int n_call)
145   { int i, j, flag;
146 	  int off_set ;                      /* memory off set       */
147 	  int h_now = now_data/2 ;
148 	  int r_base_adr_index, r_base_adr_flag;
149 	  float f_sin, f_cos, temp_r, temp_i, tmp ;
150 
151     // if even number, do seperation 
152     if ( now_data % 2 == 0 && now_data > 2 )       /* seperate even and odd points */
153 	  {
154 	  	c_fourier (h_now, n_call) ;                     /* call even    */
155 		  c_fourier (h_now, n_call + data_size/now_data) ;    /* odd  */
156       off_set = now_data;     /* destination offset position */
157       r_base_adr_flag = 0;
158       r_base_adr_index = 0;
159 	  	if ( data_size > now_data)         /* if not first call    */
160 		  {
161 		    if (n_call >=  half_size / now_data ) r_base_adr_flag = 1; /* choose memory odd*/
162         r_base_adr_index += off_set;
163 		  }
164 		  off_set = now_data / 2 ;           /* lower step offset, get source*/
165 
166 		  f_cos = (float)1.0 ;
167 		  f_sin = (float)0.0 ;
168 	  	for (i=0; i<h_now; i++)               /* calculate f1k, f2k   */
169 	  	{
170 		    if (i != 0 )
171 		    {
172 			    f_sin = (float) Math.PI * i / h_now ;                        /* shift angle  */
173           if ( reverseFFT != 0 ) f_sin = -f_sin;
174 			    f_cos = (float) Math.cos ( (double) f_sin) ;
175 			    f_sin = (float) Math.sin ( (double) f_sin) ;
176 		    }
177 		    temp_r =  adr1_r [i+ off_set] * f_cos + adr1_i [i+ off_set] * f_sin ;
178 		    temp_i = -adr1_r [i+ off_set] * f_sin + adr1_i [i+ off_set] * f_cos ;
179 
180         if ( r_base_adr_flag == 0 )
181         {
182           adr0_r [i+r_base_adr_index] = (adr0_r [i+off_set] + temp_r) / 2 ;  /* first half f(k)       */
183           adr0_i [i+r_base_adr_index] = (adr0_i [i+ off_set] + temp_i) /2 ;
184           if ( reverseFFT != 0 )
185           {
186             adr0_r [i+r_base_adr_index] *= 2;
187             adr0_i [i+r_base_adr_index] *= 2;
188           }
189         }
190         else
191         {
192           adr1_r [i+r_base_adr_index] = (adr0_r [i+off_set] + temp_r) /2 ;  /* first half f(k)       */
193           adr1_i [i+r_base_adr_index] = (adr0_i [i+ off_set] + temp_i) /2 ;
194           if ( reverseFFT != 0 )
195           {
196             adr1_r [i+r_base_adr_index] *= 2;
197             adr1_i [i+r_base_adr_index] *= 2;
198           }
199         }
200         if ( data_size != now_data || reverseFFT != 0 )  // not first call, do f(n+k)
201         {
202           if ( r_base_adr_flag == 0 )
203           {
204             adr0_r [i+r_base_adr_index + h_now ] =  (adr0_r [i+off_set] - temp_r) / 2 ; // f(n+k)
205             if ( reverseFFT != 0 ) adr0_r [i+r_base_adr_index + h_now ] *= 2;
206           }
207           else
208           {
209             adr1_r [i+r_base_adr_index + h_now ] =  (adr0_r [i+off_set] - temp_r) / 2 ; // f(n+k)
210             if ( reverseFFT != 0 ) adr1_r [i+r_base_adr_index + h_now ] *= 2;
211           }
212 
213           if ( data_size != now_data) // no imaginay for first call
214           {
215             if ( r_base_adr_flag == 0 )
216             {
217               adr0_i [i+r_base_adr_index + h_now ] =  (adr0_i [i+off_set] - temp_i) / 2 ;
218               if ( reverseFFT != 0 ) adr0_i [i+r_base_adr_index + h_now ] *= 2;
219             }
220             else
221             {
222               adr1_i [i+r_base_adr_index + h_now ] =  (adr0_i [i+off_set] - temp_i) / 2 ;
223               if ( reverseFFT != 0 ) adr1_i [i+r_base_adr_index + h_now ] *= 2;
224             }
225           }
226         }
227 		  }
228 
229       if ( data_size == now_data)    /*  if first call, do middle point   */
230       {
231         adr0_r[half_size] = 0;   // init to 0
232         adr0_i[half_size] = 0;
233 
234         for (i=0; i<data_size; i++)
235         {
236           if ( reverseFFT != 0 && i>h_now) j=data_size -i;
237           else j=i;
238           if ( i%2 == 0 ) adr0_r[half_size] += q_array[j];
239           else adr0_r[half_size] -= q_array[j];
240         }
241         if ( reverseFFT == 0 ) adr0_r[half_size] /= data_size;
242       }
243 
244 	  }
245 	  else
246     {                        /* only last two points, or odd points */
247       if ( n_call * now_data * 2 < data_size ) r_base_adr_flag = 0;
248       else r_base_adr_flag = 1;
249       r_base_adr_index = 0;
250        if ( data_size != now_data ) r_base_adr_index = now_data; // add offset
251 
252        if ( data_size != now_data && now_data == 2 ) // after seperation to 2 points
253        {
254          if ( reverseFFT == 0 )
255          {
256            if ( r_base_adr_flag == 0 )
257            {
258              adr0_r[r_base_adr_index] = (float) ( q_array[n_call]+ q_array[n_call + half_size]) / 2;
259              adr0_r[r_base_adr_index + 1] = (float) ( q_array[n_call]- q_array[n_call + half_size]) / 2;
260            }
261            else
262            {
263              adr1_r[r_base_adr_index] = (float) ( q_array[n_call]+ q_array[n_call + half_size]) / 2;
264              adr1_r[r_base_adr_index + 1] = (float) ( q_array[n_call]- q_array[n_call + half_size]) / 2;
265            }
266          }
267          else // reverseFFT, get imaginary part
268          {
269            i = n_call+ half_size;
270            if ( i>half_size )  // second half
271            {
272              i = data_size - i ;
273              flag=1;
274            }
275            else flag=0;
276 
277            if ( r_base_adr_flag == 0 )
278            {
279              adr0_r[r_base_adr_index] = (float) ( q_array[n_call]+ q_array[i]);
280              adr0_r[r_base_adr_index + 1] = (float) ( q_array[n_call]- q_array[i]);
281            }
282            else
283            {
284              adr1_r[r_base_adr_index] = (float) ( q_array[n_call]+ q_array[i]);
285              adr1_r[r_base_adr_index + 1] = (float) ( q_array[n_call]- q_array[i]);
286            }
287            
288            if ( i + half_size < data_size )
289            {
290              if(flag==0)
291              {
292                if ( r_base_adr_flag == 0 )
293                {
294                  adr0_i[r_base_adr_index] = (float) q_array[i+ half_size] ;
295                  adr0_i[r_base_adr_index + 1] = (float) -q_array[i+ half_size] ;
296                }
297                else
298                {
299                  adr1_i[r_base_adr_index] = (float) q_array[i+ half_size] ;
300                  adr1_i[r_base_adr_index + 1] = (float) -q_array[i+ half_size] ;
301                }
302              }
303              else
304              {
305                if ( r_base_adr_flag == 0 )
306                {
307                  adr0_i[r_base_adr_index] = (float) -q_array[i+ half_size] ;
308                  adr0_i[r_base_adr_index + 1] = (float) q_array[i+ half_size] ;
309                }
310                else
311                {
312                  adr1_i[r_base_adr_index] = (float) -q_array[i+ half_size] ;
313                  adr1_i[r_base_adr_index + 1] = (float) q_array[i+ half_size] ;
314                }
315              }
316            }
317 
318            if ( n_call != 0 ) // first entry no imaginary
319            {
320              if ( r_base_adr_flag == 0 )
321              {
322                adr0_i[r_base_adr_index] += q_array[n_call+half_size] ;
323                adr0_i[r_base_adr_index+1] += q_array[n_call+half_size] ;
324              }
325              else
326              {
327                adr1_i[r_base_adr_index] += q_array[n_call+half_size] ;
328                adr1_i[r_base_adr_index+1] += q_array[n_call+half_size] ;
329              }
330            }
331          }
332        }
333        else
334        {
335          if ( data_size == now_data )
336          {
337            if ( reverseFFT != 0 ) h_now = data_size;
338            else h_now += 1 ;  // complete normal FT (odd number)
339          }
340          else h_now = now_data;       // odd data
341 
342          f_cos = (float)1.0 ; //init
343          f_sin = (float)0.0 ;
344          for (i=0; i<h_now; i++)  // odd number data, half of whole data
345          {
346            if ( r_base_adr_flag == 0 )
347            {
348              adr0_r[r_base_adr_index+i] = 0;   // init
349              adr0_i[r_base_adr_index+i] = 0;
350            }
351            else
352            {
353              adr1_r[r_base_adr_index+i] = 0;   // init
354              adr1_i[r_base_adr_index+i] = 0;
355            }
356            for (j=0; j<now_data; j++) //FT
357            {
358              if (j==0)  //init
359              {
360                f_cos = (float)1.0 ;
361                f_sin = (float)0.0 ;
362              }
363              else if ( i != 0 && j != 0 ) //get sin and cos
364              {
365                f_sin = (float) (-2 * Math.PI * i * j / now_data) ; /* shift angle  */
366                if ( reverseFFT != 0 ) f_sin = -f_sin;
367                f_cos = (float) Math.cos ( (double) f_sin) ;
368                f_sin = (float) Math.sin ( (double) f_sin) ;
369              }
370              off_set = n_call + data_size / now_data * j;  // offset
371              flag=0;
372              if ( reverseFFT != 0 && off_set > half_size )
373              {
374                off_set = data_size - off_set ;
375                flag=1;
376              }
377              if (i==0 || j==0)
378              {
379                if ( r_base_adr_flag == 0 )
380                  adr0_r[r_base_adr_index+i] += (float)(q_array[off_set]);   // init
381                else
382                  adr1_r[r_base_adr_index+i] += (float)(q_array[off_set]);   // init               
383              }
384              else
385              {
386                if ( r_base_adr_flag == 0 )
387                  adr0_r[r_base_adr_index+i] += (float)(q_array[off_set] * f_cos );
388                else
389                  adr1_r[r_base_adr_index+i] += (float)(q_array[off_set] * f_cos );
390              }
391              if ( ( reverseFFT == 0 && now_data != 2 ) || (reverseFFT != 0 && data_size != now_data) )
392              {
393                if ( r_base_adr_flag == 0 )
394                  adr0_i[r_base_adr_index+i] += (float)(q_array[off_set] * f_sin );
395                else
396                  adr1_i[r_base_adr_index+i] += (float)(q_array[off_set] * f_sin );
397              }
398              if ( reverseFFT == 0 ) continue;
399              if ( off_set == 0 || (off_set == half_size && data_size %2 == 0) )continue;
400              //imaginary has different sign
401              if (flag==0)
402              {
403                if ( r_base_adr_flag == 0 )
404                  adr0_r[r_base_adr_index+i] -= (float)(q_array[off_set + half_size] * f_sin );
405                else
406                  adr1_r[r_base_adr_index+i] -= (float)(q_array[off_set + half_size] * f_sin );
407              }
408              else
409              {
410                if ( r_base_adr_flag == 0 )
411                  adr0_r[r_base_adr_index+i] += (float)(q_array[off_set + half_size] * f_sin );
412                else
413                  adr1_r[r_base_adr_index+i] += (float)(q_array[off_set + half_size] * f_sin );
414              }
415              if (data_size != now_data)
416              {
417                if (flag==0)
418                {
419                  if ( r_base_adr_flag == 0 )
420                    adr0_i[r_base_adr_index+i] += (float)(q_array[off_set + half_size] * f_cos );
421                  else
422                    adr1_i[r_base_adr_index+i] += (float)(q_array[off_set + half_size] * f_cos );
423                }
424                else
425                {
426                  if ( r_base_adr_flag == 0 )
427                    adr0_i[r_base_adr_index+i] -= (float)(q_array[off_set + half_size] * f_cos );
428                  else
429                    adr1_i[r_base_adr_index+i] -= (float)(q_array[off_set + half_size] * f_cos );
430                }
431              }
432            }
433            if (reverseFFT == 0)
434            {
435              if ( r_base_adr_flag == 0 )
436                adr0_r[r_base_adr_index+i] /= now_data;
437              else
438                adr1_r[r_base_adr_index+i] /= now_data;
439              if (now_data!=2)
440              {
441                if ( r_base_adr_flag == 0 )
442                  adr0_i[r_base_adr_index+i] /= now_data;
443                else
444                  adr1_i[r_base_adr_index+i] /= now_data;
445              }
446            }
447          }
448        }
449      }
450 
451   }
452 
453 }