promotional, or sales literature without prior written consent from
MIT in each case. */
-/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/array.c,v 9.35 1989/03/29 04:34:06 pas Exp $ */
-
-/* ARRAY =
- sequence of REAL(float or double numbers) with a tag on the front */
-/*___________________________________________________________________*/
-/* contents
- Scheme_Array datatype
- constructors, selectors, operators
- procedures for converting between C_Array, and Scheme_Vector
- basic and advanced array primitive operations (see manual.scm) */
-/*___________________________________________________________________*/
-
-/* See array.h for definition using NM_VECTOR,
- and for many useful EXTERN */
+/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/array.c,v 9.36 1989/06/22 21:50:18 pas Exp $ */
\f
#include "scheme.h"
#include <math.h>
#include <values.h>
/* <values.h> contains some math constants */
-/*__________________________________________________________*/
+/* .
+ ARRAY (as a scheme object)
+ is a usual array (in C) containing REAL numbers (float/double)
+ and tagged as a NM_VECTOR
+
+ Basic contents:
+ constructors, selectors, arithmetic operations,
+ conversion routines between C_Array, and Scheme_Vector
+
+ see array.h for macros, NM_VECTOR, and extern
+ */
-/* First some utilities */
+
+/* first some utilities */
int Scheme_Number_To_REAL(Arg, Cell) Pointer Arg; REAL *Cell;
/* 0 means conversion ok, 1 means too big, 2 means not a number */
}
return (0);
}
-\f
+
int Scheme_Number_To_Double(Arg, Cell) Pointer Arg; double *Cell;
/* 0 means conversion ok, 1 means too big, 2 means not a number */
{ long Value;
return (0);
}
-/*__________________begin__________________*/
+/* c */
-/* I think this is not needed, can be done at s-code ...
-DEFINE_PRIMITIVE ("ARRAY?", Prim_array_predicate, 1, 1, 0)
-{ Primitive_1_Args();
- if (Type_Code(Arg1)==TC_ARRAY) return SHARP_F;
- else return SHARP_F;
-}
-*/
+/* I think this is not needed, it can be done in scheme
+ DEFINE_PRIMITIVE ("ARRAY?", Prim_array_predicate, 1, 1, 0)
+ { Primitive_1_Args();
+ if (Type_Code(Arg1)==TC_ARRAY) return SHARP_F;
+ else return SHARP_F;
+ }
+ */
DEFINE_PRIMITIVE ("VECTOR->ARRAY", Prim_vector_to_array, 1, 1, 0)
{ Pointer Scheme_Vector_To_Scheme_Array();
return Scheme_Array_To_Scheme_Vector(Arg1);
}
-DEFINE_PRIMITIVE ("ARRAY-CONS", Prim_array_cons, 2, 2, 0)
-{ long Length, i, allocated_cells;
- REAL Init_Value, *Next;
- int Error_Number;
+/* array-cons = (array-allocate followed by array-initialize!)
+ The two are separated because all too often, we only need
+ array memory space. Even though the initialization is fast, it
+ happens so often that we get savings.
+ Also array-initialize! occurs via subarray-offset-scale!.
+
+ */
+DEFINE_PRIMITIVE ("ARRAY-ALLOCATE", Prim_array_allocate, 1,1, 0)
+{ long n,allocated_cells;
Pointer Result;
- Primitive_2_Args();
- Arg_1_Type(TC_FIXNUM);
- Range_Check(Length, Arg1, 0, ARRAY_MAX_LENGTH, ERR_ARG_1_BAD_RANGE);
+ PRIMITIVE_HEADER (1);
+ CHECK_ARG (1, FIXNUM_P);
- Error_Number = Scheme_Number_To_REAL(Arg2, &Init_Value);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ n = arg_nonnegative_integer(1); /* length of array to allocate */
+ if (n > ARRAY_MAX_LENGTH) error_bad_range_arg(1); /* avoid memory overflow */
- Allocate_Array(Result,Length,allocated_cells);
- Next = Scheme_Array_To_C_Array(Result);
- for (i=0; i<Length; i++) *Next++ = Init_Value;
- return Result;
+ Allocate_Array(Result,n,allocated_cells);
+ PRIMITIVE_RETURN (Result);
}
+
DEFINE_PRIMITIVE ("ARRAY-CONS-REALS", Prim_array_cons_reals, 3, 3, 0)
{ long i, Length, allocated_cells;
REAL *a;
double from, dt;
Pointer Result;
- int Error_Number;
+ int errcode;
Primitive_3_Args();
- Error_Number = Scheme_Number_To_Double(Arg1, &from); /* starting time */
- if (Error_Number == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- Error_Number = Scheme_Number_To_Double(Arg2, &dt); /* dt interval */
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ errcode = Scheme_Number_To_Double(Arg1, &from); /* starting time */
+ if (errcode == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
+ errcode = Scheme_Number_To_Double(Arg2, &dt); /* dt interval */
+ if (errcode == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
Length = Get_Integer(Arg3); /* number of points */
Allocate_Array(Result,Length,allocated_cells);
DEFINE_PRIMITIVE ("ARRAY-SET!", Prim_array_set, 3, 3, 0)
{ long Index;
REAL *Array, Old_Value;
- int Error_Number;
+ int errcode;
Primitive_3_Args();
Arg_1_Type(TC_ARRAY);
Array = Scheme_Array_To_C_Array(Arg1);
Old_Value = Array[Index];
- Error_Number = Scheme_Number_To_REAL(Arg3, &Array[Index]);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
+ errcode = Scheme_Number_To_REAL(Arg3, &Array[Index]);
+ if (errcode == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
Reduced_Flonum_Result((double) Old_Value);
}
-DEFINE_PRIMITIVE ("ARRAY-COPY", Prim_array_copy, 1, 1, 0)
-{ long Length, i, allocated_cells;
- REAL *To_Array, *From_Array;
- SCHEME_ARRAY Result;
-
- Primitive_1_Args();
- Arg_1_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
-
- Allocate_Array(Result, Length, allocated_cells);
- Array_Copy(Arg1, Result);
- return Result;
-}
-
-/*____________________data file readers___________
+/*____________________ file readers ___________
ascii and 2bint formats
- ________________________________________________*/
+ ______________________________________________*/
+
+/* Reading data from files
+ ATTENTION: for reading REAL numbers, use "lf" for double, "%f" for float
+ */
+#if (REAL_IS_DEFINED_DOUBLE == 1)
+#define REALREAD "%lf"
+#define REALREAD2 "%lf %lf"
+#else
+#define REALREAD "%f"
+#define REALREAD2 "%f %f"
+#endif
DEFINE_PRIMITIVE ("ARRAY-READ-ASCII-FILE", Prim_array_read_ascii_file, 2, 2, 0)
{ FILE *fp;
REAL *a; long N; FILE *fp;
{ long i;
for (i=0; i<N; i++) {
- if ( (fscanf(fp, "%lf", &(a[i]))) != 1)
+ if ( (fscanf(fp, REALREAD, &(a[i]))) != 1)
{ printf("Not enough values read ---\n Last Point was %d with value % .16e \n", i, a[i-1]);
return SHARP_F; }}
Close_File(fp);
Close_File(fp);
}
/* C_Array_Write_2bint_File
- is not implemented yet, don't have the time to to it now. */
+ is not implemented yet, I do not have the time to do it now. */
-\f
-DEFINE_PRIMITIVE ("SUBARRAY", Prim_subarray, 3, 3, 0)
-{ long Length, i, allocated_cells, Start, End, New_Length;
- REAL *To_Here, *From_Here;
- Pointer Result;
+/* ----- Read data from files --- end*/
- Primitive_3_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_FIXNUM);
- Arg_3_Type(TC_FIXNUM);
- Length = Array_Length(Arg1);
- Range_Check(Start, Arg2, 0, Array_Length(Arg1)-1, ERR_ARG_2_BAD_RANGE);
- Range_Check(End, Arg3, 0, Array_Length(Arg1)-1, ERR_ARG_3_BAD_RANGE);
- if (Start>End) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- New_Length = (End - Start) + 1;
- Allocate_Array(Result, New_Length, allocated_cells);
- From_Here = Nth_Array_Loc(Arg1, Start);
- To_Here = Scheme_Array_To_C_Array(Result);
+\f
+/* ARRAY-COPY! a very powerful primitive
+ See array.scm for its many applications.
+ Be Careful when source and destination are the same array.
+ */
+DEFINE_PRIMITIVE ("SUBARRAY-COPY!", Prim_subarray_copy, 5, 5, 0)
+{ long i, i1, i2;
+ long m, at1, at2;
+ REAL *a,*b;
- C_Array_Copy(From_Here, To_Here, New_Length);
- return Result;
-}
-
-DEFINE_PRIMITIVE ("ARRAY-SET-SUBARRAY!", Prim_array_set_subarray, 4, 4, 0)
-{ long Length, i, Start, End, New_Length;
- REAL *To_Here, *From_Here;
- Pointer Result;
-
- Primitive_4_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_FIXNUM);
- Arg_3_Type(TC_FIXNUM);
- Arg_4_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
- Range_Check(Start, Arg2, 0, Array_Length(Arg1)-1, ERR_ARG_2_BAD_RANGE);
- Range_Check(End, Arg3, 0, Array_Length(Arg1)-1, ERR_ARG_3_BAD_RANGE);
- if (Start>End) Primitive_Error(ERR_ARG_3_BAD_RANGE);
-
- New_Length = (End - Start) + 1;
- if (New_Length!=Array_Length(Arg4)) Primitive_Error(ERR_ARG_4_BAD_RANGE);
- From_Here = Scheme_Array_To_C_Array(Arg4);
- To_Here = Nth_Array_Loc(Arg1, Start);
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, ARRAY_P); /* source array a */
+ CHECK_ARG (2, ARRAY_P); /* destination array b */
- C_Array_Copy(From_Here, To_Here, New_Length);
- return Arg1;
-}
-
-DEFINE_PRIMITIVE ("ARRAY-APPEND", Prim_array_append, 2, 2, 0)
-{ long Length, Length1, Length2, i, allocated_cells;
- REAL *To_Here, *From_Here;
- Pointer Result;
-
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Length1 = Array_Length(Arg1);
- Length2 = Array_Length(Arg2);
- Length = Length1 + Length2;
-
- Allocate_Array(Result, Length, allocated_cells);
- To_Here = Scheme_Array_To_C_Array(Result);
- From_Here = Scheme_Array_To_C_Array(Arg1);
-
- for (i=0; i < Length1; i++) {
- *To_Here++ = *From_Here;
- From_Here++ ;
- }
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ b = Scheme_Array_To_C_Array(ARG_REF(2));
+ at1 = arg_nonnegative_integer(3); /* at1 = starting index in source array */
+ at2 = arg_nonnegative_integer(4); /* at2 = starting index in destination array */
+ m = arg_nonnegative_integer(5); /* m = number of points to copy */
+
+ if ((at1 + m) > (Array_Length(ARG_REF(1)))) error_bad_range_arg(3);
+ if ((at2 + m) > (Array_Length(ARG_REF(2)))) error_bad_range_arg(4);
+ /* These 2 checks cover all cases */
- From_Here = Scheme_Array_To_C_Array(Arg2);
- for (i=0; i < Length2; i++) {
- *To_Here++ = *From_Here;
- From_Here++ ;
- }
+ for (i=0,i1=at1,i2=at2; i<m; i++,i1++,i2++)
+ b[i2] = a[i1];
- return Result;
+ PRIMITIVE_RETURN (NIL);
}
+
DEFINE_PRIMITIVE ("ARRAY-REVERSE!", Prim_array_reverse, 1, 1, 0)
{ long Length, i,j, Half_Length;
REAL *Array, Temp;
return Arg1;
}
-DEFINE_PRIMITIVE ("ARRAY-SCALE!", Prim_array_scale, 2, 2, 0)
-{ long Length, i;
- REAL *To_Here, *From_Here, Scale;
- Pointer Result;
- int Error_Number;
+DEFINE_PRIMITIVE ("ARRAY-TIME-REVERSE!",
+ Prim_array_time_reverse, 1, 1, 0)
+{ long i, n;
+ REAL *a;
+ void C_Array_Time_Reverse();
+ PRIMITIVE_HEADER (1);
+ CHECK_ARG (1, ARRAY_P);
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ n = Array_Length(ARG_REF(1));
+
+ C_Array_Time_Reverse(a,n);
+
+ PRIMITIVE_RETURN (NIL);
+}
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
- Error_Number = Scheme_Number_To_REAL(Arg2, &Scale);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+/* time-reverse
+ x[0] remains fixed. (time=0)
+ x[i] swapped with x[n-i]. (mirror image around x[0])
+ */
+void C_Array_Time_Reverse(x,n)
+ REAL *x;
+ long n;
+{ long i, ni, n2;
+ REAL xt;
+ if ((n % 2) == 0) /* even length */
+ { n2 = (n/2);
+ for (i=1; i<n2; i++) /* i=1,2,..,n/2-1 */
+ { ni = n-i;
+ xt = x[i];
+ x[i] = x[ni];
+ x[ni] = xt; }}
+ else /* odd length */
+ { n2 = (n+1)/2; /* (n+1)/2 = (n-1)/2 + 1 */
+ for (i=1; i<n2; i++) /* i=1,2,..,(n-1)/2 */
+ { ni = n-i;
+ xt = x[i];
+ x[i] = x[ni];
+ x[ni] = xt; }}
+}
+
+/* The following is smart
+ and avoids computation when offset or scale are degenerate 0,1
+ */
+DEFINE_PRIMITIVE ("SUBARRAY-OFFSET-SCALE!",
+ Prim_subarray_offset_scale, 5, 5, 0)
+{ long i, at, m,mplus;
+ REAL *a, offset,scale;
+ int errcode;
+
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, FIXNUM_P);
+ CHECK_ARG (3, FIXNUM_P);
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ at = arg_nonnegative_integer(2); /* at = starting index */
+ m = arg_nonnegative_integer(3); /* m = number of points to change */
+
+ mplus = at + m;
+ if (mplus > (Array_Length(ARG_REF(1)))) error_bad_range_arg(3);
- Result = Arg1;
- From_Here = Scheme_Array_To_C_Array(Arg1);
- To_Here = Scheme_Array_To_C_Array(Result);
+ errcode = Scheme_Number_To_REAL(ARG_REF(4), &offset);
+ if (errcode==1) error_bad_range_arg(4); if (errcode==2) error_wrong_type_arg(4);
+ errcode = Scheme_Number_To_REAL(ARG_REF(5), &scale);
+ if (errcode==1) error_bad_range_arg(5); if (errcode==2) error_wrong_type_arg(5);
+
+ if ((offset == 0.0) && (scale == 1.0))
+ ; /* be smart */
+ else if (scale == 0.0)
+ for (i=at; i<mplus; i++) a[i] = offset;
+ else if (offset == 0.0)
+ for (i=at; i<mplus; i++) a[i] = scale * a[i];
+ else if (scale == 1.0)
+ for (i=at; i<mplus; i++) a[i] = offset + a[i];
+ else
+ for (i=at; i<mplus; i++) a[i] = offset + scale * a[i];
+
+ PRIMITIVE_RETURN (NIL);
+}
- for (i=0; i < Length; i++) {
- *To_Here++ = (Scale * (*From_Here));
- From_Here++ ;
- }
- return Result;
+
+DEFINE_PRIMITIVE ("COMPLEX-SUBARRAY-COMPLEX-SCALE!",
+ Prim_complex_subarray_complex_scale, 6,6, 0)
+{ long i, at,m,mplus;
+ REAL *a,*b; /* (a,b) = (real,imag) arrays */
+ double temp, minus_y, x, y; /* (x,y) = (real,imag) scale */
+ int errcode;
+
+ PRIMITIVE_HEADER (6);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, ARRAY_P);
+ CHECK_ARG (3, FIXNUM_P);
+ CHECK_ARG (4, FIXNUM_P);
+
+ at = arg_nonnegative_integer(3); /* at = starting index */
+ m = arg_nonnegative_integer(4); /* m = number of points to change */
+ mplus = at + m;
+ if (mplus > (Array_Length(ARG_REF(1)))) error_bad_range_arg(4);
+
+ errcode = Scheme_Number_To_Double(ARG_REF(5), &x);
+ if (errcode==1) error_bad_range_arg(5); if (errcode==2) error_wrong_type_arg(5);
+ errcode = Scheme_Number_To_Double(ARG_REF(6), &y);
+ if (errcode==1) error_bad_range_arg(6); if (errcode==2) error_wrong_type_arg(6);
+
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ b = Scheme_Array_To_C_Array(ARG_REF(2));
+ if ((Array_Length(ARG_REF(1))) != (Array_Length(ARG_REF(2)))) error_bad_range_arg(2);
+
+ if (x==0.0) /* imaginary only */
+ if (y==0.0)
+ for (i=at; i<mplus; i++)
+ { a[i] = 0.0;
+ b[i] = 0.0; }
+ else if (y==1.0)
+ for (i=at; i<mplus; i++)
+ { temp = b[i];
+ b[i] = a[i];
+ a[i] = (-temp); }
+ else if (y==-1.0)
+ for (i=at; i<mplus; i++)
+ { temp = b[i];
+ b[i] = (-a[i]);
+ a[i] = temp; }
+ else
+ { minus_y = (-y);
+ for (i=at; i<mplus; i++)
+ { temp = y * ((double) a[i]);
+ a[i] = (REAL) (minus_y * ((double) b[i]));
+ b[i] = (REAL) temp; }}
+ else if (y==0.0) /* real only */
+ if (x==1.0) ;
+ else for (i=at; i<mplus; i++)
+ { a[i] = (REAL) (x * ((double) a[i]));
+ b[i] = (REAL) (x * ((double) b[i])); }
+ else /* full complex scale */
+ for (i=at; i<mplus; i++)
+ { temp = ((double) a[i])*x - ((double) b[i])*y;
+ b[i] = (REAL) (((double) b[i])*x + ((double) a[i])*y);
+ a[i] = (REAL) temp; }
+
+ PRIMITIVE_RETURN (NIL);
}
+
+DEFINE_PRIMITIVE ("CS-ARRAY-TO-COMPLEX-ARRAY!",
+ Prim_cs_array_to_complex_array, 3, 3, 0)
+{ long n,n2,n2_1, i;
+ REAL *a, *b,*c;
+
+ PRIMITIVE_HEADER (3);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, ARRAY_P);
+ CHECK_ARG (3, ARRAY_P);
+
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ n = Array_Length(ARG_REF(1));
+ b = Scheme_Array_To_C_Array(ARG_REF(2));
+ c = Scheme_Array_To_C_Array(ARG_REF(3));
+ if (n!=(Array_Length(ARG_REF(2)))) error_bad_range_arg(2);
+ if (n!=(Array_Length(ARG_REF(3)))) error_bad_range_arg(3);
+
+ b[0] = a[0]; c[0] = 0.0;
+
+ n2 = n/2; /* integer division truncates down */
+ n2_1 = n2+1;
+
+ if (2*n2 == n) /* even length, n2 is only real */
+ { b[n2] = a[n2]; c[n2] = 0.0; }
+ else /* odd length, make the loop include the n2 index */
+ { n2 = n2+1;
+ n2_1 = n2; }
+
+ for (i=1; i<n2; i++) { b[i] = a[i];
+ c[i] = a[n-i]; }
+ for (i=n2_1; i<n; i++) { b[i] = a[n-i];
+ c[i] = (-a[i]); }
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+DEFINE_PRIMITIVE ("CS-ARRAY-MULTIPLY-INTO-SECOND-ONE!",
+ Prim_cs_array_multiply_into_second_one, 2, 2, 0)
+{ long n,n2;
+ REAL *a, *b;
+ void cs_array_multiply_into_second_one();
+ PRIMITIVE_HEADER (2);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, ARRAY_P);
+
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ n = Array_Length(ARG_REF(1));
+ b = Scheme_Array_To_C_Array(ARG_REF(2));
+ if (n!=(Array_Length(ARG_REF(2)))) error_bad_range_arg(2);
+ n2 = n/2; /* integer division truncates down */
+ cs_array_multiply_into_second_one(a,b, n,n2);
+ PRIMITIVE_RETURN (NIL);
+}
+
+void cs_array_multiply_into_second_one(a,b, n,n2)
+ REAL *a, *b; long n,n2;
+{ REAL temp;
+ long i,ni;
+ b[0] = a[0] * b[0];
+
+ if (2*n2 == n) /* even length, n2 is only real */
+ b[n2] = a[n2] * b[n2];
+ else
+ n2 = n2+1; /* odd length, make the loop include the n2 index */
+
+ for (i=1; i<n2; i++)
+ { ni = n-i;
+ temp = a[i]*b[i] - a[ni]*b[ni]; /* real part */
+ b[ni] = a[i]*b[ni] + a[ni]*b[i]; /* imag part */
+ b[i] = temp; }
+}
+
+DEFINE_PRIMITIVE ("CS-ARRAY-DIVIDE-INTO-XXX!",
+ Prim_cs_array_divide_into_xxx, 4, 4, 0)
+{ long n,n2, one_or_two;
+ REAL *a, *b, inf;
+ int errcode;
+ void cs_array_divide_into_z();
+
+ PRIMITIVE_HEADER (4);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, ARRAY_P);
+ errcode = Scheme_Number_To_REAL(ARG_REF(3), &inf);
+ if (errcode==1) error_bad_range_arg(3); if (errcode==2) error_wrong_type_arg(3);
+ CHECK_ARG (4, FIXNUM_P);
+ one_or_two = arg_nonnegative_integer(4); /* where to store result of division */
+
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ b = Scheme_Array_To_C_Array(ARG_REF(2));
+ n = Array_Length(ARG_REF(1));
+ if (n!=(Array_Length(ARG_REF(2)))) error_bad_range_arg(2);
+ n2 = n/2; /* integer division truncates down */
+
+ if (one_or_two == 1)
+ cs_array_divide_into_z(a,b, a, n,n2, inf);
+ else if (one_or_two == 2)
+ cs_array_divide_into_z(a,b, b, n,n2, inf);
+ else
+ error_bad_range_arg(4);
+ PRIMITIVE_RETURN (NIL);
+}
+
+void cs_array_divide_into_second_one(a,b, n,n2,inf) /* used in image.c */
+ REAL *a,*b, inf; long n,n2;
+{ void cs_array_divide_into_z();
+ cs_array_divide_into_z(a,b, b, n,n2,inf);
+}
+
+void cs_array_divide_into_z(a,b, z, n,n2, inf) /* z can be either a or b */
+ REAL *a,*b,*z, inf; long n,n2;
+{ long i,ni;
+ REAL temp, radius;
+
+ if (b[0] == 0.0)
+ if (a[0] == 0.0) z[0] = 1.0;
+ else z[0] = a[0] * inf;
+ else z[0] = a[0] / b[0];
+
+ if (2*n2 == n) /* even length, n2 is only real */
+ if (b[n2] == 0.0)
+ if (a[n2] == 0.0) z[n2] = 1.0;
+ else z[n2] = a[n2] * inf;
+ else z[n2] = a[n2] / b[n2];
+ else
+ n2 = n2+1; /* odd length, make the loop include the n2 index */
+
+ for (i=1; i<n2; i++)
+ { ni = n-i;
+ radius = b[i]*b[i] + b[ni]*b[ni]; /* b^2 denominator = real^2 + imag^2 */
+
+ if (radius == 0.0) {
+ if (a[i] == 0.0) z[i] = 1.0;
+ else z[i] = a[i] * inf;
+ if (a[ni] == 0.0) z[ni] = 1.0;
+ else z[ni] = a[ni] * inf; }
+ else {
+ temp = a[i]*b[i] + a[ni]*b[ni];
+ z[ni] = (a[ni]*b[i] - a[i]*b[ni]) / radius; /* imag part */
+ z[i] = temp / radius; /* real part */
+ }}
+}
+
+
+
+
/* ARRAY-UNARY-FUNCTION!
- applies a unary function on each element of an array.
- The name of this proc comes from "(array-unary-function! array function)"
+ apply unary-function elementwise on array
+
+ Available functions :
*/
-/* available functions
- */
-
void REALabs(a,b) REAL *a,*b;
{ (*b) = ( (REAL) fabs( (double) (*a)) );
}
}
void REALround(a,b) REAL *a,*b; /* towards nearest integer */
{ double integral_part, modf();
- if ((*a) >= 0.0) /* It may be faster to look at the sign of mantissa and dispatch */
+ if ((*a) >= 0.0) /* It may be faster to look at the sign
+ of mantissa, and dispatch */
modf( ((double) ((*a)+0.5)), &integral_part);
else
modf( ((double) ((*a)-0.5)), &integral_part);
(*b) = ( (REAL) yn(((int) order), ((double) (*a))) );
}
-/* Table to store the available functions for transforming arrays.
- It also stores the corresponding numofargs (whether unary or binary function).
+/* Table to store the available unary-functions.
+ Also some binary functions at the end -- not available right now.
+ The (1 and 2)s denote the numofargs (1 for unary 2 for binary)
*/
struct array_func_table {
return Result;
}
-/* The following is accumulate of + and *
- code numbers are 0 1
+
+/* Accumulate
+ using combinators + or *
+ corresponding type codes 0 1
*/
-DEFINE_PRIMITIVE ("ARRAY-ACCUMULATE", Prim_array_accumulate, 2, 2, 0)
-{ long Length, i;
- REAL *a, result;
- long functc;
+DEFINE_PRIMITIVE ("SUBARRAY-ACCUMULATE", Prim_subarray_accumulate, 4,4, 0)
+{ long at,m,mplus, tc, i;
+ REAL *a;
+ double result;
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_FIXNUM);
- Length = Array_Length(Arg1);
- Range_Check(functc, Arg2, 0, 1, ERR_ARG_2_BAD_RANGE);
- a = Scheme_Array_To_C_Array(Arg1);
+ PRIMITIVE_HEADER (4);
+ CHECK_ARG (1, ARRAY_P); /* a = input array */
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ tc = arg_nonnegative_integer(2); /* tc = type code 0 or 1 */
+ at = arg_nonnegative_integer(3); /* at = starting index */
+ m = arg_nonnegative_integer(4); /* m = number of points to process */
- if (functc==0)
+ mplus = at + m;
+ if (mplus > (Array_Length(ARG_REF(1)))) error_bad_range_arg(4);
+
+ if (tc==0)
{ result = 0.0;
- for (i=0;i<Length;i++) result = result + a[i];
- }
- else if (functc==1)
+ for (i=at;i<mplus;i++) result = result + ((double) a[i]); }
+ else if (tc==1)
{ result = 1.0;
- for (i=0;i<Length;i++) result = result * a[i];
- }
- Reduced_Flonum_Result((double) result);
+ for (i=at;i<mplus;i++) result = result * ((double) a[i]); }
+ else
+ error_bad_range_arg(2);
+
+ Reduced_Flonum_Result(result);
}
+
/* The following searches for value within tolerance
starting from index=from in array.
Returns first index where match occurs. -- (useful for finding zeros)
{ long Length, from, i;
REAL *a, value; /* value to search for */
double tolerance; /* tolerance allowed */
- int Error_Number;
+ int errcode;
Primitive_4_Args();
Arg_1_Type(TC_ARRAY);
a = Scheme_Array_To_C_Array(Arg1); Length = Array_Length(Arg1);
- Error_Number = Scheme_Number_To_REAL(Arg2, &value);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
- Error_Number = Scheme_Number_To_Double(Arg3, &tolerance);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
+ errcode = Scheme_Number_To_REAL(Arg2, &value);
+ if (errcode == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ errcode = Scheme_Number_To_Double(Arg3, &tolerance);
+ if (errcode == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
Arg_4_Type(TC_FIXNUM);
Range_Check(from, Arg4, 0, Length-1, ERR_ARG_4_BAD_RANGE);
return SHARP_F;
}
-DEFINE_PRIMITIVE ("ARRAY-MIN-MAX-INDEX", Prim_array_min_max_index, 1, 1, 0)
-{ long Length, nmin, nmax;
+DEFINE_PRIMITIVE ("SUBARRAY-MIN-MAX-INDEX", Prim_subarray_min_max_index, 3,3, 0)
+{ long at,m,mplus;
+ long nmin, nmax;
Pointer Result, *Orig_Free;
- REAL *Array;
-
- Primitive_1_Args();
- Arg_1_Type(TC_ARRAY);
- Array= Scheme_Array_To_C_Array(Arg1);
- Length = Array_Length(Arg1);
- C_Array_Find_Min_Max(Array, Length, &nmin, &nmax);
+ REAL *a;
+
+ PRIMITIVE_HEADER (3);
+ CHECK_ARG (1, ARRAY_P);
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ at = arg_nonnegative_integer(2); /* at = starting index */
+ m = arg_nonnegative_integer(3); /* m = number of points to process */
+
+ mplus = at + m;
+ if (mplus > (Array_Length(ARG_REF(1)))) error_bad_range_arg(3);
+
+ C_Array_Find_Min_Max ( &(a[at]), m, &nmin, &nmax);
+ nmin = nmin + at; /* offset appropriately */
+ nmax = nmax + at;
+
Primitive_GC_If_Needed(4);
Result = Make_Pointer(TC_LIST, Free);
Orig_Free = Free;
*Orig_Free++ = Make_Pointer(TC_LIST, Orig_Free+1);
*Orig_Free++ = Make_Non_Pointer(TC_FIXNUM, nmax);
*Orig_Free=EMPTY_LIST;
- return Result;
+
+ PRIMITIVE_RETURN (Result);
}
+
void C_Array_Find_Min_Max(x, n, nmin, nmax) REAL *x; long n, *nmax, *nmin;
{ REAL *xold = x;
register REAL xmin, xmax;
}
-/* The following becomes obsolete.
- Done using array-reduce + divide by array-length
+/* array-average
+ can be done with (array-reduce +) and division by array-length.
+ But there is also this C primitive. Keep it around, may be useful someday.
*/
+
DEFINE_PRIMITIVE ("ARRAY-AVERAGE", Prim_array_find_average, 1, 1, 0)
{ long Length; REAL average;
+ void C_Array_Find_Average();
Primitive_1_Args();
Arg_1_Type(TC_ARRAY);
Length = Array_Length(Arg1);
C_Array_Find_Average( Scheme_Array_To_C_Array(Arg1), Length, &average);
Reduced_Flonum_Result((double) average);
}
+
/* Computes the average in pieces, so as to reduce
roundoff smearing in cumulative sum.
example= first huge positive numbers, then small nums, then huge negative numbers.
*/
+
void C_Array_Find_Average(Array, Length, pAverage)
long Length; REAL *Array, *pAverage;
{ long i;
}
DEFINE_PRIMITIVE ("ARRAY-CLIP-MIN-MAX!", Prim_array_clip_min_max, 3, 3, 0)
-{ long Length, i; /* , allocated_cells; */
+{ long Length, i;
REAL *To_Here, *From_Here, xmin, xmax;
- Pointer Result;
- int Error_Number;
+ int errcode;
Primitive_3_Args();
Arg_1_Type(TC_ARRAY);
- Error_Number=Scheme_Number_To_REAL(Arg2, &xmin);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
- Error_Number=Scheme_Number_To_REAL(Arg3, &xmax);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
+ errcode=Scheme_Number_To_REAL(Arg2, &xmin);
+ if (errcode == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ errcode=Scheme_Number_To_REAL(Arg3, &xmax);
+ if (errcode == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
Length = Array_Length(Arg1);
- Result = Arg1;
+
From_Here = Scheme_Array_To_C_Array(Arg1);
- To_Here = Scheme_Array_To_C_Array(Result);
+ To_Here = Scheme_Array_To_C_Array(Arg1);
if (xmin>xmax) Primitive_Error(ERR_ARG_3_BAD_RANGE);
for (i=0; i < Length; i++) {
else if ((*From_Here)>xmax) *To_Here++ = xmax;
else *To_Here++ = *From_Here;
From_Here++ ; }
- return Result;
+ return SHARP_F;
}
-DEFINE_PRIMITIVE ("ARRAY-MAKE-POLAR!", Prim_array_make_polar, 2, 2, 0)
+DEFINE_PRIMITIVE ("COMPLEX-ARRAY-TO-POLAR!",
+ Prim_complex_array_to_polar, 2, 2, 0)
{ long Length, i;
REAL *To_Here_Mag, *To_Here_Phase;
REAL *From_Here_Real, *From_Here_Imag;
- Pointer Result_Mag, Result_Phase, answer;
-
+ Pointer Result_Mag, Result_Phase;
+
Primitive_2_Args();
Arg_1_Type(TC_ARRAY);
Arg_2_Type(TC_ARRAY);
To_Here_Phase = Scheme_Array_To_C_Array(Result_Phase);
for (i=0; i < Length; i++)
- {
- C_Make_Polar(*From_Here_Real, *From_Here_Imag, *To_Here_Mag, *To_Here_Phase);
+ { C_Make_Polar(*From_Here_Real, *From_Here_Imag, *To_Here_Mag, *To_Here_Phase);
From_Here_Real++ ; From_Here_Imag++ ;
- To_Here_Mag++ ; To_Here_Phase++ ;
- }
+ To_Here_Mag++ ; To_Here_Phase++ ; }
- Primitive_GC_If_Needed(4);
- answer = Make_Pointer(TC_LIST, Free);
- *Free++ = Result_Mag;
- *Free = Make_Pointer(TC_LIST, Free+1);
- Free += 1;
- *Free++ = Result_Phase;
- *Free++ = EMPTY_LIST;
- return answer;
+ return SHARP_F;
}
-DEFINE_PRIMITIVE ("ARRAY-FIND-MAGNITUDE", Prim_array_find_magnitude, 2, 2, 0)
-{ long Length, i, allocated_cells;
- REAL *From_Here_Real, *From_Here_Imag, *To_Here;
- Pointer Result;
+DEFINE_PRIMITIVE ("COMPLEX-ARRAY-MAGNITUDE!", Prim_complex_array_magnitude, 3, 3, 0)
+{ long n, i;
+ REAL *a, *b, *mg;
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
- if (Length != Array_Length(Arg2)) Primitive_Error(ERR_ARG_1_BAD_RANGE);
+ PRIMITIVE_HEADER (3);
+ CHECK_ARG (1, ARRAY_P); /* input array -- n real part */
+ CHECK_ARG (2, ARRAY_P); /* input array -- n imag part */
+ CHECK_ARG (3, ARRAY_P); /* ouput array -- n magnitude */
+
+ n = Array_Length(ARG_REF(1));
+ if (n != Array_Length(ARG_REF(2))) error_bad_range_arg(2);
+ if (n != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+
+ a = Scheme_Array_To_C_Array(ARG_REF(1)); /* real part */
+ b = Scheme_Array_To_C_Array(ARG_REF(2)); /* imag part */
+ mg = Scheme_Array_To_C_Array(ARG_REF(3)); /* magnitude */
+
+ for (i=0; i<n; i++)
+ mg[i] = (REAL) sqrt( (double) a[i]*a[i] + b[i]*b[i] );
+
+ PRIMITIVE_RETURN (NIL);
+}
- Allocate_Array(Result, Length, allocated_cells);
- To_Here = Scheme_Array_To_C_Array(Result);
- From_Here_Real = Scheme_Array_To_C_Array(Arg1);
- From_Here_Imag = Scheme_Array_To_C_Array(Arg2);
- for (i=0; i<Length; i++) {
- C_Find_Magnitude(*From_Here_Real, *From_Here_Imag, *To_Here);
- From_Here_Real++ ;
- From_Here_Imag++ ;
- To_Here++ ;
+DEFINE_PRIMITIVE ("CS-ARRAY-MAGNITUDE!", Prim_cs_array_magnitude, 1, 1, 0)
+{ long n, i;
+ REAL *a;
+ void cs_array_magnitude();
+ PRIMITIVE_HEADER (1);
+ CHECK_ARG (1, ARRAY_P);
+ a = Scheme_Array_To_C_Array(ARG_REF(1)); /* input cs-array */
+ n = Array_Length(ARG_REF(1)); /* becomes a standard array on return */
+
+ cs_array_magnitude(a,n);
+ PRIMITIVE_RETURN (NIL);
+}
+
+/* result is a standard array (even signal, real data)
+ */
+void cs_array_magnitude(a,n)
+ REAL *a; long n;
+{ long i, n2, ni;
+ n2 = n/2; /* integer division truncates down */
+
+ a[0] = (REAL) fabs((double) a[0]); /* imag=0 */
+
+ if (2*n2 == n) /* even length, n2 is only real */
+ a[n2] = (REAL) fabs((double) a[n2]); /* imag=0 */
+ else
+ n2 = n2+1; /* odd length, make the loop include the n2 index */
+
+ for (i=1; i<n2; i++)
+ { ni = n-i;
+ a[i] = (REAL) sqrt( (double) a[i]*a[i] + (double) a[ni]*a[ni] );
+ a[ni] = a[i]; /* even signal */
}
- return Result;
}
-/* ATTENTION: To1,To2 SHOULD BE Length1-1, and Length2-2 RESPECTIVELY ! */
+/* Rectangular and Polar cs-arrays
+
+ cs-arrays have Even magnitude and (almost) Odd angle
+ hence we store the magnitude into the real part of cs-array
+ and the angle into the imag part of cs-array
+
+ Except for index 0 and index n2(when n even) are only real
+ Hence the angle can be either 0 or pi
+ ---> information is lost when going from rect to polar (compromise to save space)
+
+ To invert from polar to rect, assume the signal is a continuous curve
+ and choose sign (from angle 0 or pi) at index i=0 same sign as i=1,
+ and i=n2 same sign as i=n2+1
+ */
+
+DEFINE_PRIMITIVE ("CS-ARRAY-TO-POLAR!", Prim_cs_array_to_polar, 1,1, 0)
+{ long n, i;
+ REAL *a;
+ void cs_array_to_polar();
+
+ PRIMITIVE_HEADER (1);
+ CHECK_ARG (1, ARRAY_P); /* input and output array -- both cs-arrays */
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ n = Array_Length(ARG_REF(1));
+
+ cs_array_to_polar(a,n);
+ PRIMITIVE_RETURN (NIL);
+}
+
+void cs_array_to_polar(a,n)
+ REAL *a; long n;
+{ long i, n2;
+ double real, imag; /* temporary variables */
+ n2 = n/2; /* integer division truncates down */
+
+ a[0] = (REAL) fabs((double) a[0]); /* implicitly angle = 0, but it could be pi */
+ if (2*n2 == n) /* even length, n2 is only real */
+ a[n2] = (REAL) fabs((double) a[n2]); /* implicitly angle = 0, but it could be pi */
+ else
+ n2 = n2+1; /* odd length, make the loop include the n2 index */
+
+ for (i=1; i<n2; i++)
+ { real = (double) a[i];
+ imag = (double) a[n-i];
+ a[i] = (REAL) sqrt( real*real + imag*imag );
+ if (a[i] == 0.0)
+ a[n-i] = 0.0;
+ else
+ a[n-i] = (REAL) atan2( imag, real ); }
+}
+
+DEFINE_PRIMITIVE ("CS-ARRAY-TO-RECTANGULAR!",
+ Prim_cs_array_to_rectangular, 1,1, 0)
+{ long n,n2, i;
+ double magn,angl; /* temporary variables */
+ REAL *a;
+
+ PRIMITIVE_HEADER (1);
+ CHECK_ARG (1, ARRAY_P); /* input and output array -- both cs-arrays */
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ n = Array_Length(ARG_REF(1));
+ n2 = n/2; /* integer division truncates down */
+
+ if (a[1] > 0.0)
+ a[0] = a[0]; /* assume angle = 0 */
+ else
+ a[0] = (-a[0]); /* assume angle = pi */
+
+ if (2*n2 == n) /* even length, n2 is real only */
+ if (a[n2+1] > 0.0)
+ a[n2] = a[n2]; /* assume angle = 0 */
+ else
+ a[n2] = (-a[n2]); /* assume angle = pi */
+ else
+ n2 = n2+1; /* odd length, make the loop include the n2 index */
+
+ for (i=1; i<n2; i++)
+ { magn = (double) a[i];
+ angl = (double) a[n-i];
+ a[i] = (REAL) magn * cos(angl);
+ a[n-i] = (REAL) magn * sin(angl); }
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+
+/* Convolution in the Time-Domain
+ */
+
+/* In the following macro
+ To1 and To2 should be (Length1-1) and (Length2-1) respectively.
+ */
#define C_Convolution_Point_Macro(X, Y, To1, To2, N, Result) \
{ long Min_of_N_To1=min((N),(To1)); \
long mi, N_minus_mi; \
Reduced_Flonum_Result(C_Result);
}
-DEFINE_PRIMITIVE ("ARRAY-CONVOLUTION", Prim_array_convolution, 2, 2, 0)
-{ long Endpoint1, Endpoint2, allocated_cells, i;
- /* ASSUME A SIGNAL FROM INDEX 0 TO ENDPOINT=LENGTH-1 */
- long Resulting_Length;
- REAL *Array1, *Array2, *To_Here;
- Pointer Result;
+DEFINE_PRIMITIVE ("ARRAY-CONVOLUTION-IN-TIME!",
+ Prim_array_convolution_in_time, 3, 3, 0)
+{ long n,m,l, n_1,m_1, i;
+ REAL *a,*b,*c;
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Endpoint1 = Array_Length(Arg1) - 1;
- Endpoint2 = Array_Length(Arg2) - 1;
- Resulting_Length = Endpoint1 + Endpoint2 + 1;
- Array1 = Scheme_Array_To_C_Array(Arg1);
- Array2 = Scheme_Array_To_C_Array(Arg2);
-
- allocated_cells = (Resulting_Length * REAL_SIZE) + ARRAY_HEADER_SIZE;
- Primitive_GC_If_Needed(allocated_cells);
- Result = Make_Pointer(TC_ARRAY, Free);
- Free[ARRAY_HEADER] = Make_Non_Pointer(TC_MANIFEST_ARRAY, allocated_cells-1);
- Free[ARRAY_LENGTH] = Resulting_Length;
- Free += allocated_cells;
- To_Here = Scheme_Array_To_C_Array(Result);
+ PRIMITIVE_HEADER (3);
+ CHECK_ARG (1, ARRAY_P); /* input array a -- length n */
+ CHECK_ARG (2, ARRAY_P); /* input array b -- length m */
+ CHECK_ARG (3, ARRAY_P); /* ouput array c -- length l = (n + m - 1) */
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ b = Scheme_Array_To_C_Array(ARG_REF(2));
+ c = Scheme_Array_To_C_Array(ARG_REF(3));
- for (i=0; i<Resulting_Length; i++) {
- C_Convolution_Point_Macro(Array1, Array2, Endpoint1, Endpoint2, i, *To_Here);
- To_Here++;
- }
- return Result;
+ n = Array_Length(ARG_REF(1));
+ m = Array_Length(ARG_REF(2));
+ l = n+m-1; /* resulting length */
+ if (l != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+
+ n_1 = n-1; m_1 = m-1;
+ for (i=0; i<l; i++)
+ { C_Convolution_Point_Macro(a, b, n_1, m_1, i, c[i]); }
+
+ PRIMITIVE_RETURN (NIL);
}
-DEFINE_PRIMITIVE ("ARRAY-MULTIPLICATION-INTO-SECOND-ONE!", Prim_array_multiplication_into_second_one, 2, 2, 0)
+DEFINE_PRIMITIVE ("ARRAY-MULTIPLY-INTO-SECOND-ONE!",
+ Prim_array_multiply_into_second_one, 2, 2, 0)
{ long Length, i;
REAL *To_Here;
REAL *From_Here_1, *From_Here_2;
return Result;
}
-DEFINE_PRIMITIVE ("ARRAY-COMPLEX-MULTIPLICATION-INTO-SECOND-ONE!", Prim_array_complex_multiplication_into_second_one, 4, 4, 0)
+DEFINE_PRIMITIVE ("COMPLEX-ARRAY-MULTIPLY-INTO-SECOND-ONE!",
+ Prim_complex_array_multiply_into_second_one, 4, 4, 0)
{ long Length, i;
REAL *To_Here_1, *To_Here_2;
REAL *From_Here_1, *From_Here_2, *From_Here_3, *From_Here_4;
}
return SHARP_F;
}
-void C_Array_Complex_Multiply_Into_First_One(a,b,c,d, length)
+
+void C_Array_Complex_Multiply_Into_First_One(a,b,c,d, length) /* used in fft.c */
REAL *a,*b,*c,*d; long length;
{ long i;
REAL temp;
}
-DEFINE_PRIMITIVE ("ARRAY-DIVISION-INTO-FIRST-ONE!", Prim_array_division_into_first_one, 3, 3, 0)
-{ long Length, i;
- SCHEME_ARRAY scheme_result;
- REAL *x,*y,*result;
- REAL infinity;
- int Error_Number;
+DEFINE_PRIMITIVE ("ARRAY-DIVIDE-INTO-XXX!",
+ Prim_array_divide_into_xxx, 4,4, 0)
+{ long n, i, one_or_two;
+ REAL *x,*y,*z, inf;
+ int errcode;
+ void array_divide_into_z();
- Primitive_3_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
- if (Length != Array_Length(Arg2)) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- Error_Number = Scheme_Number_To_REAL(Arg3, &infinity); /* User-Provided Infinity */
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
+ PRIMITIVE_HEADER (4);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, ARRAY_P);
+ errcode = Scheme_Number_To_REAL(ARG_REF(3), &inf);
+ if (errcode==1) error_bad_range_arg(3); if (errcode==2) error_wrong_type_arg(3);
+ CHECK_ARG (4, FIXNUM_P);
+ one_or_two = arg_nonnegative_integer(4); /* where to store result of division */
- scheme_result = Arg1;
- result = Scheme_Array_To_C_Array(scheme_result);
- x = Scheme_Array_To_C_Array(Arg1);
- y = Scheme_Array_To_C_Array(Arg2);
+ x = Scheme_Array_To_C_Array(ARG_REF(1));
+ y = Scheme_Array_To_C_Array(ARG_REF(2));
+ n = Array_Length(ARG_REF(1));
+ if (n!=(Array_Length(ARG_REF(2)))) error_bad_range_arg(2);
- for (i=0; i < Length; i++) {
- if (y[i] == 0.0) {
- if (x[i] == 0.0) /* zero/zero */
- result[i] = 1.0;
- else
- result[i] = infinity * x[i];
- }
- else
- result[i] = x[i] / y[i];
- }
- return scheme_result;
+ if (one_or_two == 1)
+ array_divide_into_z( x,y, x, n, inf);
+ else if (one_or_two == 2)
+ array_divide_into_z( x,y, y, n, inf);
+ else
+ error_bad_range_arg(4);
+ PRIMITIVE_RETURN (NIL);
}
-DEFINE_PRIMITIVE ("ARRAY-DIVISION-INTO-SECOND-ONE!", Prim_array_division_into_second_one, 3, 3, 0)
-{ long Length, i;
- SCHEME_ARRAY scheme_result;
- REAL *x,*y,*result;
- REAL infinity;
- int Error_Number;
-
- Primitive_3_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
- if (Length != Array_Length(Arg2)) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- Error_Number = Scheme_Number_To_REAL(Arg3, &infinity); /* User-Provided Infinity */
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
-
- scheme_result = Arg2;
- result = Scheme_Array_To_C_Array(scheme_result);
- x = Scheme_Array_To_C_Array(Arg1);
- y = Scheme_Array_To_C_Array(Arg2);
-
- for (i=0; i < Length; i++) {
+void array_divide_into_z( x,y, z, n, inf) /* z can either x or y */
+ REAL *x,*y,*z, inf; long n;
+{ long i;
+ for (i=0; i<n; i++) {
if (y[i] == 0.0) {
- if (x[i] == 0.0) /* zero/zero */
- result[i] = 1.0;
- else
- result[i] = infinity * x[i];
- }
- else
- result[i] = x[i] / y[i];
+ if (x[i] == 0.0) z[i] = 1.0;
+ else z[i] = inf * x[i]; }
+ else z[i] = x[i] / y[i];
}
- return scheme_result;
}
-DEFINE_PRIMITIVE ("ARRAY-COMPLEX-DIVISION-INTO-FIRST-ONE!", Prim_array_complex_multiplication_into_first_one, 5, 5, 0)
-{ long Length, i;
- SCHEME_ARRAY scheme_result_r, scheme_result_i;
- REAL *x_r,*x_i, *y_r,*y_i, *result_r,*result_i;
- register REAL Temp, radius;
- REAL infinity;
- int Error_Number;
+DEFINE_PRIMITIVE ("COMPLEX-ARRAY-DIVIDE-INTO-XXX!",
+ Prim_complex_array_divide_into_xxx, 6,6, 0)
+{ long n, i, one_or_two;
+ REAL inf, *xr,*xi, *yr,*yi;
+ void complex_array_divide_into_z();
+ int errcode;
- Primitive_5_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Arg_3_Type(TC_ARRAY);
- Arg_4_Type(TC_ARRAY);
- Length = Array_Length(Arg1);
- if (Length != Array_Length(Arg2)) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Length != Array_Length(Arg3)) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Length != Array_Length(Arg4)) Primitive_Error(ERR_ARG_4_BAD_RANGE);
- Error_Number = Scheme_Number_To_REAL(Arg5, &infinity); /* User-Provided Infinity */
- if (Error_Number == 1) Primitive_Error(ERR_ARG_5_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_5_WRONG_TYPE);
-
- scheme_result_r = Arg1;
- scheme_result_i = Arg2;
- result_r = Scheme_Array_To_C_Array(scheme_result_r);
- result_i = Scheme_Array_To_C_Array(scheme_result_i);
- x_r = Scheme_Array_To_C_Array(Arg1);
- x_i = Scheme_Array_To_C_Array(Arg2);
- y_r = Scheme_Array_To_C_Array(Arg3);
- y_i = Scheme_Array_To_C_Array(Arg4);
+ PRIMITIVE_HEADER (6);
+ CHECK_ARG (1, ARRAY_P); /* real numerator */
+ CHECK_ARG (2, ARRAY_P); /* imag numerator */
+ CHECK_ARG (3, ARRAY_P); /* real denominator */
+ CHECK_ARG (4, ARRAY_P); /* imag denominator */
+ CHECK_ARG (6, FIXNUM_P); /* one_or_two = where to store the result */
+ n = Array_Length(ARG_REF(1));
+ if (n != Array_Length(ARG_REF(2))) error_bad_range_arg(2);
+ if (n != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+ if (n != Array_Length(ARG_REF(4))) error_bad_range_arg(4);
+ errcode = Scheme_Number_To_REAL(ARG_REF(5), &inf); /* User-Provided Infinity */
+ if (errcode==1) error_bad_range_arg(5); if (errcode==2) error_wrong_type_arg(5);
- for (i=0; i < Length; i++) {
- Temp = (x_r[i] * y_r[i]) + (x_i[i] * y_i[i]);
- radius = (y_r[i] * y_r[i]) + (y_i[i] * y_i[i]);
-
+ one_or_two = arg_nonnegative_integer(6);
+ xr = Scheme_Array_To_C_Array(ARG_REF(1));
+ xi = Scheme_Array_To_C_Array(ARG_REF(2));
+ yr = Scheme_Array_To_C_Array(ARG_REF(3));
+ yi = Scheme_Array_To_C_Array(ARG_REF(4));
+
+ if (one_or_two == 1)
+ complex_array_divide_into_z(xr,xi, yr,yi, xr,xi, n, inf);
+ else if (one_or_two == 2)
+ complex_array_divide_into_z(xr,xi, yr,yi, yr,yi, n, inf);
+ else
+ error_bad_range_arg(6);
+ PRIMITIVE_RETURN (NIL);
+}
+
+void complex_array_divide_into_z(xr,xi, yr,yi, zr,zi, n, inf) /* z can be either x or y */
+ REAL *xr,*xi, *yr,*yi, *zr,*zi, inf; long n;
+{ long i;
+ register REAL temp, radius;
+
+ for (i=0; i<n; i++)
+ { radius = (yr[i] * yr[i]) + (yi[i] * yi[i]); /* denominator */
if (radius == 0.0) {
- if (x_r[i] == 0.0) result_r[i] = 1.0;
- else result_r[i] = infinity * x_r[i];
- if (x_i[i] == 0.0) result_i[i] = 1.0;
- else result_i[i] = infinity * x_i[i];
- }
+ if (xr[i] == 0.0) zr[i] = 1.0;
+ else zr[i] = inf * xr[i];
+ if (xi[i] == 0.0) zi[i] = 1.0;
+ else zi[i] = inf * xi[i]; }
else {
- result_i[i] = ( (x_i[i] * y_r[i]) - (x_r[i] * y_i[i]) ) / radius;
- result_r[i] = Temp / radius;
- }
- }
- return SHARP_F;
+ temp = xr[i] * yr[i] + xi[i] * yi[i];
+ zi[i] = (xi[i] * yr[i] - xr[i] * yi[i]) / radius;
+ zr[i] = temp / radius;
+ }}
}
-DEFINE_PRIMITIVE ("ARRAY-LINEAR-SUPERPOSITION-INTO-SECOND-ONE!", Prim_array_linear_superposition_into_second_one, 4, 4, 0)
-{ long Length, i;
+
+DEFINE_PRIMITIVE ("ARRAY-LINEAR-SUPERPOSITION-INTO-SECOND-ONE!",
+ Prim_array_linear_superposition_into_second_one, 4, 4, 0)
+{ long n, i;
REAL *To_Here, Coeff1, Coeff2;
REAL *From_Here_1, *From_Here_2;
Pointer Result;
- int Error_Number;
+ int errcode;
Primitive_4_Args();
- Error_Number = Scheme_Number_To_REAL(Arg1, &Coeff1);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
+ errcode = Scheme_Number_To_REAL(Arg1, &Coeff1);
+ if (errcode == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
Arg_2_Type(TC_ARRAY);
- Error_Number = Scheme_Number_To_REAL(Arg3, &Coeff2);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
+ errcode = Scheme_Number_To_REAL(Arg3, &Coeff2);
+ if (errcode == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
Arg_4_Type(TC_ARRAY);
- Length = Array_Length(Arg2);
- if (Length != Array_Length(Arg4)) Primitive_Error(ERR_ARG_4_BAD_RANGE);
+ n = Array_Length(Arg2);
+ if (n != Array_Length(Arg4)) Primitive_Error(ERR_ARG_4_BAD_RANGE);
Result = Arg4;
From_Here_2 = Scheme_Array_To_C_Array(Arg4);
To_Here = Scheme_Array_To_C_Array(Result);
- for (i=0; i < Length; i++) {
+ for (i=0; i < n; i++) {
*To_Here++ = (Coeff1 * (*From_Here_1)) + (Coeff2 * (*From_Here_2));
From_Here_1++ ;
From_Here_2++ ;
double twopi = 6.28318530717958;
Pointer Result, Pfunction_number, Psignal_frequency;
Pointer Pfunction_Number;
- int Error_Number;
+ int errcode;
REAL *To_Here;
double unit_square_wave(), unit_triangle_wave();
Arg_4_Type(TC_FIXNUM);
Range_Check(Function_Number, Arg1, 0, 10, ERR_ARG_1_BAD_RANGE); /* fix this */
- Error_Number = Scheme_Number_To_Double(Arg2, &Signal_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ errcode = Scheme_Number_To_Double(Arg2, &Signal_Frequency);
+ if (errcode == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
if (Signal_Frequency == 0) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- Error_Number = Scheme_Number_To_Double(Arg3, &Sampling_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
+ errcode = Scheme_Number_To_Double(Arg3, &Sampling_Frequency);
+ if (errcode == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
if (Sampling_Frequency == 0) Primitive_Error(ERR_ARG_3_BAD_RANGE);
Range_Check(N, Arg4, 0, ARRAY_MAX_LENGTH, ERR_ARG_4_BAD_RANGE);
double Sampling_Frequency, DT, DTi;
double twopi = 6.28318530717958;
Pointer Result;
- int Error_Number;
+ int errcode;
REAL *To_Here, twopi_dt;
Primitive_3_Args();
Arg_3_Type(TC_FIXNUM);
Range_Check(Function_Number, Arg1, 0, 6, ERR_ARG_1_BAD_RANGE);
- Error_Number = Scheme_Number_To_Double(Arg2, &Sampling_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ errcode = Scheme_Number_To_Double(Arg2, &Sampling_Frequency);
+ if (errcode == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
if (Sampling_Frequency == 0) Primitive_Error(ERR_ARG_2_BAD_RANGE);
Range_Check(N, Arg3, 0, ARRAY_MAX_LENGTH, ERR_ARG_3_BAD_RANGE);
{ Pointer *From_Here;
REAL *To_Here;
long Length, i;
- int Error_Number;
+ int errcode;
From_Here = Nth_Vector_Loc(Scheme_Vector, VECTOR_DATA);
To_Here = Array;
Length = Vector_Length(Scheme_Vector);
for (i=0; i < Length; i++, From_Here++) {
- Error_Number = Scheme_Number_To_REAL(*From_Here, To_Here);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
+ errcode = Scheme_Number_To_REAL(*From_Here, To_Here);
+ if (errcode == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
+ if (errcode == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
To_Here++; /* this gets incremented by REAL_SIZE ! */
}
promotional, or sales literature without prior written consent from
MIT in each case. */
-/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/array.h,v 9.29 1989/02/19 17:51:28 jinx Exp $ */
+/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/array.h,v 9.30 1989/06/22 21:51:57 pas Rel $ */
\f
#define REAL float
+#define REAL_IS_DEFINED_DOUBLE 0
+/*
+ When REAL is float, set = 0
+ When REAL is double, set = 1
+ This is used by #ifdef in some places like "fscanf"
+ */
+
#define REAL_SIZE ((sizeof(Pointer)+sizeof(REAL)-1)/ sizeof(Pointer))
-/* REAL should be either double or float.
- */
-/* Scheme_Arrays are implemented as NON_MARKED_VECTOR's.
- Do not forget to include object.h
- */
+/* Scheme_Arrays are implemented as NON_MARKED_VECTOR
+ Do not forget to include object.h */
+
+#define ARRAY_P NON_MARKED_VECTOR_P
+/* This is used in places like "CHECK_ARG(1, ARRAY_P)" */
#define TC_ARRAY TC_NON_MARKED_VECTOR
#define TC_MANIFEST_ARRAY TC_MANIFEST_NM_VECTOR
#define Linear_Map(slope,offset,From,To) { (To) = (((slope)*(From))+offset); }
-#define C_Find_Magnitude(Real, Imag, Mag_Cell) \
-{ double double_Real=((double) Real), double_Imag=((double) Imag); \
- Mag_Cell = (REAL) sqrt((double_Real*double_Real)+(double_Imag*double_Imag)); \
-}
-
#define mabs(x) (((x)<0) ? -(x) : (x))
#define max(x,y) (((x)<(y)) ? (y) : (x))
#define min(x,y) (((x)<(y)) ? (x) : (y))
-/* FROM ARRAY.C */
+/* From array.c
+ */
extern int Scheme_Number_To_REAL();
extern int Scheme_Number_To_Double();
-extern void C_Array_Find_Min_Max(); /* Find the index of the minimum (*nmin), maximum (*nmax). */
-extern void C_Array_Find_Average();
+extern void C_Array_Find_Min_Max();
extern void C_Array_Make_Histogram(); /* REAL *Array,*Histogram; long Length,npoints */
extern void C_Array_Complex_Multiply_Into_First_One();
/* Pointer Scheme_Vector; REAL *Array;
*/
-/* FROM BOB-XT.C */
+/* From bob-xt.c */
extern void Find_Offset_Scale_For_Linear_Map();
/* REAL Min,Max, New_Min,New_Max, *Offset,*Scale;
*/
promotional, or sales literature without prior written consent from
MIT in each case. */
-/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/fft.c,v 9.25 1989/02/19 18:09:08 jinx Exp $ */
+/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/fft.c,v 9.26 1989/06/22 21:52:02 pas Rel $ */
-/* Fourier Transforms (pas)
- 1. DFT (FFT),
- 2. CZT (chirp-z-transform).
- 3. 2D DFT
- */
+/* Time-Frequency Transforms (pas) */
#include "scheme.h"
#include "prims.h"
#include "array.h"
#include "image.h"
-/* The following implement the DFT as defined in Seibert's Book (6003).
- FORWARD DFT ---- Negative exponent and division by N
- BACKWARD DFT --- Positive exponent
- Note: Seibert's Forward DFT is Oppenheim's Backward DFT.
+/* SUMMARY
+ - pas_cft (complex data, DIF, split-radix)
+ - pas_rft (real data, DIT, split-radix) output is conjugate-symmetric
+ - pas_csft (cs data, DIF, split-radix) output is real
+ - pas_cft
+ - pas_rft2d
+ - pas_csft2d
+
+
+ Stuff before 4-15-1989
+ - C_Array_FFT (complex data, radix=2, NOT-in-place)
+ - CZT (chirp-z-transform) uses the old cft (hence slow).
+ - 2d DFT
+ */
+
+/* The DFT is as defined in Siebert 6003 book,
+ i.e.
+ forward DFT = Negative exponent and division by N
+ backward DFT = Positive exponent
+ (note Seibert forward DFT is Oppenheim backward DFT)
*/
+/* mathematical constants */
#ifdef PI
#undef PI
#endif
-
#define PI 3.141592653589793238462643
#define TWOPI 6.283185307179586476925287
+#define SQRT_2 1.4142135623730950488
+#define ONE_OVER_SQRT_2 .7071067811865475244
/* Abramowitz and Stegun */
+
+
+DEFINE_PRIMITIVE ("PAS-CFT!", Prim_pas_cft, 5, 5, 0)
+{ long i, length, power, flag;
+ REAL *f1,*f2, *wcos,*w3cos,*w3sin;
+ void pas_cft();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* flag forward-backward */
+ CHECK_ARG (2, ARRAY_P); /* real part */
+ CHECK_ARG (3, ARRAY_P); /* imag part */
+ CHECK_ARG (4, ARRAY_P); /* twiddle tables, total length = 3*(length/4) */
+ CHECK_ARG (5, FIXNUM_P); /* (1)=tables precomputed, else recompute */
+
+ flag = Get_Integer(ARG_REF(1));
+ length = Array_Length(ARG_REF(2));
+ if (length != (Array_Length(ARG_REF(3)))) error_bad_range_arg(2);
+
+ for (power=0, i=length; i>1; power++)
+ { if ( (i % 2) == 1) error_bad_range_arg(2);
+ i=i/2; }
+
+ f1 = Scheme_Array_To_C_Array(ARG_REF(2));
+ f2 = Scheme_Array_To_C_Array(ARG_REF(3));
+ if (f1==f2) error_wrong_type_arg(2);
+
+ wcos = Scheme_Array_To_C_Array(ARG_REF(4)); /* twiddle tables */
+ if (Array_Length(ARG_REF(4)) != (3*length/4)) error_bad_range_arg(4);
+ w3cos = wcos + length/4;
+ w3sin = w3cos + length/4;
+ if ((arg_nonnegative_integer(5)) == 1)
+ pas_cft(1, flag, f1,f2, length, power, wcos,w3cos,w3sin);
+ else
+ pas_cft(0, flag, f1,f2, length, power, wcos,w3cos,w3sin);
+ /* 1 means tables are already made
+ 0 means compute new tables */
+
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+DEFINE_PRIMITIVE ("PAS-CFT-MAKE-TWIDDLE-TABLES!",
+ Prim_pas_cft_make_twiddle_tables, 2, 2, 0)
+{ long length, power, i;
+ REAL *wcos,*w3cos,*w3sin;
+ void pas_cft_make_twiddle_tables_once();
+ PRIMITIVE_HEADER (2);
+
+ length = arg_nonnegative_integer(1); /* length of cft that we intend to compute */
+ CHECK_ARG (2, ARRAY_P); /* storage for twiddle tables */
+ if (Array_Length(ARG_REF(2)) != (3*length/4)) error_bad_range_arg(2);
+
+ power=0;
+ for (power=0, i=length; i>1; power++)
+ { if ( (i % 2) == 1) error_bad_range_arg(1);
+ i=i/2; }
+
+ wcos = Scheme_Array_To_C_Array(ARG_REF(2)); /* twiddle tables */
+ w3cos = wcos + length/4;
+ w3sin = w3cos + length/4;
+ pas_cft_make_twiddle_tables_once(length, power, wcos,w3cos,w3sin);
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+
+/*
+ C COMPLEX FOURIER TRANSFORM (Split-Radix, Decimation-in-frequency)
+ C (adapted and optimized from Sorensen,et.al. ASSP-34 no.1 page 152, February 1986)
+ */
+
+/* Twiddle Tables for PAS_CFT;
+ (tables for forward transform only)
+ Inverse transform === forward CFT (without 1/N scaling) followed by time-reversal.
+ /
+ The tables contain (2pi/N)*i for i=0,1,2,..,N/4
+ (except i=0 is ignored, never used)
+ /
+ Table for wsin[i] is not needed because wsin[i]=wcos[n4-i].
+ Table for w3sin[i] is needed however. The previous relationship does not work for w3sin.
+ */
+
+/* There are two routines for making twiddle tables:
+ a fast one, and a slower one but more precise.
+ The differences in speed and accuracy are actually rather small, but anyway.
+ Use the slow one for making permanent tables.
+ */
+
+void pas_cft_make_twiddle_tables(n,m, wcos,w3cos,w3sin) /* efficient version */
+ REAL *wcos, *w3cos, *w3sin;
+ long n,m;
+{ long i, n4;
+ double tm;
+ REAL costm,sintm;
+ n4 = n/4;
+ for (i=1; i<n4; i++) /* start from table entry 1 */
+ { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
+ wcos[i] = (REAL) cos(tm);
+ }
+ for (i=1; i<n4; i++)
+ { costm = wcos[i];
+ sintm = wcos[n4-i];
+ w3cos[i] = costm * (1 - 4*sintm*sintm); /* see my notes */
+ w3sin[i] = sintm * (4*costm*costm - 1);
+ }
+}
+
+void pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin) /* slow version, more accurate */
+ REAL *wcos, *w3cos, *w3sin;
+ long n,m;
+{ long i, n4;
+ double tm;
+ REAL costm,sintm;
+ n4 = n/4;
+ for (i=1; i<n4; i++) /* start from table entry 1 */
+ { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
+ wcos[i] = (REAL) cos(tm);
+ tm = tm * 3.0; /* this is more precise (in the 16th decimal) than */
+ w3cos[i] = (REAL) cos(tm); /* the more efficient version. (I tested by for/backward) */
+ w3sin[i] = (REAL) sin(tm);
+ }
+}
+
+void pas_cft(tables_ok,flag, x,y,n,m, wcos,w3cos,w3sin)
+ REAL *x,*y, *wcos,*w3cos,*w3sin;
+ long n,m, flag, tables_ok;
+{ REAL scale;
+ long i,j;
+ void pas_cft_make_twiddle_tables();
+ void C_Array_Time_Reverse();
+ void pas_cft_forward_loop();
+
+ if (tables_ok != 1) /* 1 means = tables already made */
+ pas_cft_make_twiddle_tables(n,m, wcos,w3cos,w3sin);
+
+ if (flag == 1) /* forward cft */
+ { pas_cft_forward_loop(x,y,n,m, wcos,w3cos,w3sin);
+ scale = (REAL) (1.0 / ((double) n));
+ for (i=0; i<n; i++)
+ { x[i] = x[i]*scale;
+ y[i] = y[i]*scale; }}
+ else /* backward cft */
+ { for (j=0; j<n; j++) y[j] = (-y[j]); /* conjugate before */
+ pas_cft_forward_loop(x,y,n,m, wcos,w3cos,w3sin);
+ for (j=0; j<n; j++) y[j] = (-y[j]); /* conjugate after */
+ }
+}
+
+void pas_cft_forward_loop(x,y,n,m, wcos,w3cos,w3sin) /* n >= 4 */
+ REAL *x,*y, *wcos,*w3cos,*w3sin;
+ long n,m;
+{ /* REAL a,a3,e; no need anymore, use tables */
+ REAL r1,r2,s1,s2,s3, xt, cc1,cc3,ss1,ss3;
+ long n1,n2,n4, i,j,k, is,id, i0,i1,i2,i3;
+ long windex0, windex, windex_n4; /* indices for twiddle tables */
+ /********** fortran indices start from 1,... **/
+ x = x-1; /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */
+ y = y-1;
+ /********** fortran indices start from 1,... **/
+ /* c */
+ /* c-----first M-1 stages of transform */
+ /* c */
+ windex_n4 = n/4; /* need for indexing sin via wcos twiddle table */
+ n2 = 2*n;
+ for (k=1; k<m; k++) /* DO 10 K = 1, M-1 */
+ { n2 = n2>>1; /* n2 = n2/2; */
+ n4 = n2>>2; /* n4 = n2/4; */
+ /* e = 6.283185307179586476925287 / ((REAL) n2); no need anymore, use tables */
+ /* a = 0.0; */
+ {
+ /* j=1; */
+ /* The first iteration in the loop "DO 20 J = 1, N4"
+ is done specially to save operations involving sin=0, cos=1 */
+ /* a = j*e; no need anymore, use tables */
+ is = 1; /* is = j; */
+ id = 2*n2;
+ label40first:
+ for (i0=is; i0<n; i0=i0+id) /* 40 DO 30 I0 = IS,N-1,ID */
+ { i1 = i0 + n4;
+ i2 = i1 + n4;
+ i3 = i2 + n4;
+ /* c */
+ r1 = x[i0] - x[i2];
+ x[i0] = x[i0] + x[i2];
+ r2 = x[i1] - x[i3];
+ x[i1] = x[i1] + x[i3];
+ s1 = y[i0] - y[i2];
+ y[i0] = y[i0] + y[i2];
+ s2 = y[i1] - y[i3];
+ y[i1] = y[i1] + y[i3];
+ /* c */
+ s3 = r1 - s2;
+ r1 = r1 + s2;
+ s2 = s1 - r2; /* original used to be s2 = r2 - s1; */
+ r2 = r2 + s1;
+ x[i2] = r1;
+ y[i2] = s2; /* used to be y[i2] = (-s2); */
+ x[i3] = s3;
+ y[i3] = r2;
+ /* x[i2] = r1*cc1 + s2*ss1; used to be, see below
+ y[i2] = s2*cc1 - r1*ss1; used to be, see below, inside the DO 20 J=1,N4
+ x[i3] = s3*cc3 + r2*ss3;
+ y[i3] = r2*cc3 - s3*ss3; */
+ } /* 30 CONTINUE */
+ is = 2*id - n2 + 1; /* is = 2*id - n2 + j; */
+ id = 4*id;
+ if (is < n) goto label40first; /* IF (IS.LT.N) GOTO 40 */
+ }
+ /* c */
+ windex0 = 1<<(k-1);
+ windex = windex0;
+ for (j=2; j<=n4; j++) /* DO 20 J = 1, N4 */
+ {
+ /* windex = (j-1)*(1<<(k-1)); -- done with trick to avoid (j-1) and 1<<(k-1) */
+ cc1 = wcos[windex];
+ ss1 = wcos[windex_n4 - windex]; /* see my notes */
+ cc3 = w3cos[windex];
+ ss3 = w3sin[windex]; /* sin-from-cos trick does not work here */
+ windex = j*windex0; /* same trick as "a = j*e" */
+ /* a3 = 3*a;
+ cc1 = cos(a);
+ ss1 = sin(a);
+ cc3 = cos(a3);
+ ss3 = sin(a3);
+ a = j*e;*/
+ is = j;
+ id = 2*n2;
+ label40:
+ for (i0=is; i0<n; i0=i0+id) /* 40 DO 30 I0 = IS,N-1,ID */
+ { i1 = i0 + n4;
+ i2 = i1 + n4;
+ i3 = i2 + n4;
+ /* c */
+ r1 = x[i0] - x[i2];
+ x[i0] = x[i0] + x[i2];
+ r2 = x[i1] - x[i3];
+ x[i1] = x[i1] + x[i3];
+ s1 = y[i0] - y[i2];
+ y[i0] = y[i0] + y[i2];
+ s2 = y[i1] - y[i3];
+ y[i1] = y[i1] + y[i3];
+ /* c */
+ s3 = r1 - s2;
+ r1 = r1 + s2;
+ s2 = s1 - r2; /* original used to be s2 = r2 - s1; */
+ r2 = r2 + s1;
+ x[i2] = r1*cc1 + s2*ss1; /* used to be x[i2] = r1*cc1 - s2*ss1; */
+ y[i2] = s2*cc1 - r1*ss1; /* used to be y[i2] = (-s2*cc1 - r1*ss1); */
+ x[i3] = s3*cc3 + r2*ss3;
+ y[i3] = r2*cc3 - s3*ss3;
+ } /* 30 CONTINUE */
+ is = 2*id - n2 + j;
+ id = 4*id;
+ if (is < n) goto label40; /* IF (IS.LT.N) GOTO 40 */
+ } /* 20 CONTINUE */
+ } /* 10 CONTINUE */
+ /* c
+ c-----------last-stage, length-2 butterfly ----------------c
+ c */
+ is = 1;
+ id = 4;
+ label50:
+ for (i0=is; i0<=n; i0=i0+id) /* 50 DO 60 I0 = IS, N, ID */
+ { i1 = i0 + 1;
+ r1 = x[i0];
+ x[i0] = r1 + x[i1];
+ x[i1] = r1 - x[i1];
+ r1 = y[i0];
+ y[i0] = r1 + y[i1];
+ y[i1] = r1 - y[i1];
+ } /* 60 CONTINUE */
+ is = 2*id - 1;
+ id = 4*id;
+ if (is < n) goto label50; /* IF (IS.LT.N) GOTO 50 */
+ /*
+ c
+ c-----------bit-reverse-counter---------------c
+ */
+ label100:
+ j = 1;
+ n1 = n - 1;
+ for (i=1; i<=n1; i++) /* DO 104 I = 1, N1 */
+ { if (i >= j) goto label101; /* if (i .ge. j) goto 101 */
+ xt = x[j];
+ x[j] = x[i];
+ x[i] = xt;
+ xt = y[j];
+ y[j] = y[i];
+ y[i] = xt;
+ label101: k = n>>1; /* k = n/2; */
+ label102: if (k>=j) goto label103;
+ j = j - k;
+ k = k>>1; /* k = k/2; */
+ goto label102;
+ label103: j = j + k;
+ } /* 104 CONTINUE */
+ /* c-------------------------------------*/
+ /* c */
+} /* RETURN \r END */
+
+
+\f
+DEFINE_PRIMITIVE ("PAS-RFT-CSFT!", Prim_pas_rft_csft, 5, 5, 0)
+{ long i, length, power, flag, ft_type;
+ REAL *f1, *wcos,*w3cos,*w3sin;
+ void pas_rft(), pas_csft();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* flag 1=forward, else backward transform */
+ CHECK_ARG (2, ARRAY_P); /* Input data (real or cs) */
+ CHECK_ARG (3, ARRAY_P); /* Twiddle tables, total length = 4*(length/8) */
+ CHECK_ARG (4, FIXNUM_P); /* (1)=tables precomputed, else recompute */
+ CHECK_ARG (5, FIXNUM_P); /* ft_type = 1 or 3
+ 1 means compute rft, 3 means compute csft */
+ flag = Get_Integer(ARG_REF(1));
+ f1 = Scheme_Array_To_C_Array(ARG_REF(2));
+ length = Array_Length(ARG_REF(2));
+ for (power=0, i=length; i>1; power++)
+ { if ( (i % 2) == 1) error_bad_range_arg(2);
+ i=i/2; }
+
+ wcos = Scheme_Array_To_C_Array(ARG_REF(3)); /* twiddle tables */
+ if (Array_Length(ARG_REF(3)) != (4*length/8)) error_bad_range_arg(3);
+ w3cos = wcos + (length/4);
+ w3sin = w3cos + (length/8);
+
+ ft_type = (arg_nonnegative_integer(5)); /* rft or csft */
+ if (ft_type == 1) {
+ if ((arg_nonnegative_integer(4)) == 1)
+ pas_rft (1, flag, f1, length, power, wcos,w3cos,w3sin);
+ else pas_rft (0, flag, f1, length, power, wcos,w3cos,w3sin);
+ }
+ else if (ft_type == 3) {
+ if ((arg_nonnegative_integer(4)) == 1)
+ pas_csft (1, flag, f1, length, power, wcos,w3cos,w3sin);
+ else pas_csft (0, flag, f1, length, power, wcos,w3cos,w3sin);
+ /* 1 means tables are already made
+ 0 means compute new tables */
+ }
+ else error_bad_range_arg(5);
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+DEFINE_PRIMITIVE ("PAS-REALDATA-MAKE-TWIDDLE-TABLES!",
+ Prim_pas_realdata_make_twiddle_tables, 2, 2, 0)
+{ long length, power, i;
+ REAL *wcos,*w3cos,*w3sin;
+ void pas_realdata_make_twiddle_tables_once();
+ PRIMITIVE_HEADER (2);
+
+ length = arg_nonnegative_integer(1); /* length of rft that we intend to compute */
+ CHECK_ARG (2, ARRAY_P); /* storage for twiddle tables */
+ if (Array_Length(ARG_REF(2)) != (4*length/8)) error_bad_range_arg(2);
+
+ power=0;
+ for (power=0, i=length; i>1; power++)
+ { if ( (i % 2) == 1) error_bad_range_arg(1);
+ i=i/2; }
+
+ wcos = Scheme_Array_To_C_Array(ARG_REF(2)); /* twiddle tables */
+ w3cos = wcos + length/4;
+ w3sin = w3cos + length/8;
+ pas_realdata_make_twiddle_tables_once(length, power, wcos,w3cos,w3sin);
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+/*
+ C REAL FOURIER TRANSFORM (Split-Radix, Decimation-in-time)
+ C (adapted from Sorensen,et.al. ASSP-35 no.6 page 849, October 1986)
+ C
+ C the output is [Re(0),Re(1),...,Re(n/2), Im(n/2-1),...,Im(1)]
+ */
+
+/* Twiddle Tables for PAS_RFT and PAS_CSFT
+ are identical. -> pas_realdata_make_twiddle_tables
+ (but they are indexed differently in each case)
+ /
+ The tables contain (2pi/N)*i where i=0,1,2,..,N/4 for wcos
+ and (2pi/N)*i where i=0,1,2,..,N/8 for w3cos w3sin
+ (except i=0 is ignored, never used)
+ /
+ Table for wsin[i] is not needed because wsin[i]=wcos[n4-i].
+ Table for w3sin[i] is needed however. The previous relationship does not work for w3sin.
+ /
+ Instead of getting sin() from a wsin[i] table with i=1,..,N/8
+ we get it from wcos[n4-i].
+ This way we can use a CFT table which goes up to N/4
+ for RFT CSFT also. We do so in image-processing (rft2d-csft2d).
+ */
+
+/* There are two routines for making twiddle tables:
+ a fast one, and a slower one but more precise.
+ The differences in speed and accuracy are actually rather small, but anyway.
+ Use the slow one for making tables that stay around.
+ */
+
+void pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin) /* efficient version */
+ REAL *wcos, *w3cos, *w3sin;
+ long n,m;
+{ long i, n4, n8;
+ double tm;
+ REAL costm,sintm;
+ n4 = n/4;
+ n8 = n/8;
+ for (i=1; i<n4; i++) /* start from table entry 1 */
+ { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
+ wcos[i] = (REAL) cos(tm);
+ }
+ for (i=1; i<n8; i++)
+ { costm = wcos[i];
+ sintm = wcos[n4-i];
+ w3cos[i] = costm * (1 - 4*sintm*sintm); /* see my notes */
+ w3sin[i] = sintm * (4*costm*costm - 1);
+ }
+}
+
+void pas_realdata_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin) /* slow version, more accurate */
+ REAL *wcos, *w3cos, *w3sin;
+ long n,m;
+{ long i, n4, n8;
+ double tm;
+ REAL costm,sintm;
+ n4 = n/4;
+ n8 = n/8;
+ for (i=1; i<n8; i++) /* start from table entry 1 */
+ { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
+ wcos[i] = (REAL) cos(tm);
+ tm = tm * 3.0; /* this is more precise (in the 16th decimal) than */
+ w3cos[i] = (REAL) cos(tm); /* the more efficient version. (I tested by for/backward) */
+ w3sin[i] = (REAL) sin(tm);
+ }
+ for (i=n8; i<n4; i++)
+ { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
+ wcos[i] = (REAL) cos(tm);
+ }
+}
+
+void pas_rft(tables_ok,flag, x,n,m, wcos,w3cos,w3sin)
+ REAL *x, *wcos,*w3cos,*w3sin;
+ long n,m, flag, tables_ok;
+{ REAL scale;
+ long i;
+ void pas_realdata_make_twiddle_tables();
+ void pas_rft_forward_loop();
+
+ if (tables_ok != 1) /* 1 means = tables already made */
+ pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin);
+
+ pas_rft_forward_loop(x,n,m, wcos,w3cos,w3sin);
+
+ if (flag == 1) /* forward rft */
+ { scale = (REAL) (1.0 / ((double) n));
+ for (i=0; i<n; i++) x[i] = x[i] * scale; }
+ else /* backward rft */
+ for (i=((n/2)+1); i<n; i++) x[i] = (-x[i]); /* time-reverse cs-array */
+}
+
+/* rft
+ forward transform === forward_loop + 1/N scaling
+ inverse transform === forward_loop + time-reversal (without 1/N scaling)
+ */
+
+/* wcos must be length n/4
+ w3cos, w3sin must be length n/8
+ (greater than n/8 is fine also, e.g. use cft tables)
+ */
+
+void pas_rft_forward_loop(x,n,m, wcos,w3cos,w3sin)
+ REAL *x, *wcos,*w3cos,*w3sin;
+ long n,m;
+{ /* REAL a,a3,e; no need anymore, use tables */
+ REAL r1, xt, cc1,cc3,ss1,ss3, t1,t2,t3,t4,t5,t6;
+ long n1,n2,n4,n8, i,j,k, is,id, i0,i1,i2,i3,i4,i5,i6,i7,i8;
+ long windex0, windex, windex_n4; /* indices for twiddle tables */
+ /********** fortran indices start from 1,... **/
+ x = x-1; /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */
+ /********** fortran indices start from 1,... **/
+ /* c */
+ windex_n4 = n/4; /* need for indexing sin via wcos twiddle table */
+ /* c
+ c-----------bit-reverse-counter---------------c
+ */
+ label100:
+ j = 1;
+ n1 = n - 1;
+ for (i=1; i<=n1; i++) /* DO 104 I = 1, N1 */
+ { if (i >= j) goto label101; /* if (i .ge. j) goto 101 */
+ xt = x[j];
+ x[j] = x[i];
+ x[i] = xt;
+ label101: k = n>>1; /* k = n/2; */
+ label102: if (k>=j) goto label103;
+ j = j - k;
+ k = k>>1; /* k = k/2; */
+ goto label102;
+ label103: j = j + k;
+ } /* 104 CONTINUE */
+ /* c-------------------------------------*/
+ /* c */
+ /* c ----length-two-butterflies----------- */
+ is = 1;
+ id = 4;
+ label70:
+ for (i0=is; i0<=n; i0=i0+id) /* 70 DO 60 I0 = IS,N,ID */
+ { i1 = i0 + 1;
+ r1 = x[i0];
+ x[i0] = r1 + x[i1];
+ x[i1] = r1 - x[i1];
+ } /* 60 CONTINUE */
+ is = 2*id - 1;
+ id = 4*id;
+ if (is < n) goto label70; /* IF (IS.LT.N) GOTO 70 */
+ /* C
+ C -------L-shaped-butterflies-------- */
+ n2 = 2;
+ for (k=2; k<=m; k++) /* DO 10 K = 2,M */
+ { n2 = n2 * 2;
+ n4 = n2>>2; /* n4 = n2/4; */
+ n8 = n2>>3; /* n8 = n2/8; */
+ /* e = 6.283185307179586476925287 / ((REAL) n2); no need anymore, use tables */
+ is = 0;
+ id = n2 * 2;
+ label40:
+ for (i=is; i<n; i=i+id) /* 40 DO 38 I = IS,N-1,ID */
+ { i1 = i + 1;
+ i2 = i1 + n4;
+ i3 = i2 + n4;
+ i4 = i3 + n4;
+ t1 = x[i4] + x[i3];
+ x[i4] = x[i4] - x[i3];
+ x[i3] = x[i1] - t1;
+ x[i1] = x[i1] + t1;
+ if (n4 == 1) goto label38; /* IF (N4.EQ.1) GOTO 38 */
+ i1 = i1 + n8;
+ i2 = i2 + n8;
+ i3 = i3 + n8;
+ i4 = i4 + n8;
+ /* t1 = (x[i3] + x[i4]) / sqrt(2.0); -- this is more precise, it uses extended
+ t2 = (x[i3] - x[i4]) / sqrt(2.0); -- precision inside 68881, but slower */
+ t1 = (x[i3] + x[i4]) * ONE_OVER_SQRT_2;
+ t2 = (x[i3] - x[i4]) * ONE_OVER_SQRT_2;
+ x[i4] = x[i2] - t1;
+ x[i3] = -x[i2] - t1;
+ x[i2] = x[i1] - t2;
+ x[i1] = x[i1] + t2;
+ label38: /* 38 CONTINUE */
+ ;
+ }
+ is = 2*id - n2;
+ id = 4*id;
+ if (is < n) goto label40; /* IF (IS.LT.N) GOTO 40 */
+ /* a = e; */
+ windex0 = 1<<(m-k);
+ windex = windex0;
+ for (j=2; j<=n8; j++) /* DO 32 J = 2,N8 */
+ {
+ /* windex = (j-1)*(1<<(m-k)); -- done with trick to avoid (j-1) and 1<<(m-k) */
+ cc1 = wcos[windex];
+ ss1 = wcos[windex_n4 - windex]; /* sin-from-cos trick: see my notes */
+ cc3 = w3cos[windex];
+ ss3 = w3sin[windex]; /* sin-from-cos trick does not work here */
+ windex = j*windex0; /* same trick as "a = j*e" */
+ /* a3 = 3*a;
+ cc1 = cos(a);
+ ss1 = sin(a);
+ cc3 = cos(a3);
+ ss3 = sin(a3);
+ a = j*e;*/
+ is = 0;
+ id = 2*n2;
+ label36: /* 36 DO 30 I = IS,N-1,ID */
+ for (i=is; i<n; i=i+id)
+ { i1 = i + j;
+ i2 = i1 + n4;
+ i3 = i2 + n4;
+ i4 = i3 + n4;
+ i5 = i + n4 - j + 2;
+ i6 = i5 + n4;
+ i7 = i6 + n4;
+ i8 = i7 + n4;
+ t1 = x[i3]*cc1 + x[i7]*ss1;
+ t2 = x[i7]*cc1 - x[i3]*ss1;
+ t3 = x[i4]*cc3 + x[i8]*ss3;
+ t4 = x[i8]*cc3 - x[i4]*ss3;
+ t5 = t1 + t3;
+ t6 = t2 + t4;
+ t3 = t3 - t1; /* t3 = t1 - t3; */
+ t4 = t2 - t4;
+ x[i8] = x[i6] + t6;
+ x[i3] = t6 - x[i6];
+ x[i4] = x[i2] + t3; /* x[i4] = x[i2] - t3; */
+ x[i7] = t3 - x[i2]; /* x[i7] = -x[i2] - t3; */
+ x[i6] = x[i1] - t5;
+ x[i1] = x[i1] + t5;
+ x[i2] = x[i5] + t4;
+ x[i5] = x[i5] - t4;
+ } /* 30 CONTINUE */
+ is = 2*id - n2;
+ id = 4*id;
+ if (is < n) goto label36; /* IF (IS.LT.N) GOTO 36 */
+ } /* 32 CONTINUE */
+ } /* 10 CONTINUE */
+} /* RETURN \r END */
+
+
+/*
+ C CONJUGATE SYMMETRIC FOURIER TRANSFORM (Split-Radix, Decimation-in-time)
+ C (adapted from Sorensen,et.al. ASSP-35 no.6 page 849, October 1986)
+ C
+ C input is [Re(0),Re(1),...,Re(n/2), Im(n/2-1),...,Im(1)]
+ C output is real
+ */
+
+/* twiddle tables identical with rft
+ for comments see rft */
+
+void pas_csft(tables_ok,flag, x,n,m, wcos,w3cos,w3sin)
+ REAL *x, *wcos,*w3cos,*w3sin;
+ long n,m, flag, tables_ok;
+{ REAL scale;
+ long i,n2;
+ void pas_realdata_make_twiddle_tables();
+ void C_Array_Time_Reverse();
+ void pas_csft_backward_loop();
+
+ if (tables_ok != 1) /* 1 means = tables already made */
+ pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin);
+
+ if (flag == 1) /* forward csft */
+ { n2 = n/2;
+ scale = (REAL) (1.0 / ((double) n));
+ for (i=0; i<=n2; i++) x[i] = x[i]*scale;
+ scale = (-scale);
+ for (i=n2+1; i<n; i++) x[i] = x[i]*scale; /* scale and conjugate cs-array */
+ pas_csft_backward_loop(x,n,m, wcos,w3cos,w3sin);
+ }
+ else /* backward csft */
+ pas_csft_backward_loop(x,n,m, wcos,w3cos,w3sin);
+}
+
+/* csft
+ forward transform === backward_loop + 1/N scaling + time-reversal
+ inverse transform === backward_loop
+ */
+
+/* wcos must be length n/4
+ w3cos, w3sin must be length n/8
+ (greater than n/8 is fine also, e.g. use cft tables)
+ */
+
+void pas_csft_backward_loop(x,n,m, wcos,w3cos,w3sin)
+ REAL *x, *wcos,*w3cos,*w3sin;
+ long n,m;
+{ /* REAL a,a3,e; no need anymore, use tables */
+ REAL r1, xt, cc1,cc3,ss1,ss3, t1,t2,t3,t4,t5;
+ long n1,n2,n4,n8, i,j,k, is,id, i0,i1,i2,i3,i4,i5,i6,i7,i8;
+ long windex0, windex, windex_n4; /* indices for twiddle tables */
+ /********** fortran indices start from 1,... **/
+ x = x-1; /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */
+ /********** fortran indices start from 1,... **/
+ /* c */
+ windex_n4 = n/4; /* need for indexing sin via wcos twiddle table */
+ /* c */
+ /* c */
+ /* c
+ c -------L-shaped-butterflies-------- */
+ n2 = 2*n;
+ for (k=1; k<m; k++) /* do 10 k = 1,m-1 */
+ { is = 0;
+ id = n2;
+ n2 = n2>>1; /* n2 = n2/2; */
+ n4 = n2>>2; /* n4 = n2/4; */
+ n8 = n4>>1; /* n8 = n4/2; */
+ /* e = 6.283185307179586476925287 / ((REAL) n2); no need anymore, use tables */
+ label17:
+ for (i=is; i<n; i=i+id) /* 17 do 15 i = is,(n-1),id */
+ { i1 = i + 1;
+ i2 = i1 + n4;
+ i3 = i2 + n4;
+ i4 = i3 + n4;
+ t1 = x[i1] - x[i3];
+ x[i1] = x[i1] + x[i3];
+ x[i2] = 2*x[i2];
+ x[i3] = t1 - 2*x[i4];
+ x[i4] = t1 + 2*x[i4];
+ if (n4 == 1) goto label15; /* if (n4.eq.1) goto 15 */
+ i1 = i1 + n8;
+ i2 = i2 + n8;
+ i3 = i3 + n8;
+ i4 = i4 + n8;
+ t1 = (x[i2] - x[i1]) * ONE_OVER_SQRT_2;
+ t2 = (-(x[i4] + x[i3])) * ONE_OVER_SQRT_2;
+ /* t1 = (x[i2] - x[i1])/sqrt(2.0);
+ t2 = (x[i4] + x[i3])/sqrt(2.0); */
+ x[i1] = x[i1] + x[i2];
+ x[i2] = x[i4] - x[i3];
+ x[i3] = 2 * (t2-t1); /* x[i3] = 2 * (-t2-t1); */
+ x[i4] = 2 * (t2+t1); /* x[i4] = 2 * (-t2+t1); */
+ label15:
+ ;
+ } /* 15 continue */
+ is = 2*id - n2;
+ id = 4*id;
+ if (is < (n-1)) goto label17; /* if (is.lt.(n-1)) goto 17 */
+ /* a = e; */
+ windex0 = 1<<(k-1); /* see my notes */
+ windex = windex0;
+ for (j=2; j<=n8; j++) /* do 20 j=2,n8 */
+ {
+ /* windex = (j-1)*(1<<(k-1)); -- done with trick to avoid (j-1) and 1<<(k-1) */
+ cc1 = wcos[windex];
+ ss1 = wcos[windex_n4 - windex]; /* sin-from-cos trick: see my notes */
+ cc3 = w3cos[windex];
+ ss3 = w3sin[windex]; /* sin-from-cos trick does not work here */
+ windex = j*windex0; /* same trick as "a = j*e" */
+ /* a3 = 3*a;
+ cc1 = cos(a);
+ ss1 = sin(a);
+ cc3 = cos(a3);
+ ss3 = sin(a3);
+ a = j*e; */
+ is = 0;
+ id = 2*n2;
+ label40:
+ for (i=is; i<n; i=i+id) /* 40 do 30 i = is,(n-1),id */
+ { i1 = i + j;
+ i2 = i1 + n4;
+ i3 = i2 + n4;
+ i4 = i3 + n4;
+ i5 = i + n4 - j + 2;
+ i6 = i5 + n4;
+ i7 = i6 + n4;
+ i8 = i7 + n4;
+ t1 = x[i1] - x[i6];
+ x[i1] = x[i1] + x[i6];
+ t2 = x[i5] - x[i2];
+ x[i5] = x[i5] + x[i2];
+ t3 = x[i8] + x[i3];
+ x[i6] = x[i8] - x[i3];
+ t4 = x[i4] + x[i7];
+ x[i2] = x[i4] - x[i7];
+ t5 = t1 - t4;
+ t1 = t1 + t4;
+ t4 = t2 - t3;
+ t2 = t2 + t3;
+ x[i3] = t5*cc1 + t4*ss1;
+ x[i7] = t5*ss1 - t4*cc1; /* x[i7] = (-t4*cc1) + t5*ss1; */
+ x[i4] = t1*cc3 - t2*ss3;
+ x[i8] = t2*cc3 + t1*ss3;
+ } /* 30 continue */
+ is = 2*id - n2;
+ id = 4*id;
+ if (is < (n-1)) goto label40; /* if (is.lt.(n-1)) goto 40 */
+ } /* 20 continue */
+ } /* 10 continue */
+ /* c */
+ /* c ----length-two-butterflies----------- */
+ is = 1;
+ id = 4;
+ label70:
+ for (i0=is; i0<=n; i0=i0+id) /* 70 DO 60 I0 = IS,N,ID */
+ { i1 = i0 + 1;
+ r1 = x[i0];
+ x[i0] = r1 + x[i1];
+ x[i1] = r1 - x[i1];
+ } /* 60 CONTINUE */
+ is = 2*id - 1;
+ id = 4*id;
+ if (is < n) goto label70; /* IF (IS.LT.N) GOTO 70 */
+ /* c */
+ /* c-----------bit-reverse-counter---------------c */
+ label100:
+ j = 1;
+ n1 = n - 1;
+ for (i=1; i<=n1; i++) /* DO 104 I = 1, N1 */
+ { if (i >= j) goto label101; /* if (i .ge. j) goto 101 */
+ xt = x[j];
+ x[j] = x[i];
+ x[i] = xt;
+ label101: k = n>>1; /* k = n/2; */
+ label102: if (k>=j) goto label103;
+ j = j - k;
+ k = k>>1; /* k = k/2; */
+ goto label102;
+ label103: j = j + k;
+ } /* 104 CONTINUE */
+ /* c */
+} /* RETURN \r END */
+
+
+\f
+
+/* Image processing only for square images
+ (old stuff handles non-square but is slow)
+ For 2d FTs precomputed tables or not, make almost no difference in total time.
+ */
+
+DEFINE_PRIMITIVE ("PAS-CFT2D!", Prim_pas_cft2d, 5,5, 0)
+{ long i, length, power, flag, rows,rowpower;
+ REAL *f1,*f2, *wcos,*w3cos,*w3sin;
+ void pas_cft2d();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* flag forward-backward */
+ CHECK_ARG (2, ARRAY_P); /* real part */
+ CHECK_ARG (3, ARRAY_P); /* imag part */
+ CHECK_ARG (4, ARRAY_P); /* twiddle tables, length = 3*(rows/4) */
+ CHECK_ARG (5, FIXNUM_P); /* (1)=tables precomputed, else recompute */
+
+ flag = Get_Integer(ARG_REF(1));
+ length = Array_Length(ARG_REF(2));
+ if (length != (Array_Length(ARG_REF(3)))) error_bad_range_arg(2);
+
+ for (power=0, i=length; i>1; power++) /* length must be power of 2 */
+ { if ( (i % 2) == 1) error_bad_range_arg(2);
+ i=i/2; }
+
+ if ((power % 2) == 1) error_bad_range_arg(2);
+ rowpower = (power/2);
+ rows = (1<<rowpower); /* square image */
+
+ f1 = Scheme_Array_To_C_Array(ARG_REF(2));
+ f2 = Scheme_Array_To_C_Array(ARG_REF(3));
+ if (f1==f2) error_wrong_type_arg(2);
+
+ wcos = Scheme_Array_To_C_Array(ARG_REF(4)); /* twiddle tables */
+ if (Array_Length(ARG_REF(4)) != (3*rows/4)) error_bad_range_arg(4);
+ w3cos = wcos + rows/4;
+ w3sin = w3cos + rows/4;
+ if ((arg_nonnegative_integer(5)) == 1)
+ pas_cft2d(1, flag, f1,f2, rows, rowpower, wcos,w3cos,w3sin);
+ else
+ pas_cft2d(0, flag, f1,f2, rows, rowpower, wcos,w3cos,w3sin);
+ /* 1 means tables are already made
+ 0 means compute new tables */
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+/* pas_cft2d
+ n = rows of square image, rows is power of 2
+ m = rowpower
+ Scaling (1/n) is done all-at-once at the end.
+ Time-Reversing is done intermediately, it is more efficient.
+ */
+void pas_cft2d(tables_ok,flag, x,y,n,m, wcos,w3cos,w3sin)
+ REAL *x,*y, *wcos,*w3cos,*w3sin;
+ long n,m, flag, tables_ok;
+{ REAL scale, *xrow,*yrow;
+ long i,j, rows,cols, total_length;
+ void pas_cft_make_twiddle_tables_once();
+ void C_Array_Time_Reverse();
+ void pas_cft_forward_loop();
+
+ if (tables_ok != 1) /* 1 means = tables already made */
+ pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin);
+
+ rows = n;
+ cols = rows; /* square image */
+ total_length = rows*rows;
+
+ if (flag != 1) /* backward transform */
+ for (i=0; i<total_length; i++) y[i] = (-y[i]); /* conjugate before */
+
+ xrow = x; yrow = y; /* ROW-WISE */
+ for (i=0; i<rows; i++) /* forward or backward */
+ { pas_cft_forward_loop( xrow,yrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols;
+ yrow = yrow + cols; }
+
+ Image_Fast_Transpose(x, rows); /* COLUMN-WISE */
+ Image_Fast_Transpose(y, rows);
+ xrow = x; yrow = y;
+ for (i=0; i<rows; i++) /* forward or backward */
+ { pas_cft_forward_loop( xrow,yrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols;
+ yrow = yrow + cols; }
+
+ Image_Fast_Transpose(x, rows);
+ Image_Fast_Transpose(y, rows);
+
+ if (flag == 1) /* forward : scale */
+ { scale = (REAL) (1.0 / ((double) total_length));
+ for (i=0; i<total_length; i++)
+ { x[i] = x[i]*scale;
+ y[i] = y[i]*scale; }}
+ else /* backward : conjugate after */
+ for (i=0; i<total_length; i++) y[i] = (-y[i]);
+}
+
+
+DEFINE_PRIMITIVE ("PAS-RFT2D-CSFT2D!", Prim_pas_rft2d_csft2d, 5,5, 0)
+{ long i, length, power, flag, ft_type, rows,rowpower;
+ REAL *f1, *wcos,*w3cos,*w3sin;
+ void pas_rft2d(), pas_csft2d();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* flag 1=forward, else backward transform */
+ CHECK_ARG (2, ARRAY_P); /* Input data (real or cs) */
+ CHECK_ARG (3, ARRAY_P); /* CFT twiddle tables, length = 3*(rows/4) */
+ CHECK_ARG (4, FIXNUM_P); /* (1)=tables precomputed, else recompute */
+ CHECK_ARG (5, FIXNUM_P); /* ft_type = 1 or 3
+ 1 means rft, 3 means csft */
+ flag = Get_Integer(ARG_REF(1));
+ f1 = Scheme_Array_To_C_Array(ARG_REF(2));
+ length = Array_Length(ARG_REF(2));
+ for (power=0, i=length; i>1; power++) /* length must be power of 2 */
+ { if ( (i % 2) == 1) error_bad_range_arg(2);
+ i=i/2; }
+
+ if ((power % 2) == 1) error_bad_range_arg(2);
+ rowpower = (power/2);
+ rows = (1<<rowpower); /* square image */
+
+ wcos = Scheme_Array_To_C_Array(ARG_REF(3)); /* CFT twiddle tables */
+ if (Array_Length(ARG_REF(3)) != (3*rows/4)) error_bad_range_arg(3);
+ w3cos = wcos + rows/4;
+ w3sin = w3cos + rows/4;
+
+ ft_type = (arg_nonnegative_integer(5)); /* rft2d or csft2d */
+ if (ft_type == 1) {
+ if ((arg_nonnegative_integer(4)) == 1)
+ pas_rft2d (1, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
+ else pas_rft2d (0, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
+ }
+ else if (ft_type == 3) {
+ if ((arg_nonnegative_integer(4)) == 1)
+ pas_csft2d (1, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
+ else pas_csft2d (0, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
+ /* 1 means tables are already made
+ 0 means compute new tables */
+ }
+ else error_bad_range_arg(5);
+
+ PRIMITIVE_RETURN (NIL);
+}
+
+/* c RFT2D CSFT2D
+ The frequencies are scrabled wrt what cft2d (and the old image-fft) give.
+ See cs-image-magnitude and cs-image-real which unscrable automatically.
+ c
+ c Implementation notes:
+ c conjugate in one domain is reverse and conjugate in other
+ c reverse in one domain is reverse in other
+ c
+ c reverse cs-array which is identical to conjugate cs-array (same domain)
+ c is reverse in other domain
+ c
+ c conjugate cs-array before csft is-better-than reverse real-array afterwards
+ c
+ c
+ c rft2d-csft2d use 1d-cft tables to compute rft
+ c cft tables are simply larger than realdata tables.
+ */
+
+/* pas_rft2d
+ n = rows of square image, rows is power of 2
+ m = rowpower
+ Scaling (1/n) is done all-at-once at the end.
+ Time-Reversing is done intermediately, it is more efficient.
+ */
+void pas_rft2d(tables_ok,flag, x, n,m, wcos,w3cos,w3sin)
+ REAL *x, *wcos,*w3cos,*w3sin;
+ long n,m, flag, tables_ok;
+{ REAL scale, *xrow,*yrow;
+ long i,j, rows,cols, total_length, n2;
+ void pas_cft_make_twiddle_tables_once();
+ void C_Array_Time_Reverse();
+ void pas_rft_forward_loop(), pas_cft_forward_loop();
+
+ if (tables_ok != 1) /* 1 means = tables already made */
+ pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin);
+
+ rows = n;
+ cols = rows; /* square image */
+ n2 = n/2;
+ total_length = rows*rows;
+
+ xrow = x; /* First ROW-WISE */
+ if (flag == 1) /* forward transform */
+ for (i=0; i<rows; i++)
+ { pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols; }
+ else /* backward transform */
+ for (i=0; i<rows; i++)
+ { pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ for (j=n2+1; j<n; j++) xrow[j] = (-xrow[j]); /* time-reverse cs-array */
+ xrow = xrow + cols; }
+
+ Image_Fast_Transpose(x, rows); /* COLUMN-WISE */
+
+ /* TREAT specially rows 0 and n2,
+ they are real and go into cs-arrays */
+ if (flag == 1) /* forward transform */
+ { xrow = x + 0 ;
+ pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ xrow = x + n2*cols;
+ pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin); }
+ else /* backward transform */
+ { xrow = x + 0 ;
+ pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ for (j=n2+1; j<n; j++) xrow[j] = (-xrow[j]); /* time-reverse cs-array */
+ xrow = x + n2*cols;
+ pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ for (j=n2+1; j<n; j++) xrow[j] = (-xrow[j]); /* time-reverse cs-array */
+ }
+
+ /* TREAT the rest of the rows with CFT
+ */
+ if (flag != 1) /* backward : conjugate before */
+ for (i=(n2+1)*cols; i<total_length; i++) x[i] = (-x[i]);
+
+ xrow = x + cols; /* real part */
+ yrow = x + (rows-1)*cols; /* imag part */
+ for (i=1; i<n2; i++) /* forward or backward transform */
+ { pas_cft_forward_loop(xrow,yrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols;
+ yrow = yrow - cols; }
+ /* DO NOT TRANSPOSE BACK, leave real-imag in horizontal rows, save.
+ */
+
+ if (flag == 1) /* forward : scale */
+ { scale = (REAL) (1.0 / ((double) total_length));
+ for (i=0; i<total_length; i++)
+ x[i] = x[i]*scale; }
+ else /* backward : conjugate after */
+ for (i=(n2+1)*cols; i<total_length; i++)
+ x[i] = (-x[i]);
+}
+
+
+/* pas_csft2d
+ n = rows of square image, rows is power of 2
+ m = rowpower
+ Scaling (1/n) is done all-at-once at the end.
+ Time-Reversing is done intermediately, it is more efficient.
+ */
+void pas_csft2d(tables_ok,flag, x, n,m, wcos,w3cos,w3sin)
+ REAL *x, *wcos,*w3cos,*w3sin;
+ long n,m, flag, tables_ok;
+{ REAL scale, *xrow,*yrow;
+ long i,j, rows,cols, total_length, n2;
+ void pas_cft_make_twiddle_tables_once();
+ void C_Array_Time_Reverse();
+ void pas_csft_backward_loop(), pas_cft_forward_loop();
+
+ if (tables_ok != 1) /* 1 means = tables already made */
+ pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin);
+
+ rows = n;
+ cols = rows; /* square image */
+ n2 = n/2;
+ total_length = rows*rows;
+
+ /* First ROW-WISE */
+
+ /* TREAT SPECIALLY ROWS 0 and n2, they are cs-arrays and they go into real */
+ if (flag == 1) /* forward transform */
+ { xrow = x + 0 ;
+ for (j=n2+1; j<n; j++) xrow[j]=(-xrow[j]); /* conjugate before */
+ pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ xrow = x + n2*cols;
+ for (j=n2+1; j<n; j++) xrow[j]=(-xrow[j]); /* conjugate before */
+ pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin); }
+ else /* backward transform */
+ { xrow = x + 0 ;
+ pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ xrow = x + n2*cols;
+ pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin); }
+
+ /* TREAT the rest of the rows with CFT
+ */
+ if (flag != 1) /* backward : conjugate before */
+ for (i=(n2+1)*cols; i<total_length; i++) x[i] = (-x[i]);
+
+ xrow = x + cols; /* real part */
+ yrow = x + (rows-1)*cols; /* imag part */
+ for (i=1; i<n2; i++) /* forward or backward transform */
+ { pas_cft_forward_loop(xrow,yrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols;
+ yrow = yrow - cols; }
+
+ if (flag != 1) /* backward : conjugate after */
+ for (i=(n2+1)*cols; i<total_length; i++) x[i] = (-x[i]);
+
+ Image_Fast_Transpose(x, rows);
+ /* Second COLUMN-WISE
+ Everything should be cs-arrays now */
+
+ xrow = x;
+ if (flag == 1) /* forward transform */
+ for (i=0; i<rows; i++)
+ { for (j=n2+1; j<n; j++) xrow[j]=(-xrow[j]); /* conjugate before */
+ pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols; }
+ else /* backward transform */
+ for (i=0; i<rows; i++)
+ { pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
+ xrow = xrow + cols; }
+
+ if (flag == 1) /* forward : scale */
+ { scale = (REAL) (1.0 / ((double) total_length));
+ for (i=0; i<total_length; i++)
+ x[i] = x[i] * scale; }
+}
+
+
+\f
+/* STUFF BEFORE 4-15-1989
+ */
+
void Make_FFT_Tables(w1, w2, n, flag)
REAL *w1, *w2; long n, flag; /* n = length of data */
{ long m, n2=n/2; /* n2 = length of w1,w2 */
The answer is left in f1, f2.
*/
+/* C_Array_FFT
+ complex data, radix=2, not-in-place.
+ (adapted from an fft program I got from Yekta)
+ */
void C_Array_FFT(flag, power, n, f1, f2, g1,g2,w1,w2)
long flag, power, n; REAL f1[], f2[], g1[], g2[], w1[], w2[];
{ long n2=n>>1, a;
}
\f
-/* CHIRP-Z-TRANSFORM (for complex data)
+/* CHIRP-Z-TRANSFORM (complex data)
*/
-#define take_modulo_one(x, answer) \
-{ long ignore_integral_part; \
- double modf(); \
- answer = (REAL) modf( ((double) x), &ignore_integral_part); }
-
-#define power_of_2_to_power(n,power) \
-{ long i; \
- for (power=0,i=n; i>1; power++) \
- { if ((i%2) == 1) { printf("\nNot a power of 2--- %d\n",n); exit(1); } \
- i=i/2; }}
-
-long smallest_power_of_2_ge(n)
- long n;
-{ long i,power;
- if (n<0) { printf("\nsmallest_pwr_of_2_ge negative argument--- %d\n", n); exit(1); }
- power=0; i=1;
- while (i<n)
- { power++; i=i*2; }
- return(power);
-}
-
-/* C_Array_CZT ------------------ Generalization of DFT.
+/* C_Array_CZT Generalization of DFT
+ ;;
Frequency is scaled as an L/2-point DFT of the input data (zero padded to L/2).
- ARGUMENTS and ASSUMPTIONS:
+ ;;
phi = starting point (on unit circle) -- Range 0,1 (covers 0,2pi like DFT angle)
rho = resolution (angular frequency spacing) -- Range 0,1 (maps 0,2pi like DFT angle)
N = input data length
f1,f2,fo1,fo2,g1,g2 must be of length L
fft_w1,fft_w2 must be of length L/2
czt_w1,czt_w2 must be of length max(M,N) ----
+ ;;
RESULT is left on f1,f2 (M complex numbers).
*/
C_Array_CZT(phi,rho, N,M,log2_L, f1,f2,fo1,fo2, g1,g2, fft_w1,fft_w2,czt_w1,czt_w2)
void CZT_Pre_Multiply(phi,rho, f1,f2, N,L) /* phi = starting point */
double phi,rho; REAL *f1,*f2; long N,L; /* this proc is two complex multiplication */
{ long i;
- double double_i, tmp, A1, A2;
+ double tmp, A1, A2;
rho = rho*.5; /* To make 1/2 in exponent "(n^2)/2" */
for (i=0;i<N;i++)
- { double_i = ((double) i);
- tmp = ((rho*double_i)+phi) * double_i;
- take_modulo_one(tmp,tmp); /* allows more decimal places */
- tmp = TWOPI * tmp;
+ { tmp = ((double) i);
+ tmp = TWOPI * ((phi + (rho*tmp)) * tmp);
A1 = cos(tmp);
A2 = sin(tmp);
tmp = A1*f1[i] - A2*f2[i];
rho = rho*.5; /* the 1/2 in the "(n^2)/2" exponent */
for (i=0;i<maxMN;i++)
{ tmp = ((double) i);
- tmp = tmp * (rho*tmp);
- take_modulo_one(tmp,tmp); /* allows more decimal places */
- tmp = TWOPI * tmp;
+ tmp = TWOPI * (tmp * (rho*tmp));
czt_w1[i] = (REAL) cos(tmp);
- czt_w2[i] = (REAL) sin(tmp);
- }}
+ czt_w2[i] = (REAL) sin(tmp); }
+}
+
-/* not used:
+long smallest_power_of_2_ge(n)
+ long n;
+{ long i,power;
+ if (n<0) { printf("\n ABORT program! smallest_pwr_of_2_ge negative argument--- %d\n", n); fflush(stdout); }
+ power=0; i=1;
+ while (i<n)
+ { power++; i=i*2; }
+ return(power);
+}
+
+/* stuff not currently used
+
void CZT_Post_Multiply(f1,f2,czt_w1,czt_w2,M)
REAL *f1,*f2,*czt_w1,*czt_w2; long M;
{ long i;
f2[i] = 2.0 * (f1[i]*czt_w2[i] + f2[i]*czt_w1[i]);
f1[i] = 2.0 * tmp;
}}
+
+#define take_modulo_one(x, answer) \
+{ long ignore_integral_part; \
+ double modf(); \
+ answer = (double) modf( ((double) x), &ignore_integral_part); }
+ ^ this only works when answer is double
+
*/
+
\f
/* 2D DFT ---------------- row-column decomposition
(3D not working yet)
Image_Transpose(Temp_Array, Real_Array, ncols, nrows); /* TRANSPOSE BACK: order of frequencies. */
}
}
-\f
+
Square_Image_2D_FFT_In_Scheme_Heap(flag, nrows, Real_Array, Imag_Array)
long flag,nrows; REAL *Real_Array, *Imag_Array;
{ REAL *f1,*f2,*g1,*g2,*w1,*w2, *Work_Here;
Image_Fast_Transpose(Real_Array, nrows); /* TRANSPOSE BACK: order of frequencies. */
Image_Fast_Transpose(Imag_Array, nrows);
}
-\f
+
C_Array_3D_FFT_In_Scheme_Heap(flag, ndeps, nrows, ncols, Real_Array, Imag_Array)
long flag, ndeps, nrows, ncols; REAL *Real_Array, *Imag_Array;
{ long l, m, n;
printf("aborted\n.");
}
}
-\f
+
Cube_Space_3D_FFT_In_Scheme_Heap(flag, ndeps, Real_Array, Imag_Array)
long flag, ndeps; REAL *Real_Array, *Imag_Array;
{ register long l, m, n;
Pointer answer;
REAL *f1,*f2,*g1,*g2,*w1,*w2;
REAL *Work_Here;
- Primitive_3_Args();
- Arg_1_Type(TC_FIXNUM); /* flag */
- Arg_2_Type(TC_ARRAY); /* real */
- Arg_3_Type(TC_ARRAY); /* imag */
- Set_Time_Zone(Zone_Math);
- flag = Get_Integer(Arg1);
- length = Array_Length(Arg2);
- if (length != (Array_Length(Arg3))) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- power=0;
+ PRIMITIVE_HEADER (4);
+ flag = arg_nonnegative_integer(1); /* forward or backward */
+ CHECK_ARG (2, ARRAY_P); /* input real */
+ CHECK_ARG (3, ARRAY_P); /* input imag */
+
+ length = Array_Length(ARG_REF(2));
+ if (length != (Array_Length(ARG_REF(3)))) error_bad_range_arg(2);
+
for (power=0, i=length; i>1; power++) {
- if ( (i % 2) == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
+ if ( (i % 2) == 1) error_bad_range_arg(2);
i=i/2; }
- f1 = Scheme_Array_To_C_Array(Arg2);
- f2 = Scheme_Array_To_C_Array(Arg3);
- if (f1==f2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
+ f1 = Scheme_Array_To_C_Array(ARG_REF(2));
+ f2 = Scheme_Array_To_C_Array(ARG_REF(3));
+ if (f1==f2) error_wrong_type_arg(2);
Primitive_GC_If_Needed(length*3*REAL_SIZE);
Work_Here = (REAL *) Free;
C_Array_FFT(flag, power, length, f1,f2,g1,g2,w1,w2);
- Primitive_GC_If_Needed(4);
- answer = Make_Pointer(TC_LIST, Free);
- *Free++ = Arg2;
- *Free = Make_Pointer(TC_LIST, Free+1);
- Free += 1;
- *Free++ = Arg3;
- *Free++ = NIL;
- return answer;
+ PRIMITIVE_RETURN (NIL);
}
-DEFINE_PRIMITIVE ("ARRAY-CZT", Prim_array_czt, 5, 5, 0)
+DEFINE_PRIMITIVE ("ARRAY-CZT!", Prim_array_czt, 6,6, 0)
{ double phi,rho;
- long N,M,L;
- long log2_L,maxMN,smallest_power_of_2_ge(), allocated_cells;
- int Error_Number, Scheme_Number_To_Double();
- Pointer answer, answer_1,answer_2;
- REAL *f1,*f2,*fo1,*fo2, *g1,*g2, *fft_w1,*fft_w2,*czt_w1,*czt_w2;
- REAL *Work_Here;
- Primitive_5_Args();
- Arg_3_Type(TC_FIXNUM);
- Arg_4_Type(TC_ARRAY); /* real */
- Arg_5_Type(TC_ARRAY); /* imag */
- Set_Time_Zone(Zone_Math);
-
- Error_Number = Scheme_Number_To_Double(Arg1, &phi); /* phi=starting point [0,1]*/
- if (Error_Number == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- Error_Number = Scheme_Number_To_Double(Arg2, &rho); /* rho=resolution [0,1] */
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
-
- M = Get_Integer(Arg3); /* M=length of output data */
- N = Array_Length(Arg4); /* N=length of input data */
- if (N != (Array_Length(Arg5))) Primitive_Error(ERR_ARG_5_BAD_RANGE);
- if ((M+N-1)<4) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ long N,M,L, i;
+ long log2_L, maxMN;
+ long smallest_power_of_2_ge();
+ int errcode;
+ REAL *a,*b,*c,*d;
+ REAL *f1,*f2,*fo1,*fo2, *g1,*g2, *fft_w1,*fft_w2,*czt_w1,*czt_w2, *Work_Here;
+
+ PRIMITIVE_HEADER (6);
+ /* 1 phi=starting point [0,1]*/
+ /* 2 rho=resolution [0,1] */
+ CHECK_ARG (3, ARRAY_P); /* input real */
+ CHECK_ARG (4, ARRAY_P); /* input imag */
+ CHECK_ARG (5, ARRAY_P); /* output real */
+ CHECK_ARG (6, ARRAY_P); /* output imag */
+
+ errcode = Scheme_Number_To_Double(ARG_REF(1), &phi);
+ if (errcode==1) error_bad_range_arg(1); if (errcode==2) error_wrong_type_arg(1);
+ errcode = Scheme_Number_To_Double(ARG_REF(2), &rho);
+ if (errcode==1) error_bad_range_arg(2); if (errcode==2) error_wrong_type_arg(2);
+
+ a = Scheme_Array_To_C_Array(ARG_REF(3));
+ b = Scheme_Array_To_C_Array(ARG_REF(4));
+ c = Scheme_Array_To_C_Array(ARG_REF(5));
+ d = Scheme_Array_To_C_Array(ARG_REF(6));
+
+ N = Array_Length(ARG_REF(3)); /* N = input length */
+ M = Array_Length(ARG_REF(5)); /* M = output length */
+ if (N!=(Array_Length(ARG_REF(4)))) error_bad_range_arg(3);
+ if (M!=(Array_Length(ARG_REF(6)))) error_bad_range_arg(5);
+
+ if ((M+N-1) < 4) error_bad_range_arg(5);
log2_L = smallest_power_of_2_ge(M+N-1);
L = 1<<log2_L; /* length of intermediate computation arrays */
- maxMN = max(M,N); /* length of czt tables */
+ maxMN = (((M)<(N)) ? (N) : (M)); /* length of czt tables = maximum(M,N) */
Primitive_GC_If_Needed( ((7*L) + (2*maxMN)) * REAL_SIZE);
- Work_Here = (REAL *) Free;
- g1 = Work_Here; /* We will allocate answer arrays here, after czt is done. */
- g2 = Work_Here + L;
- fo1 = Work_Here + (2*L);
- fo2 = Work_Here + (3*L);
- fft_w1 = Work_Here + (4*L);
- fft_w2 = Work_Here + (4*L) + (L/2); /* length of fft tables = L/2*/
- czt_w1 = Work_Here + (5*L);
- czt_w2 = Work_Here + (5*L) + maxMN;
- f1 = Work_Here + (5*L) + (2*maxMN); /* CZT stores its results here */
- f2 = Work_Here + (6*L) + (2*maxMN);
-
- C_Array_Copy( (Scheme_Array_To_C_Array(Arg4)), f1, N);
- C_Array_Copy( (Scheme_Array_To_C_Array(Arg5)), f2, N);
+ g1 = (REAL *) Free;
+ g2 = g1 + L;
+ fo1 = g2 + L;
+ fo2 = fo1 + L;
+ fft_w1 = fo2 + L;
+ fft_w2 = fft_w1 + (L/2);
+ czt_w1 = fft_w2 + (L/2);
+ czt_w2 = czt_w1 + maxMN;
+ f1 = czt_w2 + maxMN; /* CZT stores its results here */
+ f2 = f1 + L;
+
+ for (i=0; i<N; i++) { f1[i] = a[i]; /* input data */
+ f2[i] = b[i]; }
+
C_Array_CZT(phi,rho, N,M,log2_L, f1,f2,fo1,fo2, g1,g2, fft_w1,fft_w2,czt_w1,czt_w2);
- Allocate_Array(answer_1, M, allocated_cells); /* Overwrite g1,g2 which started at Free */
- Allocate_Array(answer_2, M, allocated_cells);
- C_Array_Copy(f1, (Scheme_Array_To_C_Array(answer_1)), M);
- C_Array_Copy(f2, (Scheme_Array_To_C_Array(answer_2)), M);
+ for (i=0; i<M; i++) { c[i] = f1[i]; /* results */
+ d[i] = f2[i]; }
- Primitive_GC_If_Needed(4);
- answer = Make_Pointer(TC_LIST, Free);
- *Free++ = answer_1;
- *Free = Make_Pointer(TC_LIST, Free+1);
- Free += 1;
- *Free++ = answer_2;
- *Free++ = NIL;
- return answer;
+ PRIMITIVE_RETURN (NIL);
}
DEFINE_PRIMITIVE ("ARRAY-2D-FFT!", Prim_array_2d_fft, 5, 5, 0)
promotional, or sales literature without prior written consent from
MIT in each case. */
-/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/image.c,v 9.27 1988/08/15 20:49:26 cph Exp $ */
+/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/image.c,v 9.28 1989/06/22 21:52:26 pas Rel $ */
#include "scheme.h"
#include "prims.h"
}
}
-DEFINE_PRIMITIVE ("SUBIMAGE", Prim_subimage, 5, 5, 0)
-{ long Length, new_Length;
- long i,j;
- Pointer Pnrows, Pncols, Prest, Parray;
- long lrow, hrow, lcol, hcol;
- long nrows, ncols, new_nrows, new_ncols;
- REAL *Array, *To_Here;
- Pointer Result, Array_Data_Result, *Orig_Free;
- int Error_Number;
- long allocated_cells;
-
- Primitive_5_Args();
- Arg_1_Type(TC_LIST); /* image = (nrows ncols array) */
- Pnrows = Vector_Ref(Arg1, CONS_CAR);
- Prest = Vector_Ref(Arg1, CONS_CDR);
- Pncols = Vector_Ref(Prest, CONS_CAR);
- Prest = Vector_Ref(Prest, CONS_CDR);
- Parray = Vector_Ref(Prest, CONS_CAR);
- if (Vector_Ref(Prest, CONS_CDR) != NIL) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- if (Type_Code(Parray) != TC_ARRAY) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
-
- Range_Check(nrows, Pnrows, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Range_Check(ncols, Pncols, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Range_Check(lrow, Arg2, 0, nrows, ERR_ARG_2_BAD_RANGE);
- Range_Check(hrow, Arg3, lrow, nrows, ERR_ARG_3_BAD_RANGE);
- Range_Check(lcol, Arg4, 0, ncols, ERR_ARG_4_BAD_RANGE);
- Range_Check(hcol, Arg5, lcol, ncols, ERR_ARG_5_BAD_RANGE);
- new_nrows = hrow - lrow +1;
- new_ncols = hcol - lcol +1;
- new_Length = new_nrows * new_ncols;
-
- /* ALLOCATE SPACE */
- Primitive_GC_If_Needed(6);
- Orig_Free = Free;
- Free += 6;
- Result = Make_Pointer(TC_LIST, Orig_Free);
- *Orig_Free++ = Make_Non_Pointer(TC_FIXNUM, new_nrows);
- *Orig_Free = Make_Pointer(TC_LIST, Orig_Free+1);
- Orig_Free++;
- *Orig_Free++ = Make_Non_Pointer(TC_FIXNUM, new_ncols);
- *Orig_Free = Make_Pointer(TC_LIST, Orig_Free+1);
- Orig_Free++;
- Allocate_Array(Array_Data_Result, new_Length, allocated_cells);
- *Orig_Free++ = Array_Data_Result;
- *Orig_Free = NIL;
- /* END ALLOCATION */
-
- Array = Scheme_Array_To_C_Array(Parray);
- To_Here = Scheme_Array_To_C_Array(Array_Data_Result);
- for (i=lrow; i<=hrow; i++) {
- for (j=lcol; j<=hcol; j++) {
- *To_Here++ = Array[i*ncols+j]; /* A(i,j)--->Array[i*ncols+j] */
- }}
-
- return Result;
-}
-
-/* The following does not work properly, to be fixed if need.
+/* The following does not work, to be fixed.
*/
DEFINE_PRIMITIVE ("IMAGE-DOUBLE-TO-FLOAT!", Prim_image_double_to_float, 1, 1, 0)
{ long Length;
*To_Here = ((float) temp_value_cell);
To_Here++;
}
-
/* and now SIDE-EFFECT the ARRAY_HEADER */
allocated_cells = (Length *
((sizeof(Pointer)+sizeof(float)-1) / sizeof(Pointer)) +
return Arg1;
}
-DEFINE_PRIMITIVE ("IMAGE-SET-ROW!", Prim_image_set_row, 3, 3, 0)
-{ long Length, i,j;
- Pointer Pnrows, Pncols, Prest, Parray;
- long nrows, ncols, row_to_set;
- REAL *Array, *Row_Array;
-
- Primitive_3_Args();
- Arg_1_Type(TC_LIST); /* image = (nrows ncols array) */
- Pnrows = Vector_Ref(Arg1, CONS_CAR);
- Prest = Vector_Ref(Arg1, CONS_CDR);
- Pncols = Vector_Ref(Prest, CONS_CAR);
- Prest = Vector_Ref(Prest, CONS_CDR);
- Parray = Vector_Ref(Prest, CONS_CAR);
- if (Vector_Ref(Prest, CONS_CDR) != NIL) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- if (Type_Code(Parray) != TC_ARRAY) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
-
- Range_Check(nrows, Pnrows, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Range_Check(ncols, Pncols, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Arg_2_Type(TC_FIXNUM);
- Range_Check(row_to_set, Arg2, 0, (nrows-1), ERR_ARG_2_BAD_RANGE);
- Arg_3_Type(TC_ARRAY);
- Row_Array = Scheme_Array_To_C_Array(Arg3);
- if (Array_Length(Arg3)>ncols) Primitive_Error(ERR_ARG_3_BAD_RANGE);
-
- Array = Scheme_Array_To_C_Array(Parray);
- C_Image_Set_Row(Array, row_to_set, Row_Array, nrows, ncols);
- return Arg1;
+
+\f
+DEFINE_PRIMITIVE ("SUBIMAGE-COPY!",
+ Prim_subimage_copy, 12,12, 0)
+{ long r1,c1, r2,c2, at1r,at1c,at2r,at2c, mr,mc;
+ long nn;
+ REAL *x,*y;
+ void subimage_copy();
+ PRIMITIVE_HEADER (12);
+ CHECK_ARG (1, FIXNUM_P); /* rows 1 */
+ CHECK_ARG (2, FIXNUM_P); /* cols 1 */
+ CHECK_ARG (3, ARRAY_P); /* image array 1 = source array */
+ CHECK_ARG (4, FIXNUM_P); /* rows 2 */
+ CHECK_ARG (5, FIXNUM_P); /* cols 2 */
+ CHECK_ARG (6, ARRAY_P); /* image array 2 = destination array */
+
+ CHECK_ARG (7, FIXNUM_P); /* at1 row */
+ CHECK_ARG (8, FIXNUM_P); /* at1 col */
+ CHECK_ARG (9, FIXNUM_P); /* at2 row */
+ CHECK_ARG (10, FIXNUM_P); /* at2 col */
+ CHECK_ARG (11, FIXNUM_P); /* m row */
+ CHECK_ARG (12, FIXNUM_P); /* m col */
+
+ x = Scheme_Array_To_C_Array(ARG_REF(3));
+ y = Scheme_Array_To_C_Array(ARG_REF(6));
+ r1 = arg_nonnegative_integer(1);
+ c1 = arg_nonnegative_integer(2);
+ r2 = arg_nonnegative_integer(4);
+ c2 = arg_nonnegative_integer(5);
+
+ nn = r1*c1;
+ if (nn != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+ nn = r2*c2;
+ if (nn != Array_Length(ARG_REF(6))) error_bad_range_arg(6);
+
+ at1r = arg_nonnegative_integer(7);
+ at1c = arg_nonnegative_integer(8);
+ at2r = arg_nonnegative_integer(9);
+ at2c = arg_nonnegative_integer(10);
+ mr = arg_nonnegative_integer(11);
+ mc = arg_nonnegative_integer(12);
+
+ if (((at1r+mr)>r1) || ((at1c+mc)>c1)) error_bad_range_arg(7);
+ if (((at2r+mr)>r2) || ((at2c+mc)>c2)) error_bad_range_arg(9);
+
+ subimage_copy(x,y, r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc);
+
+ PRIMITIVE_RETURN (NIL);
}
-DEFINE_PRIMITIVE ("IMAGE-SET-COLUMN!", Prim_image_set_column, 3, 3, 0)
-{ long Length, i,j;
- Pointer Pnrows, Pncols, Prest, Parray;
- long nrows, ncols, col_to_set;
- REAL *Array, *Col_Array;
-
- Primitive_3_Args();
- Arg_1_Type(TC_LIST); /* image = (nrows ncols array) */
- Pnrows = Vector_Ref(Arg1, CONS_CAR);
- Prest = Vector_Ref(Arg1, CONS_CDR);
- Pncols = Vector_Ref(Prest, CONS_CAR);
- Prest = Vector_Ref(Prest, CONS_CDR);
- Parray = Vector_Ref(Prest, CONS_CAR);
- if (Vector_Ref(Prest, CONS_CDR) != NIL) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- if (Type_Code(Parray) != TC_ARRAY) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
-
- Range_Check(nrows, Pnrows, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Range_Check(ncols, Pncols, 0, 1024, ERR_ARG_1_BAD_RANGE);
+void subimage_copy(x,y, r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc)
+ REAL *x,*y; long r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc;
+{ long i,j;
+ REAL *xrow,*yrow;
- Arg_2_Type(TC_FIXNUM);
- Range_Check(col_to_set, Arg2, 0, (nrows-1), ERR_ARG_2_BAD_RANGE);
- Arg_3_Type(TC_ARRAY);
- Col_Array = Scheme_Array_To_C_Array(Arg3);
- if (Array_Length(Arg3)>ncols) Primitive_Error(ERR_ARG_3_BAD_RANGE);
+ xrow = x + at1r*c1 + at1c;
+ yrow = y + at2r*c2 + at2c; /* A(i,j)--->Array[i*ncols+j] */
- Array = Scheme_Array_To_C_Array(Parray);
- C_Image_Set_Col(Array, col_to_set, Col_Array, nrows, ncols);
- return Arg1;
+ for (i=0; i<mr; i++) {
+ for (j=0; j<mc; j++) yrow[j] = xrow[j];
+ xrow = xrow + c1;
+ yrow = yrow + c2;
+ }
}
-C_Image_Set_Row(Image_Array, row_to_set, Row_Array, nrows, ncols) REAL *Image_Array, *Row_Array;
-long nrows, ncols, row_to_set;
-{ long j;
- REAL *From_Here, *To_Here;
-
- To_Here = &Image_Array[row_to_set*ncols];
- From_Here = Row_Array;
- for (j=0;j<ncols;j++)
- *To_Here++ = *From_Here++;
-}
-C_Image_Set_Col(Image_Array, col_to_set, Col_Array, nrows, ncols) REAL *Image_Array, *Col_Array;
-long nrows, ncols, col_to_set;
-{ long i;
- REAL *From_Here, *To_Here;
- To_Here = &Image_Array[col_to_set];
- From_Here = Col_Array;
- for (i=0;i<nrows;i++) {
- *To_Here = *From_Here++;
- To_Here += nrows;
- }
-}
+/* image-operation-2
+ groups together procedures that use 2 image-arrays
+ (usually side-effecting the 2nd image, but not necessarily)
+ */
-DEFINE_PRIMITIVE ("IMAGE-LAPLACIAN", Prim_image_laplacian, 1, 1, 0)
-{ long nrows, ncols, Length;
- Pointer Pnrows, Pncols, Prest, Parray;
- REAL *Array, *To_Here;
- Pointer Result, Array_Data_Result, *Orig_Free;
- long allocated_cells;
- /* */
- Primitive_1_Args();
- Arg_1_Type(TC_LIST); /* image = (nrows ncols array) */
- Pnrows = Vector_Ref(Arg1, CONS_CAR);
- Prest = Vector_Ref(Arg1, CONS_CDR);
- Pncols = Vector_Ref(Prest, CONS_CAR);
- Prest = Vector_Ref(Prest, CONS_CDR);
- Parray = Vector_Ref(Prest, CONS_CAR);
- if (Vector_Ref(Prest, CONS_CDR) != NIL) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- if (Type_Code(Parray) != TC_ARRAY) Primitive_Error(ERR_ARG_1_WRONG_TYPE);
- /* */
- Range_Check(nrows, Pnrows, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Range_Check(ncols, Pncols, 0, 1024, ERR_ARG_1_BAD_RANGE);
- Length=nrows*ncols;
- /* */
- /* ALLOCATE SPACE */
- Primitive_GC_If_Needed(6);
- Orig_Free = Free;
- Free += 6;
- Result = Make_Pointer(TC_LIST, Orig_Free);
- *Orig_Free++ = Make_Non_Pointer(TC_FIXNUM, nrows);
- *Orig_Free = Make_Pointer(TC_LIST, Orig_Free+1);
- Orig_Free++;
- *Orig_Free++ = Make_Non_Pointer(TC_FIXNUM, ncols);
- *Orig_Free = Make_Pointer(TC_LIST, Orig_Free+1);
- Orig_Free++;
- Allocate_Array(Array_Data_Result, Length, allocated_cells);
- *Orig_Free++ = Array_Data_Result;
- *Orig_Free = NIL;
- /* END ALLOCATION */
- /* */
- Array = Scheme_Array_To_C_Array(Parray);
- C_image_laplacian(Array, (Scheme_Array_To_C_Array(Array_Data_Result)), nrows, ncols);
- PRIMITIVE_RETURN(Result);
+DEFINE_PRIMITIVE ("IMAGE-OPERATION-2!",
+ Prim_image_operation_2, 5,5, 0)
+{ long rows, cols, nn, opcode;
+ REAL *x,*y;
+ void image_laplacian();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* operation opcode */
+ CHECK_ARG (2, FIXNUM_P); /* rows */
+ CHECK_ARG (3, FIXNUM_P); /* cols */
+ CHECK_ARG (4, ARRAY_P); /* image array 1 */
+ CHECK_ARG (5, ARRAY_P); /* image array 2 */
+
+ x = Scheme_Array_To_C_Array(ARG_REF(4));
+ y = Scheme_Array_To_C_Array(ARG_REF(5));
+ rows = arg_nonnegative_integer(2);
+ cols = arg_nonnegative_integer(3);
+ nn = rows*cols;
+ if (nn != Array_Length(ARG_REF(4))) error_bad_range_arg(4);
+ if (nn != Array_Length(ARG_REF(5))) error_bad_range_arg(5);
+
+ opcode = arg_nonnegative_integer(1);
+
+ if (opcode==1)
+ image_laplacian(x,y,rows,cols); /* result in y */
+ else if (opcode==2)
+ error_bad_range_arg(1); /* illegal opcode */
+ else
+ error_bad_range_arg(1); /* illegal opcode */
+
+ PRIMITIVE_RETURN (NIL);
}
/* Laplacian form [4,-1,-1,-1,-1]/4
A(i,j) --> Array[i*ncols + j]
With no knowledge outside boundary, assume laplace(edge-point)=0.0 (no wrap-around, no artificial bndry)
*/
-C_image_laplacian(array, new_array, nrows, ncols)
- REAL *array, *new_array;
+void image_laplacian(x,y, nrows,ncols)
+ REAL *x, *y;
long nrows, ncols;
{ long i,j, nrows1, ncols1;
nrows1=nrows-1; ncols1=ncols-1;
- if ((nrows<2)||(ncols<2)) return(1); /* no need todo anything for 1-point image */
+ if ((nrows<2)||(ncols<2)) return; /* no need todo anything for 1-point image */
/* */
- i=0;j=0; new_array[i*ncols+j] = 0.0; /* NE corner */
- i=0;j=ncols1; new_array[i*ncols+j] = 0.0; /* NW corner */
- i=nrows1;j=0; new_array[i*ncols+j] = 0.0; /* SE corner */
- i=nrows1;j=ncols1; new_array[i*ncols+j] = 0.0; /* SW corner */
- i=0; for (j=1;j<ncols1;j++) new_array[i*ncols+j] = 0.0; /* NORTH row */
- i=nrows1; for (j=1;j<ncols1;j++) new_array[i*ncols+j] = 0.0; /* SOUTH row */
- j=0; for (i=1;i<nrows1;i++) new_array[i*ncols+j] = 0.0; /* EAST column */
- j=ncols1; for (i=1;i<nrows1;i++) new_array[i*ncols+j] = 0.0; /* WEST column */
+ i=0;j=0; y[i*ncols+j] = 0.0; /* NE corner */
+ i=0;j=ncols1; y[i*ncols+j] = 0.0; /* NW corner */
+ i=nrows1;j=0; y[i*ncols+j] = 0.0; /* SE corner */
+ i=nrows1;j=ncols1; y[i*ncols+j] = 0.0; /* SW corner */
+ i=0; for (j=1;j<ncols1;j++) y[i*ncols+j] = 0.0; /* NORTH row */
+ i=nrows1; for (j=1;j<ncols1;j++) y[i*ncols+j] = 0.0; /* SOUTH row */
+ j=0; for (i=1;i<nrows1;i++) y[i*ncols+j] = 0.0; /* EAST column */
+ j=ncols1; for (i=1;i<nrows1;i++) y[i*ncols+j] = 0.0; /* WEST column */
/* */
for (i=1;i<nrows1;i++)
- for (j=1;j<ncols1;j++) new_array[i*ncols+j] = /* interior of image */
- array[i*ncols+j] - (.25)*(array[i*ncols+(j-1)] + array[i*ncols+(j+1)] + array[(i-1)*ncols+j] + array[(i+1)*ncols+j]);
+ for (j=1;j<ncols1;j++) y[i*ncols+j] = /* interior of image */
+ x[i*ncols+j] - (.25)*(x[i*ncols+(j-1)] + x[i*ncols+(j+1)] + x[(i-1)*ncols+j] + x[(i+1)*ncols+j]);
}
+
DEFINE_PRIMITIVE ("IMAGE-DOUBLE-BY-INTERPOLATION", Prim_image_double_by_interpolation, 1, 1, 0)
{ long nrows, ncols, Length;
Pointer Pnrows, Pncols, Prest, Parray;
Vector_Set(Prest, CONS_CAR, Make_Pointer(TC_FIXNUM, nrows) );
return Arg1;
}
-\f
+
DEFINE_PRIMITIVE ("IMAGE-ROTATE-90CLW!", Prim_image_rotate_90clw, 1, 1, 0)
{ long Length;
Pointer Pnrows, Pncols, Prest, Parray;
Vector_Set(Prest, CONS_CAR, Make_Pointer(TC_FIXNUM, nrows) );
return Arg1;
}
-\f
+
DEFINE_PRIMITIVE ("IMAGE-ROTATE-90CCLW!", Prim_image_rotate_90cclw, 1, 1, 0)
{ long Length;
Pointer Pnrows, Pncols, Prest, Parray;
Vector_Set(Prest, CONS_CAR, Make_Pointer(TC_FIXNUM, nrows) );
return Arg1;
}
-\f
+
DEFINE_PRIMITIVE ("IMAGE-MIRROR!", Prim_image_mirror, 1, 1, 0)
{ long Length;
Pointer Pnrows, Pncols, Prest, Parray;
return Arg1;
}
-\f
-/* THE C ROUTINES THAT DO THE REAL WORK */
+
+/* C routines referred to above */
/*
IMAGE_FAST_TRANSPOSE
Array[to] = temp;
}}
}
-\f
+
/*
IMAGE_TRANSPOSE
A(i,j) -> B(j,i) .
New_Array[j*nrows + i] = Array[i*ncols + j]; /* (columns transposed-image) = nrows */
}}
}
-\f
+
/*
IMAGE_ROTATE_90CLW
A(i,j) <-> A(j, (nrows-1)-i) .
Rotated_Array[(j*nrows) + ((nrows-1)-i)] = Array[i*ncols+j]; /* (columns rotated_image) =nrows */
}}
}
-\f
+
/*
ROTATION 90degrees COUNTER-CLOCK-WISE:
A(i,j) <-> A((nrows-1)-j, i) . (minus 1 because we start from 0).
Rotated_Array[to_index] = Array[from_index];
}}
}
-\f
+
/*
IMAGE_MIRROR:
A(i,j) <-> A(i, (ncols-1)-j) [ The -1 is there because we count from 0] .
}}
}
-
-\f
/*
IMAGE_ROTATE_90CLW_MIRROR:
A(i,j) <-> A(j, i) this should be identical to image_transpose (see above).
Rotated_Array[to] = Array[from];
}}
}
-\f
-
-
-/* END */
+/* More Image Manipulation -----------------------
+ */
+DEFINE_PRIMITIVE ("SQUARE-IMAGE-TIME-REVERSE!",
+ Prim_square_image_time_reverse, 2,2, 0)
+{ long i, rows;
+ REAL *a;
+ void square_image_time_reverse();
+ PRIMITIVE_HEADER (2);
+ CHECK_ARG (1, ARRAY_P);
+ CHECK_ARG (2, FIXNUM_P);
+ a = Scheme_Array_To_C_Array(ARG_REF(1));
+ rows = arg_nonnegative_integer(2);
+ if ((rows*rows) != Array_Length(ARG_REF(1))) error_bad_range_arg(1);
+ square_image_time_reverse(a,rows);
+
+ PRIMITIVE_RETURN (NIL);
+}
+/* Square Image Time reverse
+ is combination of one-dimensional time-reverse
+ row-wise and column-wise.
+ It can be done slightly more efficiently than below.
+ */
+void square_image_time_reverse(x,rows)
+ REAL *x;
+ long rows;
+{ long i,cols;
+ REAL *xrow, *yrow;
+ void C_Array_Time_Reverse();
+ cols = rows; /* square image */
+
+ xrow = x;
+ for (i=0; i<rows; i++) /* row-wise */
+ { C_Array_Time_Reverse(xrow,cols);
+ xrow = xrow + cols; }
+
+ Image_Fast_Transpose(x, rows);
+
+ xrow = x;
+ for (i=0; i<rows; i++) /* column-wise */
+ { C_Array_Time_Reverse(xrow,cols);
+ xrow = xrow + cols; }
+
+ Image_Fast_Transpose(x, rows);
+}
+\f
+/* cs-images
+ */
+/* operation-1
+ groups together procedures that operate on 1 cs-image-array
+ (side-effecting the image)
+ */
-/*
-\f
-DEFINE_PRIMITIVE ("SAMPLE-PERIODIC-2D-FUNCTION", Prim_sample_periodic_2d_function, 4, 4, 0)
-{ long N, i, allocated_cells, Function_Number;
- REAL Signal_Frequency, Sampling_Frequency, DT, DTi;
- REAL twopi = 6.28318530717958, twopi_f_dt;
- Pointer Result, Pfunction_number, Psignal_frequency;
- Pointer Pfunction_Number;
- int Error_Number;
- REAL *To_Here, unit_square_wave(), unit_triangle_wave();
-
- Primitive_4_Args();
- Arg_1_Type(TC_FIXNUM);
- Arg_4_Type(TC_FIXNUM);
- Range_Check(Function_Number, Arg1, 0, 10, ERR_ARG_1_BAD_RANGE); / * fix this * /
-
- Error_Number = Scheme_Number_To_REAL(Arg2, &Signal_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
- if (Signal_Frequency == 0) Primitive_Error(ERR_ARG_2_BAD_RANGE);
-
- Error_Number = Scheme_Number_To_REAL(Arg3, &Sampling_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
- if (Sampling_Frequency == 0) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- DT = (1 / Sampling_Frequency);
- twopi_f_dt = twopi * Signal_Frequency * DT;
-
- Range_Check(N, Arg4, 0, ARRAY_MAX_LENGTH, ERR_ARG_4_BAD_RANGE);
-
- allocated_cells = (N*REAL_SIZE) + ARRAY_HEADER_SIZE;
- Primitive_GC_If_Needed(allocated_cells);
-
- Result = Make_Pointer(TC_ARRAY, Free);
- Free[ARRAY_HEADER] = Make_Non_Pointer(TC_MANIFEST_ARRAY, allocated_cells-1);
- Free[ARRAY_LENGTH] = N;
- To_Here = Scheme_Array_To_C_Array(Result);
- Free = Free+allocated_cells;
-
- DT = twopi_f_dt;
- if (Function_Number == 0)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = cos(DTi);
- else if (Function_Number == 1)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = sin(DTi);
- else if (Function_Number == 2)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = unit_square_wave(DTi);
- else if (Function_Number == 3)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = unit_triangle_wave(DTi);
+DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-1!",
+ Prim_cs_image_operation_1, 3,3, 0)
+{ long rows, opcode;
+ REAL *a;
+ void cs_image_magnitude(), cs_image_real_part(), cs_image_imag_part();
+ PRIMITIVE_HEADER (3);
+ CHECK_ARG (1, FIXNUM_P); /* operation opcode */
+ CHECK_ARG (2, FIXNUM_P); /* rows */
+ CHECK_ARG (3, ARRAY_P); /* input and output image array */
+
+ a = Scheme_Array_To_C_Array(ARG_REF(3));
+ rows = arg_nonnegative_integer(2); /* square images only */
+ if ((rows*rows) != Array_Length(ARG_REF(3))) error_bad_range_arg(1);
+ opcode = arg_nonnegative_integer(1);
+
+ if (opcode==1)
+ cs_image_magnitude(a,rows);
+ else if (opcode==2)
+ cs_image_real_part(a,rows);
+ else if (opcode==3)
+ cs_image_imag_part(a,rows);
else
- Primitive_Error(ERR_ARG_1_BAD_RANGE);
+ error_bad_range_arg(3); /* illegal opcode */
- return Result;
+ PRIMITIVE_RETURN (NIL);
}
-*/
-/* END IMAGE PROCESSING */
-
-\f
-/* Note for the macro: To1 and To2 must BE Length1-1, and Length2-2 RESPECTIVELY ! */
-/*
-#define C_Convolution_Point_Macro(X, Y, To1, To2, N, Result) \
-{ long Min_of_N_To1=min((N),(To1)); \
- long mi, N_minus_mi; \
- REAL Sum=0.0; \
- for (mi=max(0,(N)-(To2)), N_minus_mi=(N)-mi; mi <= Min_of_N_To1; mi++, N_minus_mi--) \
- Sum += (X[mi] * Y[N_minus_mi]); \
- (Result)=Sum; \
-}
-\f
-DEFINE_PRIMITIVE ("CONVOLUTION-POINT", Prim_convolution_point, 3, 3, 0)
-{ long Length1, Length2, N;
- REAL *Array1, *Array2;
- REAL C_Result;
-
- Primitive_3_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Arg_3_Type(TC_FIXNUM);
- Length1 = Array_Length(Arg1);
- Length2 = Array_Length(Arg2);
- N = Get_Integer(Arg3);
- Array1 = Scheme_Array_To_C_Array(Arg1);
- Array2 = Scheme_Array_To_C_Array(Arg2);
- C_Convolution_Point_Macro(Array1, Array2, Length1-1, Length2-1, N, C_Result);
- Reduced_Flonum_Result(C_Result);
+void cs_array_real_part(a,n)
+ REAL *a; long n;
+{ long i,n2; /* works both for even and odd length */
+ n2 = n/2;
+ for (i=n2+1;i<n;i++) a[i] = a[n-i]; /* copy real values into place */
+ /* even signal */
}
-\f
-DEFINE_PRIMITIVE ("ARRAY-CONVOLUTION", Prim_array_convolution, 2, 2, 0)
-{ long Endpoint1, Endpoint2, allocated_cells, i;
- / * ASSUME A SIGNAL FROM INDEX 0 TO ENDPOINT=LENGTH-1 * /
- long Resulting_Length;
- REAL *Array1, *Array2, *To_Here;
- Pointer Result;
-
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_ARRAY);
- Endpoint1 = Array_Length(Arg1) - 1;
- Endpoint2 = Array_Length(Arg2) - 1;
- Resulting_Length = Endpoint1 + Endpoint2 + 1;
- Array1 = Scheme_Array_To_C_Array(Arg1);
- Array2 = Scheme_Array_To_C_Array(Arg2);
-
- allocated_cells = (Resulting_Length * REAL_SIZE) + ARRAY_HEADER_SIZE;
- Primitive_GC_If_Needed(allocated_cells);
- Result = Make_Pointer(TC_ARRAY, Free);
- Free[ARRAY_HEADER] = Make_Non_Pointer(TC_MANIFEST_ARRAY, allocated_cells-1);
- Free[ARRAY_LENGTH] = Resulting_Length;
- Free += allocated_cells;
- To_Here = Scheme_Array_To_C_Array(Result);
-
- for (i=0; i<Resulting_Length; i++) {
- C_Convolution_Point_Macro(Array1, Array2, Endpoint1, Endpoint2, i, *To_Here);
- To_Here++;
- }
- return Result;
+
+void cs_array_imag_part(a,n)
+ REAL *a; long n;
+{ long i,n2;
+ n2 = n/2; /* integer division truncates down */
+ for (i=n2+1; i<n; i++) /* works both for even and odd length */
+ { a[n-i] = a[i]; /* copy imaginary values into place */
+ a[i] = (-a[i]); } /* odd signal */
+ a[0] = 0.0;
+ if (2*n2 == n) /* even length, n2 is real only */
+ a[n2] = 0.0;
}
-*/
-/* m_pi = 3.14159265358979323846264338327950288419716939937510; */
-/*
-DEFINE_PRIMITIVE ("SAMPLE-PERIODIC-FUNCTION", Prim_sample_periodic_function, 4, 4, 0)
-{ long N, i, allocated_cells, Function_Number;
- REAL Signal_Frequency, Sampling_Frequency, DT, DTi;
- REAL twopi = 6.28318530717958, twopi_f_dt;
- Pointer Result, Pfunction_number, Psignal_frequency;
- Pointer Pfunction_Number;
- int Error_Number;
- REAL *To_Here, unit_square_wave(), unit_triangle_wave();
-
- Primitive_4_Args();
- Arg_1_Type(TC_FIXNUM);
- Arg_4_Type(TC_FIXNUM);
- Range_Check(Function_Number, Arg1, 0, 10, ERR_ARG_1_BAD_RANGE); / * fix this * /
-
- Error_Number = Scheme_Number_To_REAL(Arg2, &Signal_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
- if (Signal_Frequency == 0) Primitive_Error(ERR_ARG_2_BAD_RANGE);
-
- Error_Number = Scheme_Number_To_REAL(Arg3, &Sampling_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE);
- if (Sampling_Frequency == 0) Primitive_Error(ERR_ARG_3_BAD_RANGE);
- DT = (1 / Sampling_Frequency);
- twopi_f_dt = twopi * Signal_Frequency * DT;
-
- Range_Check(N, Arg4, 0, ARRAY_MAX_LENGTH, ERR_ARG_4_BAD_RANGE);
-
- allocated_cells = (N*REAL_SIZE) + ARRAY_HEADER_SIZE;
- Primitive_GC_If_Needed(allocated_cells);
-
- Result = Make_Pointer(TC_ARRAY, Free);
- Free[ARRAY_HEADER] = Make_Non_Pointer(TC_MANIFEST_ARRAY, allocated_cells-1);
- Free[ARRAY_LENGTH] = N;
- To_Here = Scheme_Array_To_C_Array(Result);
- Free = Free+allocated_cells;
-
- DT = twopi_f_dt;
- if (Function_Number == 0)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = cos(DTi);
- else if (Function_Number == 1)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = sin(DTi);
- else if (Function_Number == 2)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = unit_square_wave(DTi);
- else if (Function_Number == 3)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = unit_triangle_wave(DTi);
- else
- Primitive_Error(ERR_ARG_1_BAD_RANGE);
-
- return Result;
-}
-\f
-REAL hamming(t, length) REAL t, length;
-{ REAL twopi = 6.28318530717958;
- REAL pi = twopi/2.;
- REAL t_bar = cos(twopi * (t / length));
- if ((t<length) && (t>0.0)) return(.08 + .46 * (1 - t_bar));
- else return (0);
-}
-\f
-REAL hanning(t, length) REAL t, length;
-{ REAL twopi = 6.28318530717958;
- REAL pi = twopi/2.;
- REAL t_bar = cos(twopi * (t / length));
- if ((t<length) && (t>0.0))
- return(.5 * (1 - t_bar));
- else return (0);
-}
-\f
-REAL unit_square_wave(t) REAL t;
-{ REAL twopi = 6.28318530717958;
- REAL fmod(), fabs();
- REAL pi = twopi/2.;
- REAL t_bar = fabs(fmod(t, twopi));
- if (t_bar < pi) return(1);
- else return(0);
+
+/* From now on (below), assume that cs-images (rows=cols) have always EVEN LENGTH
+ which is true when they come from FFTs
+ */
+
+
+
+/* In the following 3 time-reverse the bottom half rows
+ is done to match the frequencies of complex-images
+ coming from cft2d.
+ Also transpose is needed to match frequencies identically
+
+ #|
+ ;; Scrabling of frequencies in cs-images
+
+ ;; start from real image 4x4
+
+ ;; rft2d is a cs-image
+ (3.5 .375 -2.75 1.875 -.25 0. 0. -.25 -.25 -.125 0. .125 .25 .25 0. 0.)
+
+ ;; cft2d transposed
+ ;; real
+ 3.5 .375 -2.75 .375
+ -.25 0. 0. -.25 ; same as cs-image
+ -.25 -.125 0. -.125
+ -.25 -.25 0. 0. ; row3 = copy 1 + time-reverse
+ ;; imag
+ 0. 1.875 0. -1.875
+ .25 .25 0. 0. ; same as cs-image
+ 0. .125 0. -.125
+ -.25 0. 0. -.25 ; row 3 = copy 1 + negate + time-reverse
+ |#
+
+ */
+
+void cs_image_magnitude(x,rows)
+ REAL *x;
+ long rows;
+{ long i,j, cols, n,n2, nj; /* result = real ordinary image */
+ REAL *xrow, *yrow;
+ cols = rows; /* input cs-image is square */
+ n = rows;
+ n2 = n/2;
+
+ xrow = x;
+ cs_array_magnitude(xrow, n); /* row 0 is cs-array */
+ xrow = x + n2*cols;
+ cs_array_magnitude(xrow, n); /* row n2 is cs-array */
+
+ xrow = x + cols; /* real part */
+ yrow = x + (rows-1)*cols; /* imag part */
+ for (i=1; i<n2; i++) {
+ xrow[ 0] = (REAL) sqrt((double) xrow[ 0]*xrow[ 0] + yrow[ 0]*yrow[ 0]);
+ xrow[n2] = (REAL) sqrt((double) xrow[n2]*xrow[n2] + yrow[n2]*yrow[n2]);
+ yrow[ 0] = xrow[ 0];
+ yrow[n2] = xrow[n2];
+ for (j=1; j<n2; j++) {
+ nj = n-j;
+ xrow[ j] = (REAL) sqrt((double) xrow[ j]*xrow[ j] + yrow[ j]*yrow[ j]);
+ xrow[nj] = (REAL) sqrt((double) xrow[nj]*xrow[nj] + yrow[nj]*yrow[nj]);
+ yrow[j] = xrow[nj];
+ yrow[nj] = xrow[ j]; /* Bottom rows: copy (even) and time-reverse */
+ }
+ xrow = xrow + cols;
+ yrow = yrow - cols; }
+ Image_Fast_Transpose(x, n);
}
-\f
-REAL unit_triangle_wave(t) REAL t;
-{ REAL twopi = 6.28318530717958;
- REAL pi = twopi/2.;
- REAL t_bar = fabs(fmod(t, twopi));
- if (t_bar < pi) return( t_bar / pi );
- else return( (twopi - t_bar) / pi );
+
+
+void cs_image_real_part(x,rows)
+ REAL *x;
+ long rows;
+{ long i,j,cols, n,n2;
+ REAL *xrow, *yrow;
+ void cs_array_real_part();
+ cols = rows; /* square image */
+ n = rows;
+ n2 = n/2;
+
+ xrow = x;
+ cs_array_real_part(xrow, n); /* row 0 is cs-array */
+ xrow = x + n2*cols;
+ cs_array_real_part(xrow, n); /* row n2 is cs-array */
+
+ xrow = x + cols; /* real part */
+ yrow = x + (rows-1)*cols; /* imag part */
+ for (i=1; i<n2; i++) {
+ yrow[0] = xrow[0]; /* copy real part into imaginary's place (even) */
+ for (j=1; j<n; j++)
+ yrow[j] = xrow[n-j]; /* Bottom rows: copy and time-reverse */
+ xrow = xrow + cols;
+ yrow = yrow - cols; }
+ Image_Fast_Transpose(x, n);
}
-\f
-DEFINE_PRIMITIVE ("SAMPLE-APERIODIC-FUNCTION", Prim_sample_aperiodic_function, 3, 3, 0)
-{ long N, i, allocated_cells, Function_Number;
- REAL Sampling_Frequency, DT, DTi;
- REAL twopi = 6.28318530717958;
- Pointer Result;
- int Error_Number;
- REAL *To_Here, twopi_dt;
- Primitive_3_Args();
- Arg_1_Type(TC_FIXNUM);
- Arg_3_Type(TC_FIXNUM);
- Range_Check(Function_Number, Arg1, 0, 6, ERR_ARG_1_BAD_RANGE);
-
- Error_Number = Scheme_Number_To_REAL(Arg2, &Sampling_Frequency);
- if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE);
- if (Sampling_Frequency == 0) Primitive_Error(ERR_ARG_2_BAD_RANGE);
- DT = (1 / Sampling_Frequency);
- twopi_dt = twopi * DT;
+void cs_image_imag_part(x,rows)
+ REAL *x;
+ long rows;
+{ long i,j,cols, n,n2, nj;
+ REAL *xrow, *yrow;
+ void cs_array_imag_part();
+ cols = rows; /* square image */
+ n = rows;
+ n2 = n/2;
+
+ xrow = x;
+ cs_array_imag_part(xrow, n); /* row 0 is cs-array */
+ xrow = x + n2*cols;
+ cs_array_imag_part(xrow, n); /* row n2 is cs-array */
+
+ xrow = x + cols; /* real part */
+ yrow = x + (rows-1)*cols; /* imag part */
+ for (i=1; i<n2; i++) {
+ xrow[0] = yrow[0]; /* copy the imaginary part into real's place */
+ xrow[n2] = yrow[n2];
+ yrow[0] = (-yrow[0]); /* negate (odd) */
+ yrow[n2] = (-yrow[n2]);
+ for (j=1;j<n2; j++) {
+ nj = n-j;
+ xrow[j] = yrow[j]; /* copy the imaginary part into real's place */
+ xrow[nj] = yrow[nj];
+ yrow[j] = (-xrow[nj]); /* Bottom rows: negate (odd) and time-reverse */
+ yrow[nj] = (-xrow[j]); }
+ xrow = xrow + cols;
+ yrow = yrow - cols; }
+ Image_Fast_Transpose(x, n);
+}
- Range_Check(N, Arg3, 0, ARRAY_MAX_LENGTH, ERR_ARG_3_BAD_RANGE);
+\f
+/* cs-image-operation-2
+ groups together procedures that use 2 cs-image-arrays
+ (usually side-effecting the 2nd image, but not necessarily)
+ */
- allocated_cells = (N*REAL_SIZE) + ARRAY_HEADER_SIZE;
- Primitive_GC_If_Needed(allocated_cells);
-
- Result = Make_Pointer(TC_ARRAY, Free);
- Free[ARRAY_HEADER] = Make_Non_Pointer(TC_MANIFEST_ARRAY, allocated_cells-1);
- Free[ARRAY_LENGTH] = N;
- To_Here = Scheme_Array_To_C_Array(Result);
- Free = Free+allocated_cells;
-
- DT = twopi_dt;
- if (Function_Number == 0)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = rand();
- else if (Function_Number == 1)
- { REAL length=DT*N;
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = hanning(DTi, length);
- }
- else if (Function_Number == 2)
- { REAL length=DT*N;
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = hamming(DTi, length);
- }
- else if (Function_Number == 3)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = sqrt(DTi);
- else if (Function_Number == 4)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = log(DTi);
- else if (Function_Number == 5)
- for (i=0, DTi=0.0; i < N; i++, DTi += DT)
- *To_Here++ = exp(DTi);
+DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-2!",
+ Prim_cs_image_operation_2, 4,4, 0)
+{ long rows, nn, opcode;
+ REAL *x,*y;
+ void cs_image_multiply_into_second_one();
+ PRIMITIVE_HEADER (4);
+ CHECK_ARG (1, FIXNUM_P); /* operation opcode */
+ CHECK_ARG (2, FIXNUM_P); /* rows */
+ CHECK_ARG (3, ARRAY_P); /* image array 1 */
+ CHECK_ARG (4, ARRAY_P); /* image array 2 */
+
+ x = Scheme_Array_To_C_Array(ARG_REF(3));
+ y = Scheme_Array_To_C_Array(ARG_REF(4));
+ rows = arg_nonnegative_integer(2); /* square images only */
+ nn = rows*rows;
+ if (nn != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+ if (nn != Array_Length(ARG_REF(4))) error_bad_range_arg(4);
+
+ opcode = arg_nonnegative_integer(1);
+
+ if (opcode==1)
+ cs_image_multiply_into_second_one(x,y,rows); /* result in y */
+ else if (opcode==2)
+ error_bad_range_arg(1); /* illegal opcode */
else
- Primitive_Error(ERR_ARG_1_BAD_RANGE);
+ error_bad_range_arg(1); /* illegal opcode */
- return Result;
+ PRIMITIVE_RETURN (NIL);
}
-\f
-DEFINE_PRIMITIVE ("ARRAY-PERIODIC-DOWNSAMPLE", Prim_array_periodic_downsample, 2, 2, 0)
-{ long Length, Pseudo_Length, Sampling_Ratio;
- REAL *Array, *To_Here;
- Pointer Result;
- long allocated_cells, i, array_index;
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_FIXNUM);
- Length = Array_Length(Arg1);
- Sign_Extend(Arg2, Sampling_Ratio); / * Sampling_Ratio = integer ratio of sampling_frequencies * /
- Sampling_Ratio = Sampling_Ratio % Length; / * periodicity * /
- if (Sampling_Ratio < 1) Primitive_Error(ERR_ARG_2_BAD_RANGE);
-
- Array = Scheme_Array_To_C_Array(Arg1);
- Allocate_Array(Result, Length, allocated_cells);
- To_Here = Scheme_Array_To_C_Array(Result);
-
- Pseudo_Length = Length * Sampling_Ratio;
- for (i=0; i<Pseudo_Length; i += Sampling_Ratio) { / * new Array has the same Length by assuming periodicity * /
- array_index = i % Length;
- *To_Here++ = Array[array_index];
+void cs_image_multiply_into_second_one(x,y, rows)
+ REAL *x,*y;
+ long rows;
+{ long i,j,cols, n,n2;
+ REAL *xrow,*yrow, *xrow_r, *xrow_i, *yrow_r, *yrow_i, temp;
+ cols = rows; /* square image */
+ n = rows;
+ n2 = n/2;
+
+ xrow= x; yrow= y;
+ cs_array_multiply_into_second_one(xrow,yrow, n,n2); /* row 0 */
+
+ xrow= x+n2*cols; yrow= y+n2*cols;
+ cs_array_multiply_into_second_one(xrow,yrow, n,n2); /* row n2 */
+
+ xrow_r= x+cols; yrow_r= y+cols;
+ xrow_i= x+(n-1)*cols; yrow_i= y+(n-1)*cols;
+ for (i=1; i<n2; i++) {
+ for (j=0; j<n; j++) {
+ temp = xrow_r[j]*yrow_r[j] - xrow_i[j]*yrow_i[j]; /* real part */
+ yrow_i[j] = xrow_r[j]*yrow_i[j] + xrow_i[j]*yrow_r[j]; /* imag part */
+ yrow_r[j] = temp; }
+ xrow_r= xrow_r+cols; yrow_r= yrow_r+cols;
+ xrow_i= xrow_i-cols; yrow_i= yrow_i-cols;
}
-
- return Result;
}
-\f
-DEFINE_PRIMITIVE ("ARRAY-PERIODIC-SHIFT", Prim_array_periodic_shift, 2, 2, 0)
-{ long Length, Shift;
- REAL *Array, *To_Here;
- Pointer Result;
- long allocated_cells, i, array_index;
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_FIXNUM);
- Length = Array_Length(Arg1);
- Sign_Extend(Arg2, Shift);
- Shift = Shift % Length; / * periodic waveform, same sign as dividend * /
- Array = Scheme_Array_To_C_Array(Arg1);
- Allocate_Array(Result, Length, allocated_cells);
- To_Here = Scheme_Array_To_C_Array(Result);
-
- for (i=0; i<Length; i++) { / * new Array has the same Length by assuming periodicity * /
- array_index = (i+Shift) % Length;
- if (array_index<0) array_index = Length + array_index; / * wrap around * /
- *To_Here++ = Array[array_index];
- }
+/*
+ cs-image-operation-2x! is just like cs-image-operation-2!
+ but takes an additional flonum argument.
+ */
+
+DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-2x!",
+ Prim_cs_image_operation_2x, 5,5, 0)
+{ long rows, nn, opcode;
+ REAL *x,*y, flonum_arg;
+ int errcode;
+ void cs_image_divide_into_z();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* operation opcode */
+ CHECK_ARG (2, FIXNUM_P); /* rows */
+ CHECK_ARG (3, ARRAY_P); /* image array 1 */
+ CHECK_ARG (4, ARRAY_P); /* image array 2 */
+
+ errcode = Scheme_Number_To_REAL(ARG_REF(5), &flonum_arg); /* extra argument */
+ if (errcode==1) error_bad_range_arg(5); if (errcode==2) error_wrong_type_arg(5);
+
+ x = Scheme_Array_To_C_Array(ARG_REF(3));
+ y = Scheme_Array_To_C_Array(ARG_REF(4));
+ rows = arg_nonnegative_integer(2); /* square images only */
+ nn = rows*rows;
+ if (nn != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+ if (nn != Array_Length(ARG_REF(4))) error_bad_range_arg(4);
+
+ opcode = arg_nonnegative_integer(1);
+
+ if (opcode==1)
+ cs_image_divide_into_z( x,y, x, rows, flonum_arg); /* result in x */
+ else if (opcode==2)
+ cs_image_divide_into_z( x,y, y, rows, flonum_arg); /* result in y */
+ else
+ error_bad_range_arg(1); /* illegal opcode */
- return Result;
+ PRIMITIVE_RETURN (NIL);
}
-\f
-/ * this should really be done in SCHEME using ARRAY-MAP ! * /
-DEFINE_PRIMITIVE ("ARRAY-APERIODIC-DOWNSAMPLE", Prim_array_aperiodic_downsample, 2, 2, 0)
-{ long Length, New_Length, Sampling_Ratio;
- REAL *Array, *To_Here;
- Pointer Result;
- long allocated_cells, i, array_index;
- Primitive_2_Args();
- Arg_1_Type(TC_ARRAY);
- Arg_2_Type(TC_FIXNUM);
- Length = Array_Length(Arg1);
- Range_Check(Sampling_Ratio, Arg2, 1, Length, ERR_ARG_2_BAD_RANGE);
-
- Array = Scheme_Array_To_C_Array(Arg1);
- New_Length = Length / Sampling_Ratio;
- / * greater than zero * /
- Allocate_Array(Result, New_Length, allocated_cells);
- To_Here = Scheme_Array_To_C_Array(Result);
-
- for (i=0; i<Length; i += Sampling_Ratio) {
- *To_Here++ = Array[i];
+/* The convention for inf values in division 1/0
+ is just like in arrays
+ */
+
+void cs_image_divide_into_z(x,y, z, rows, inf) /* z can be either x or y */
+ REAL *x,*y,*z, inf;
+ long rows;
+{ long i,j,cols, n,n2;
+ REAL temp, radius;
+ REAL *ar_,*ai_, *br_,*bi_, *zr_,*zi_; /* Letters a,b correspond to x,y */
+ REAL *xrow,*yrow,*zrow;
+ cols = rows; /* square image */
+ n = rows;
+ n2 = n/2;
+
+ xrow= x; yrow= y; zrow= z;
+ cs_array_divide_into_z( xrow,yrow, zrow, n,n2, inf); /* row 0 */
+
+ xrow= x+n2*cols; yrow= y+n2*cols; zrow= z+n2*cols;
+ cs_array_divide_into_z( xrow,yrow, zrow, n,n2, inf); /* row n2 */
+
+ ar_= x+cols; br_= y+cols; zr_= z+cols;
+ ai_= x+(n-1)*cols; bi_= y+(n-1)*cols; zi_= z+(n-1)*cols;
+ for (i=1; i<n2; i++) {
+ for (j=0; j<n; j++) {
+ radius = br_[j]*br_[j] + bi_[j]*bi_[j]; /* b^2 denominator = real^2 + imag^2 */
+
+ if (radius == 0.0) {
+ if (ar_[j] == 0.0) zr_[j] = 1.0;
+ else zr_[j] = ar_[j] * inf;
+ if (ai_[j] == 0.0) zi_[j] = 1.0;
+ else zi_[j] = ai_[j] * inf; }
+ else {
+ temp = ar_[j]*br_[j] + ai_[j]*bi_[j];
+ zi_[j] = (ai_[j]*br_[j] - ar_[j]*bi_[j]) / radius; /* imag part */
+ zr_[j] = temp / radius; /* real part */
+ }}
+ ar_= ar_+cols; br_= br_+cols; zr_= zr_+cols;
+ ai_= ai_-cols; bi_= bi_-cols; zi_= zi_-cols;
}
-
- return Result;
}
-\f
-/ * ARRAY-APERIODIC-SHIFT can be done in scheme using subarray, and array-append * /
+\f
+/* operation-3
+ groups together procedures that use 3 cs-image-arrays
+ (usually side-effecting the 3rd image, but not necessarily)
+ */
-for UPSAMPLING
-if ((Length % Sampling_Ratio) != 0) Primitive_Error(ERR_ARG_2_BAD_RANGE);
-UNIMPLEMENTED YET
+DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-3!",
+ Prim_cs_image_operation_3, 5,5, 0)
+{ long rows, nn, opcode;
+ REAL *x,*y,*z;
+ void tr_complex_image_to_cs_image();
+ PRIMITIVE_HEADER (5);
+ CHECK_ARG (1, FIXNUM_P); /* operation opcode */
+ CHECK_ARG (2, FIXNUM_P); /* rows */
+ CHECK_ARG (3, ARRAY_P); /* image array 1 */
+ CHECK_ARG (4, ARRAY_P); /* image array 2 */
+ CHECK_ARG (5, ARRAY_P); /* image array 3 */
+
+ x = Scheme_Array_To_C_Array(ARG_REF(3));
+ y = Scheme_Array_To_C_Array(ARG_REF(4));
+ z = Scheme_Array_To_C_Array(ARG_REF(5));
+ rows = arg_nonnegative_integer(2); /* square images only */
+ nn = rows*rows;
+ if (nn != Array_Length(ARG_REF(3))) error_bad_range_arg(3);
+ if (nn != Array_Length(ARG_REF(4))) error_bad_range_arg(4);
+ if (nn != Array_Length(ARG_REF(5))) error_bad_range_arg(5);
+
+ opcode = arg_nonnegative_integer(1);
+
+ if (opcode==1)
+ tr_complex_image_to_cs_image(x,y, z,rows); /* result in z */
+ else if (opcode==2)
+ error_bad_range_arg(1); /* illegal opcode */
+ else
+ error_bad_range_arg(1); /* illegal opcode */
+
+ PRIMITIVE_RETURN (NIL);
+}
-*/
+/* x and y must be ALREADY TRANSPOSED real and imaginary parts
+ */
+void tr_complex_image_to_cs_image(x,y, z,rows)
+ REAL *x,*y,*z;
+ long rows;
+{ long i,j,cols, n,n2, n2_1_n;
+ REAL *xrow, *yrow, *zrow;
+ cols = rows; /* square image */
+ n = rows;
+ n2 = n/2;
+
+ xrow= x; yrow= y; zrow= z;
+ for (j=0; j<=n2; j++) zrow[j] = xrow[j]; /* real part of row 0 (cs-array) */
+ for (j=n2+1; j<n; j++) zrow[j] = yrow[n-j]; /* imag part of row 0 */
+
+ xrow= x+n2*cols; yrow= y+n2*cols; zrow= z+n2*cols;
+ for (j=0; j<=n2; j++) zrow[j] = xrow[j]; /* real part of row n2 (cs-array) */
+ for (j=n2+1; j<n; j++) zrow[j] = yrow[n-j]; /* imag part of row n2 */
+
+ xrow= x+cols; zrow= z+cols; n2_1_n = (n2-1)*cols;
+ for (j=0; j<n2_1_n; j++) zrow[j] = xrow[j]; /* real rows 1,2,..,n2-1 */
+
+ yrow= y+(n2-1)*cols; zrow= z+(n2+1)*cols; /* imag rows n2+1,n2+2,... */
+ for (i=1; i<n2; i++) {
+ for (j=0; j<n; j++) zrow[j] = yrow[j];
+ zrow = zrow + cols;
+ yrow = yrow - cols;
+ }
+}
-/* END OF FILE */
projects / mit-scheme.git / commitdiff