/* ** by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993 ** Rewritten into Xanadu by H. Xi, 23 May 1999 */ /********************************************************/ /* A Duhamel-Hollman split-radix dif fft */ /* Ref: Electronics Letters, Jan. 5, 1984 */ /* Complex input and output data in arrays x and y */ /* Length is n. */ /********************************************************/ global pi: float = 3.1415926536 {size:int | size >= 2} int dfft(px[size+1]: float, py[size+1]: float, size: int(size)) { var: [i:int | i > 0] int(i) i,j,is,id,m,n; [i:nat] int(i) k,n2,n4; int i0,i1,i2,i3,tmp; float a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt;; i = 2; m = 1; while (i < size) { i = i+i; m = m+1; } n = i; assert (n == size); n2 = 2 * n; for (k = 1; k <= m-1; k = k+1) { /* n2 = n2 / 2; n4 = n2 / 4; */ n4 = n2 / 8; n2 = 4 * n4; e = 2.0 *. pi /. float_of_int(n2); a = 0.0; for (j = 1; j <= n4 ; j = j+1) { a3 = 3.0 *. a; cc1 = cos(a); ss1 = sin(a); cc3 = cos(a3); ss3 = sin(a3); a = e *. float_of_int(j); is = j; id = 2 * n2; while (is < n) { i0 = is; i1 = i0 + n4; i2 = i1 + n4; i3 = i2 + n4; invariant: [i0:int,i1:int,i2:int,i3:int | 0 < i0 <= i1 <= i2 <= i3] (i0: int(i0), i1: int(i1), i2: int(i2), i3: int(i3)) while (i3 <= n) { r1 = px[i0] -. px[i2]; px[i0] = px[i0] +. px[i2]; r2 = px[i1] -. px[i3]; px[i1] = px[i1] +. px[i3]; s1 = py[i0] -. py[i2]; py[i0] = py[i0] +. py[i2]; s2 = py[i1] -. py[i3]; py[i1] = py[i1] +. py[i3]; s3 = r1 -. s2; r1 = r1 +. s2; s2 = r2 -. s1; r2 = r2 +. s1; px[i2] = r1 *. cc1 -. s2 *. ss1; py[i2] = -.s2 *. cc1 -. r1 *. ss1; px[i3] = s3 *.cc3 +. r2 *. ss3; py[i3] = r2 *. cc3 -. s3 *. ss3; i0 = i0 + n4; i1 = i1 + n4; i2 = i2 + n4; i3 = i3 + n4; } tmp = 2 * id - n2 + j; assert (tmp > 0); is = tmp; id = 4 * id; } } } /************************************/ /* Last stage, length=2 butterfly */ /************************************/ is = 1; id = 4; while (is < n) { invariant: [i0:int | 0 < i0] (i0: int(i0)) for (i0 = is; i0 < n; i0 = i0 + id) { i1 = i0 + 1; r1 = px[i0]; px[i0] = r1 +. px[i1]; px[i1] = r1 -. px[i1]; r1 = py[i0]; py[i0] = r1 +. py[i1]; py[i1] = r1 -. py[i1]; } is = 2 * id - 1; id = 4 * id; } /*************************/ /* Bit reverse counter */ /*************************/ j = 1; invariant: [j:nat | j <= size] (j: int(j)) for (i = 1; i < n; i = i+1) { if (i < j) { xt = px[j]; px[j] = px[i]; px[i] = xt; xt = py[j]; py[j] = py[i]; py[i] = xt; } k = n / 2; invariant: [k:nat | 2 * k <= size] /* {metric: max(j-k)} */ (k: int(k)) while (k < j) { /* assert (k > 0); */ j = j - k; k = k / 2; } j = j + k; } /* for (i = 0; i < 16; i = i + 1) printf("%d %g %gn",i, px[i], py[i]); */ return(n); }