2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/lino_czt.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 *********************************************************
12 * CHIRP-Z TRANSFORM ALGORITHM *
14 * THIS ROUTINE IMPLEMENTS THE CHIRP-Z TRANSFORM *
15 * ALGORITHM. USE OF THIS ALGORITHM PERMITS *
16 * COMPUTATION OF THE DFT OF A GIVEN *
17 * INPUT SEQUENCE ALONG A MORE GENERAL CONTOUR IN *
20 * REFERENCE: RABINER, SCHAFER, AND RADER, *
21 * IEEE TRANSACTIONS ON *
22 * AUDIO AND ELECTROACOUSTICS, VOL. *
25 * VERSION 0.0 : DAVID A. ROBERTS, 7/87 *
27 *********************************************************
29 INPUT - X: ARRAY CONTAINING (COMPLEX) SAMPLES
30 DOLD: REAL SAMPLE-SPACING
31 DNEW: DESIRED FOURIER SPACING
32 2*NOLD+1: NUMBER OF INPUT SAMPLES
33 2*NNEW+1: NUMBER OF FOURIER SAMPLES
35 OUTPUT - ODATA: ARRAY CONTAINING FOURIER SAMPLES
37 NOTE: ALL COMPLEX QUANTITIES ARE STORED IN (real,imag) FORM.
40 // parameter (pi = 3.141592654)
50 void lino_class::czt(REAL* x, REAL* odata, REAL dold, REAL dnew, INTEGER nold, INTEGER nnew)
52 // bug 195 - correct array sizes - swr - 1/17/06
67 pdd = Consts.pi * dold * dnew;
70 while (!(L >= (3 * nold + 1)))
79 // FORM Y(N) SEQUENCE - L VALUES
84 // bug 195 - correct array sizes - swr - 1/17/06
89 for (k = 0; k <= nold; k++)
92 rl[k] = (REAL) cos(phase);
93 im[k] = (REAL) sin(phase);
94 y[i] = rl[k] * x[ii] - im[k] * x[ii + 1];
95 y[i + 1] = im[k] * x[ii] + rl[k] * x[ii + 1];
100 for (k = 1; k <= 2 * (L - 2 * nold - 1); k++)
107 for (k = -nold; k <= -1; k++)
111 y[i] = rl1 * x[ii] - im1 * x[ii + 1];
112 y[i + 1] = im1 * x[ii] + rl1 * x[ii + 1];
117 // COMPUTE L-POINT FFT OF Y(N)
121 // FORM V(N) SEQUENCE - L VALUES
123 REAL v[2 * L]; // bug 195 - correct array sizes - swr - 1/17/06
126 for (k = 0; k <= L / 2; k++)
129 v[i] = (REAL) cos(phase);
130 v[i + 1] = (REAL) sin(phase);
135 for (k = -(L / 2) + 1; k <= -1; k++)
138 v[i + 1] = v[ii + 1];
143 // COMPUTE L-POINT FFT OF V(N)
147 // MULTIPLY THE ARRAYS Y AND V, POINT-BY-POINT; CALL THE RESULT G(N)
149 REAL g[2 * L]; // bug 195 - correct array sizes - swr - 1/17/06
151 for (k = 0; k < L; k++)
153 g[i] = v[i] * y[i] - v[i + 1] * y[i + 1];
154 g[i + 1] = v[i + 1] * y[i] + v[i] * y[i + 1];
158 // COMPUTE L-POINT INVERSE FFT OF G(N)
162 // MULTIPLY G(N) BY (dold)*exp(-i*pi*dold*dnew*k*k)
163 // MAKE SURE THAT ODATA IS IN PROPER FORMAT (DC-VALUE IN MIDDLE)
167 for (k = 0; k <= nnew; k++)
171 odata[i] = g[ii] * rl1 - g[ii + 1] * im1;
172 odata[i + 1] = g[ii + 1] * rl1 + g[ii] * im1;
178 ii = 2 * L - 2 * nnew;
179 for (k = -nnew; k <= -1; k++)
183 odata[i] = g[ii] * rl1 - g[ii + 1] * im1;
184 odata[i + 1] = g[ii + 1] * rl1 + g[ii] * im1;