1 /*
2 fft.c -- Fast Fourier Transform.
3
4 Copyright (C) 1994-2003 K. Scott Hunziker.
5 Copyright (C) 1990-1994 The Boeing Company.
6
7 See the file COPYING for license, warranty, and permission details.
8 */
9
10 static char rcsid[] =
11 "$Id: fft.c,v 1.7 2003/09/21 04:47:16 ksh Exp $";
12
13 #include "fft.h"
14 #include "entity.h"
15 #include "scalar.h"
16 #include "vector.h"
17 #include "matrix.h"
18 #include "table.h"
19 #include "cast.h"
20 #include "dense.h"
21 #include "transpose.h"
22
23 static ENTITY *PROTO (modify_fft_labels, (VECTOR * p));
24
25 static char *fft_options[] =
26 {
27 "row",
28 "estimate",
29 "measure",
30 "patient",
31 "exhaustive",
32 "forget",
33 NULL,
34 };
35 enum
36 {
37 FFT_OPT_ROW,
38 FFT_OPT_ESTIMATE,
39 FFT_OPT_MEASURE,
40 FFT_OPT_PATIENT,
41 FFT_OPT_EXHAUSTIVE,
42 FFT_OPT_FORGET,
43 };
44
45 #if HAVE_FFTW
46
47 static void PROTO (hermcpy, (COMPLEX *v, int x));
48
49 ENTITY *
bi_fftw(n,p,d)50 bi_fftw (n, p, d)
51 int n;
52 ENTITY *p;
53 ENTITY *d;
54 {
55 int dim, rigor;
56 ENTITY * volatile s = NULL;
57
58 EASSERT (p, 0, 0);
59
60 if (d != NULL && d->class != table)
61 {
62 fail ("Invalid options -- not a table.");
63 delete_2_entities (p, d);
64 raise_exception ();
65 }
66
67 if (!OK_OPTIONS ((TABLE *) d, fft_options)) raise_exception ();
68
69 dim = !IN_TABLE (d, fft_options[FFT_OPT_ROW]);
70
71 if (IN_TABLE (d, fft_options[FFT_OPT_FORGET]))
72 fftw_forget_wisdom();
73
74 rigor = FFTW_ESTIMATE;
75 if (IN_TABLE (d, fft_options[FFT_OPT_MEASURE]))
76 rigor = FFTW_MEASURE;
77 if (IN_TABLE (d, fft_options[FFT_OPT_PATIENT]))
78 rigor = FFTW_PATIENT;
79 if (IN_TABLE (d, fft_options[FFT_OPT_EXHAUSTIVE]))
80 rigor = FFTW_EXHAUSTIVE;
81
82 WITH_HANDLING
83 {
84 switch (p->class)
85 {
86 case scalar:
87 s = fftw_vector ((VECTOR *) scalar_to_vector ((SCALAR *) EAT(p)),
88 rigor);
89 break;
90
91 case vector:
92 s = fftw_vector ((VECTOR *) EAT(p), rigor);
93 break;
94
95 case matrix:
96
97 /* FFT of rows (dim=0); FFT of cols (dim=1) */
98 s = fftw_matrix ((MATRIX *) EAT(p), dim, rigor);
99 break;
100
101 default:
102 BAD_CLASS (p->class);
103 raise_exception ();
104 }
105 }
106 ON_EXCEPTION
107 {
108 delete_2_entities (p, d);
109 }
110 END_EXCEPTION;
111
112 delete_entity (d);
113 return s;
114 }
115
116 ENTITY *
fftw_vector(p,rigor)117 fftw_vector (p, rigor)
118 VECTOR *p;
119 int rigor;
120 {
121 REAL *tmpd = NULL;
122 COMPLEX *tmpc = NULL;
123 VECTOR *c = NULL;
124 fftw_plan plan = NULL;
125 int i, j, nn;
126
127 EASSERT (p, vector, 0);
128
129 WITH_HANDLING
130 {
131
132 p = (VECTOR *) dense_vector (EAT (p));
133
134 switch (p->type)
135 {
136 case integer:
137 p = (VECTOR *) cast_vector_integer_real (EAT (p));
138 /* no break fall through */
139
140 case real:
141 c = (VECTOR *) cast_vector_real_complex ((VECTOR *) copy_vector (p));
142
143 nn = p->nn;
144
145 if (nn == 0)
146 goto out; /* done if zero length */
147
148 if (rigor == FFTW_ESTIMATE)
149 {
150 /* Create the fftw plan */
151 plan = fftw_plan_dft_r2c_1d (nn, p->a.real,
152 (fftw_complex *) c->a.real, rigor);
153
154 if (plan == NULL) {
155 fail ("FFTW failed to make a plan.");
156 raise_exception ();
157 }
158
159 /* Execute the transform */
160 fftw_execute (plan);
161 }
162 else
163 {
164 tmpd = MALLOC (nn * sizeof (REAL));
165
166 /* Create the fftw plan */
167 plan = fftw_plan_dft_r2c_1d (nn, tmpd,
168 (fftw_complex *) c->a.real, rigor);
169
170 if (plan == NULL) {
171 fail ("FFTW failed to make a plan.");
172 raise_exception ();
173 }
174
175 memcpy (tmpd, p->a.real, nn * sizeof (REAL));
176
177 /* Execute the transform */
178 fftw_execute (plan);
179
180 FREE (EAT (tmpd));
181 }
182
183 /* Hermitian copy and reorder frequencies */
184 hermcpy (c->a.complex, nn);
185
186 break;
187
188 case complex:
189
190 c = (VECTOR *) dup_vector (EAT (p));
191
192 nn = c->nn;
193
194 if (nn == 0)
195 goto out; /* done if zero length */
196
197 tmpc = MALLOC (2 * nn * sizeof (REAL));
198
199 /* Create the fftw plan */
200 plan = fftw_plan_dft_1d(nn, (fftw_complex *) tmpc,
201 (fftw_complex *) tmpc, FFTW_FORWARD, rigor);
202
203 if (plan == NULL) {
204 fail ("FFTW failed to make a plan.");
205 raise_exception ();
206 }
207
208 /* Copy the data. */
209 memcpy (tmpc, c->a.complex, 2*nn*sizeof(REAL));
210
211 /* Execute the transform */
212 fftw_execute (plan);
213
214 /* copy with reorder */
215 for (i=0,j=(nn-1)/2; j<nn; i++,j++) {
216 c->a.complex[j].real = tmpc[i].real;
217 c->a.complex[j].imag = tmpc[i].imag;
218 }
219 for (i=(nn+2)/2,j=0; i<nn; i++,j++) {
220 c->a.complex[j].real = tmpc[i].real;
221 c->a.complex[j].imag = tmpc[i].imag;
222 }
223
224 FREE (EAT (tmpc));
225
226 break;
227
228 default:
229 fail ("Can't apply \"fft\" to a %s vector.", type_string[p->type]);
230 raise_exception ();
231 }
232
233 if (c->eid)
234 c->eid = modify_fft_labels ((VECTOR *) EAT (c->eid));
235 }
236
237 out: /* can't jump past the ON_EXCEPTION */
238
239 ON_EXCEPTION
240 {
241 if (plan != NULL) fftw_destroy_plan (plan);
242 TFREE (tmpc);
243 TFREE (tmpd);
244 delete_2_vectors (p, c);
245 }
246 END_EXCEPTION;
247
248 fftw_destroy_plan (plan);
249 delete_vector (p);
250 return ENT (c);
251 }
252
253 ENTITY *
fftw_matrix(p,dim,rigor)254 fftw_matrix (p, dim, rigor)
255 MATRIX *p;
256 int dim; /* FFT of rows (dim=0); FFT of cols (dim=1) */
257 int rigor;
258 {
259 REAL *tmpd = NULL;
260 COMPLEX *tmpc = NULL;
261 fftw_plan plan = NULL;
262 int i, j, l, nr, nc, nn;
263
264 EASSERT (p, matrix, 0);
265
266 WITH_HANDLING
267 {
268 p = (MATRIX *) dense_matrix (EAT (p));
269
270 switch (p->type)
271 {
272 case integer:
273 p = (MATRIX *) cast_matrix_integer_real (EAT (p));
274 /* no break fall through */
275
276 case real:
277 p = (MATRIX *) cast_matrix_real_complex (EAT (p));
278
279 if ( dim == 0 ) /* Take the FFT of each row */
280 {
281 nr = p->nr; /* Number of rows */
282 nn = p->nc; /* Number to columns (elements per row) */
283
284 if (nn == 0)
285 goto out; /* Done if columns have zero length */
286
287 if (nr)
288 {
289 /* Allocate space - transform not done in place */
290 tmpc = MALLOC (2 * nn * sizeof(REAL));
291 tmpd = MALLOC (nn * sizeof(REAL));
292
293 /* Create the fftw plan */
294 plan = fftw_plan_dft_r2c_1d(nn, tmpd,
295 (fftw_complex *) tmpc, rigor);
296
297 if (plan == NULL) {
298 fail ("FFTW failed to make a plan.");
299 raise_exception ();
300 }
301
302 for (l = 0; l < nr; l++)
303 {
304 /* Copy the data */
305 for (i=0; i<nn; i++) {
306 tmpd[i] = p->a.complex[l+i*nr].real;
307 }
308
309 /* Execute the transform */
310 fftw_execute(plan);
311
312 /* Hermitian copy and reorder frequencies */
313 hermcpy((COMPLEX *) tmpc, nn);
314
315 /* copy the data */
316 for (i=0; i<nn; i++) {
317 p->a.complex[l+i*nr].real = tmpc[i].real;
318 p->a.complex[l+i*nr].imag = tmpc[i].imag;
319 }
320
321 } /* l - done loop over each row */
322
323 /* Destroy plan and tmpc arrays */
324 FREE (EAT (tmpd));
325 FREE (EAT (tmpc));
326 }
327
328 if (p->cid)
329 p->cid = modify_fft_labels ((VECTOR *) EAT (p->cid));
330 }
331 else /* Take the FFT of each column */
332 {
333 nc = p->nc; /* Number of columns */
334 nn = p->nr; /* Number to rows (elements per column) */
335
336 if (nn == 0)
337 goto out; /* Done if rows have zero length */
338
339 if (nc)
340 {
341 /* Allocate space - transform not done in place */
342 tmpc = MALLOC (2 * nn * sizeof(REAL));
343 tmpd = MALLOC (nn * sizeof(REAL));
344
345 /* Create the fftw plan */
346 plan = fftw_plan_dft_r2c_1d(nn, tmpd,
347 (fftw_complex *) tmpc, rigor);
348
349 if (plan == NULL) {
350 fail ("FFTW failed to make a plan.");
351 raise_exception ();
352 }
353
354 for (l = 0; l < nc; l++)
355 {
356 /* Copy the data */
357 for (i=0; i<nn; i++) {
358 tmpd[i] = p->a.complex[i+l*nn].real;
359 }
360
361 /* Execute the transform */
362 fftw_execute(plan);
363
364 /* Hermitian copy and reorder frequencies */
365 hermcpy((COMPLEX *) tmpc, nn);
366
367 /* copy the data */
368 for (i=0; i<nn; i++) {
369 p->a.complex[i+l*nn].real = tmpc[i].real;
370 p->a.complex[i+l*nn].imag = tmpc[i].imag;
371 }
372
373 } /* l - done loop over each column */
374
375 FREE (EAT (tmpd));
376 FREE (EAT (tmpc));
377 }
378
379 if (p->rid)
380 p->rid = modify_fft_labels ((VECTOR *) EAT (p->rid));
381
382 }
383
384 break;
385
386 case complex:
387
388 p = (MATRIX *) dup_matrix (EAT (p));
389
390 if ( dim == 0 ) /* Take the FFT of each row */
391 {
392 nr = p->nr; /* Number of rows */
393 nn = p->nc; /* Number to columns (elements per row) */
394
395 if (nn == 0)
396 goto out; /* Done if columns have zero length */
397
398 if (nr)
399 {
400 /* Allocate space - transform done in place */
401 tmpc = MALLOC (2 * nn * sizeof(REAL));
402
403 /* Create the fftw plan */
404 plan = fftw_plan_dft_1d (nn, (fftw_complex *) tmpc,
405 (fftw_complex *) tmpc, FFTW_FORWARD,
406 rigor);
407
408 if (plan == NULL) {
409 fail ("FFTW failed to make a plan.");
410 raise_exception ();
411 }
412
413 for (l = 0; l < nr; l++)
414 {
415 /* Copy the data */
416 for (i=0; i<nn; i++) {
417 tmpc[i].real = p->a.complex[l+i*nr].real;
418 tmpc[i].imag = p->a.complex[l+i*nr].imag;
419 }
420
421 /* Execute the transform */
422 fftw_execute(plan);
423
424 /* copy with reorder */
425 for (i=0,j=l+((nn-1)/2)*nr; i<=nn/2; i++,j+=nr) {
426 p->a.complex[j].real = tmpc[i].real;
427 p->a.complex[j].imag = tmpc[i].imag;
428 }
429 for (i=(nn+2)/2,j=l; i<nn; i++,j+=nr) {
430 p->a.complex[j].real = tmpc[i].real;
431 p->a.complex[j].imag = tmpc[i].imag;
432 }
433
434 } /* l - done loop over each row */
435
436 FREE (EAT (tmpc));
437 }
438
439 if (p->cid)
440 p->cid = modify_fft_labels ((VECTOR *) EAT (p->cid));
441 }
442 else /* Take the FFT of each column */
443 {
444 nc = p->nc; /* Number of columns */
445 nn = p->nr; /* Number to rows (elements per column) */
446
447 if (nn == 0)
448 goto out; /* Done if rows have zero length */
449
450 if (nc)
451 {
452 /* Allocate space - transform done in place */
453 tmpc = MALLOC (2 * nn * sizeof(REAL));
454
455 /* Create the fftw plan */
456 plan = fftw_plan_dft_1d (nn, (fftw_complex *) tmpc,
457 (fftw_complex *) tmpc, FFTW_FORWARD,
458 rigor);
459
460 if (plan == NULL) {
461 fail ("FFTW failed to make a plan.");
462 raise_exception ();
463 }
464
465 for (l = 0; l < nc; l++)
466 {
467 /* Copy the data */
468 for (i=0; i<nn; i++) {
469 tmpc[i].real = p->a.complex[i+l*nn].real;
470 tmpc[i].imag = p->a.complex[i+l*nn].imag;
471 }
472
473 /* Execute the transform */
474 fftw_execute(plan);
475
476 /* copy with reorder */
477 for (i=0,j=(nn-1)/2+l*nn; i<=nn/2; i++,j++) {
478 p->a.complex[j].real = tmpc[i].real;
479 p->a.complex[j].imag = tmpc[i].imag;
480 }
481 for (i=(nn+2)/2,j=l*nn; i<nn; i++,j++) {
482 p->a.complex[j].real = tmpc[i].real;
483 p->a.complex[j].imag = tmpc[i].imag;
484 }
485
486 } /* l - done loop over each column */
487
488 FREE (EAT (tmpc));
489 }
490
491 if (p->rid)
492 p->rid = modify_fft_labels ((VECTOR *) EAT (p->rid));
493
494 }
495 break;
496
497 default:
498 fail ("Can't apply \"fft\" to a %s matrix.", type_string[p->type]);
499 raise_exception ();
500 }
501 }
502
503 out: /* can't jump past the ON_EXCEPTION */
504
505 ON_EXCEPTION
506 {
507 if (plan) fftw_destroy_plan (plan);
508 TFREE (tmpc);
509 TFREE (tmpd);
510 delete_matrix (p);
511 }
512 END_EXCEPTION;
513
514 fftw_destroy_plan (plan);
515 return ENT (p);
516 }
517
hermcpy(COMPLEX * v,int x)518 void hermcpy ( COMPLEX *v, int x )
519 {
520 /* Hermetian copy and frequency shift converts:
521 * [( 1, 0),( 2, 2),( 3, 3),( 4, 0),( x, x),( x, x)] to
522 * [( 3,-3),( 2,-2),( 1, 0),( 2, 2),( 3, 3),( 4, 0)]
523 * or
524 * [( 1, 0),( 2, 2),( 3, 3),( 4, 4),( x, x),( x, x), ( x, x)] to
525 * [( 4,-4),( 3,-3),( 2,-2),( 1, 0),( 2, 2),( 3, 3), ( 4, 4)]
526 */
527
528 int i, xx, xm;
529
530 xx = x/2;
531
532 for (i=xx, xm=x-1; i>=0; i--,xm--)
533 {
534 v[xm].real = v[i].real;
535 v[xm].imag = v[i].imag;
536 }
537 for (i=0,xm=x-1-!(x%2); i<(x-1)/2; i++,xm--)
538 {
539 v[i].real = v[xm].real;
540 v[i].imag = -v[xm].imag;
541 }
542 }
543 #endif
544
545 #if !HAVE_FFTW || USE_BOTH_FFT
546
547 ENTITY *
bi_fft(n,p,d)548 bi_fft (n, p, d)
549 int n;
550 ENTITY *p;
551 ENTITY *d;
552 {
553 int dim = 1;
554 ENTITY * volatile s = NULL;
555
556 EASSERT (p, 0, 0);
557
558 if (d != NULL && d->class != table)
559 {
560 fail ("Invalid options -- not a table.");
561 delete_2_entities (p, d);
562 raise_exception ();
563 }
564
565 if (!OK_OPTIONS ((TABLE *) d, fft_options)) raise_exception ();
566
567 dim = !IN_TABLE (d, fft_options[FFT_OPT_ROW]);
568
569 WITH_HANDLING
570 {
571 switch (p->class)
572 {
573 case scalar:
574 s = fft_vector ((VECTOR *) scalar_to_vector ((SCALAR *) EAT(p)));
575 break;
576
577 case vector:
578 s = fft_vector ((VECTOR *) EAT(p));
579 break;
580
581 case matrix:
582
583 /* FFT of rows (dim=0); FFT of cols (dim=1) */
584 s = fft_matrix ((MATRIX *) EAT(p), dim);
585 break;
586
587 default:
588 BAD_CLASS (p->class);
589 raise_exception ();
590 }
591 }
592 ON_EXCEPTION
593 {
594 delete_2_entities (p, d);
595 }
596 END_EXCEPTION;
597
598 delete_entity (d);
599 return s;
600 }
601
602 ENTITY *
fft_vector(p)603 fft_vector (p)
604 VECTOR *p;
605 {
606 REAL *tmp;
607 REAL s[4];
608 int i, j, k, nn;
609
610 EASSERT (p, vector, 0);
611
612 WITH_HANDLING
613 {
614
615 p = (VECTOR *) dense_vector ((VECTOR *) EAT (p));
616 p = (VECTOR *) dup_vector (EAT (p));
617
618 switch (p->type)
619 {
620 case integer:
621 p = (VECTOR *) cast_vector_integer_real (p); /* fall through */
622 case real:
623
624 /* The RFFTF routine gives kind of screwy results (sine and */
625 /* cosine coefficients) which can't be used for convolution */
626 /* in the usual way. Let's just convert to complex instead. */
627
628 p = (VECTOR *) cast_vector_real_complex ((VECTOR *) EAT (p));
629
630 nn = p->nn;
631 if (nn == 0)
632 goto out; /* done if zero length */
633
634 tmp = (REAL *) CALLOC (4 * nn + 15, sizeof (REAL));
635
636 CFFTI (&nn, tmp);
637 CFFTF (&nn, p->a.real, tmp);
638 FREE (tmp);
639
640 /* Reorder for increasing frequency from -1/2t to +1/2t */
641
642 j = 2 * ((nn - 1) / 2);
643 k = (nn % 2) ? 2 : 4;
644 for (i = j; i < j + k; i++)
645 s[i - j] = p->a.real[i];
646 for (i = 0; i < j; i++)
647 {
648 p->a.real[2 * nn - j + i - k] = p->a.real[i];
649 p->a.real[i] = p->a.real[2 * nn - j + i];
650 }
651 for (i = 2 * nn - k; i < 2 * nn; i++)
652 p->a.real[i] = s[i - 2 * nn + k];
653 break;
654
655 case complex:
656
657 nn = p->nn;
658 if (nn == 0)
659 goto out; /* done if zero length */
660
661 tmp = (REAL *) CALLOC (4 * nn + 15, sizeof (REAL));
662
663 CFFTI (&nn, tmp);
664 CFFTF (&nn, p->a.real, tmp);
665 FREE (tmp);
666
667 /* Reorder for increasing frequency from -1/2t to +1/2t */
668
669 j = 2 * ((nn - 1) / 2);
670 k = (nn % 2) ? 2 : 4;
671 for (i = j; i < j + k; i++)
672 s[i - j] = p->a.real[i];
673 for (i = 0; i < j; i++)
674 {
675 p->a.real[2 * nn - j + i - k] = p->a.real[i];
676 p->a.real[i] = p->a.real[2 * nn - j + i];
677 }
678 for (i = 2 * nn - k; i < 2 * nn; i++)
679 p->a.real[i] = s[i - 2 * nn + k];
680 break;
681
682 default:
683 fail ("Can't apply \"fft\" to a %s vector.",
684 type_string[p->type]);
685 raise_exception ();
686 }
687
688 if (p->eid)
689 p->eid = modify_fft_labels ((VECTOR *) EAT (p->eid));
690
691 }
692
693 out: /* can't jump past the ON_EXCEPTION */
694
695 ON_EXCEPTION
696 {
697 delete_vector (p);
698 }
699 END_EXCEPTION;
700
701 return (ENT (p));
702 }
703
704 ENTITY *
fft_matrix(p,dim)705 fft_matrix (p, dim)
706 MATRIX *p;
707 int dim;
708 {
709 REAL *tmp;
710 REAL s[4];
711 int i, j, k, m, nn;
712
713 EASSERT (p, matrix, 0);
714
715 /* FFT of rows (dim==0) or cols (dim==1) of matrix. */
716
717 WITH_HANDLING
718 {
719
720 p = (MATRIX *) dense_matrix ((MATRIX *) EAT (p));
721 p = (MATRIX *) dup_matrix (EAT (p));
722 p->symmetry = general;
723
724 switch (p->type)
725 {
726 case integer:
727 p = (MATRIX *) cast_matrix_integer_real (p); /* fall through */
728
729 case real:
730
731 /* The RFFTF routine gives kind of screwy results (sine and */
732 /* cosine coefficients) which can't be used for convolution */
733 /* in the usual way. Let's just convert to complex instead. */
734
735 p = (MATRIX *) cast_matrix_real_complex ((MATRIX *) EAT (p));
736
737 nn = dim ? p->nr : p->nc;
738 if (nn == 0)
739 goto out; /* done if zero length */
740
741 tmp = (REAL *) CALLOC (4 * nn + 15, sizeof (REAL));
742
743 /* CFFTF has no stride option, so must transpose if doing rows. */
744
745 if (!dim)
746 p = (MATRIX *) transpose_matrix ((MATRIX *) EAT (p));
747
748 CFFTI (&nn, tmp);
749
750 j = 2 * ((nn - 1) / 2);
751 k = (nn % 2) ? 2 : 4;
752
753 for (m=0; m<2*p->nc; m+=2)
754 {
755 CFFTF (&nn, p->a.real+m*nn, tmp);
756
757 /* Reorder for increasing frequency from -1/2t to +1/2t */
758
759 for (i = j; i < j + k; i++)
760 s[i - j] = p->a.real[i+m*nn];
761 for (i = 0; i < j; i++)
762 {
763 p->a.real[(2+m) * nn - j + i - k] = p->a.real[i+m*nn];
764 p->a.real[i+m*nn] = p->a.real[(2+m) * nn - j + i];
765 }
766 for (i = 2 * nn - k; i < 2 * nn; i++)
767 p->a.real[i+m*nn] = s[i - 2 * nn + k];
768 }
769
770 FREE (tmp);
771
772 if (!dim)
773 p = (MATRIX *) transpose_matrix ((MATRIX *) EAT (p));
774
775 break;
776
777 case complex:
778
779 nn = dim ? p->nr : p->nc;
780 if (nn == 0)
781 goto out; /* done if zero length */
782
783 tmp = (REAL *) CALLOC (4 * nn + 15, sizeof (REAL));
784
785 /* CFFTF has no stride option, so must transpose if doing rows. */
786
787 if (!dim)
788 p = (MATRIX *) transpose_matrix ((MATRIX *) EAT (p));
789
790 CFFTI (&nn, tmp);
791
792 j = 2 * ((nn - 1) / 2);
793 k = (nn % 2) ? 2 : 4;
794
795 for (m=0; m<2*p->nc; m+=2)
796 {
797 CFFTF (&nn, p->a.real+m*nn, tmp);
798
799 /* Reorder for increasing frequency from -1/2t to +1/2t */
800
801 for (i = j; i < j + k; i++)
802 s[i - j] = p->a.real[i+m*nn];
803 for (i = 0; i < j; i++)
804 {
805 p->a.real[(2+m) * nn - j + i - k] = p->a.real[i+m*nn];
806 p->a.real[i+m*nn] = p->a.real[(2+m) * nn - j + i];
807 }
808 for (i = 2 * nn - k; i < 2 * nn; i++)
809 p->a.real[i+m*nn] = s[i - 2 * nn + k];
810 }
811
812 FREE (tmp);
813
814 if (!dim)
815 p = (MATRIX *) transpose_matrix ((MATRIX *) EAT (p));
816
817 break;
818
819 default:
820 fail ("Can't apply \"fft\" to a %s matrix.", type_string[p->type]);
821 raise_exception ();
822 }
823
824 if (dim)
825 {
826 if (p->rid)
827 p->rid = modify_fft_labels ((VECTOR *) EAT (p->rid));
828 }
829 else
830 {
831 if (p->cid)
832 p->cid = modify_fft_labels ((VECTOR *) EAT (p->cid));
833 }
834
835 }
836
837 out: /* can't jump past the ON_EXCEPTION */
838
839 ON_EXCEPTION
840 {
841 delete_matrix (p);
842 }
843 END_EXCEPTION;
844
845 return ENT (p);
846 }
847 #endif
848
849 #define SAMPLE_WARN 1.0e-6
850 #define SAMPLE_ERR 1.0e-3
851
852 static ENTITY *
modify_fft_labels(p)853 modify_fft_labels (p)
854 VECTOR *p;
855 {
856 /*
857 * The given vector contains abscissas corresponding to the ordinates
858 * in another vector that is to be transformed with the FFT. The result
859 * of this routine is a vector of frequency-domain abscissas.
860 */
861
862 REAL dt, dif, t, d, f, df;
863 int n, i;
864
865 switch (p->type)
866 {
867 case integer:
868 p = (VECTOR *) cast_vector_integer_real (p);
869 break;
870 case real:
871 break;
872 case complex:
873 case character:
874 delete_vector (p);
875 return (NULL);
876 default:
877 BAD_TYPE (p->type);
878 delete_vector (p);
879 raise_exception ();
880 }
881 assert (p->type == real);
882
883 p = (VECTOR *) dup_vector ((VECTOR *) dense_vector ((VECTOR *) EAT (p)));
884
885 WITH_HANDLING
886 {
887 n = p->nn;
888 if (n > 1)
889 {
890 t = p->a.real[n - 1] - p->a.real[0];
891 dt = t / (REAL) (n - 1);
892
893 /* Go through the samples to make sure that they are evenly spaced. */
894
895 dif = 0;
896 for (i = 1; i < n; i++)
897 {
898 d = fabs (p->a.real[i] - p->a.real[i - 1] - dt);
899 if (d > dif)
900 dif = d;
901 }
902 dif /= dt;
903 if (dif > SAMPLE_ERR)
904 {
905 fail ("Unevenly spaced samples in \"fft\".");
906 raise_exception ();
907 }
908 else if (dif > SAMPLE_WARN)
909 {
910 warn ("Unevenly spaced samples in \"fft\".");
911 }
912
913 /* Now compute the new "frequency" labels. */
914
915 f = 4.0 * atan (1.0) / dt;
916 df = 2.0 * f / (REAL) n;
917 for (i = 0; i < n; i++)
918 p->a.real[n - i - 1] = df * ((n / 2) - i);
919 }
920 else
921 if (n == 1) p->a.real[0] = 0.0;
922 }
923 ON_EXCEPTION
924 {
925 delete_vector (p);
926 }
927 END_EXCEPTION;
928
929 return (ENT (p));
930 }
931