src/snark/lino_czt.cpp

   1 /*
   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) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  *********************************************************
  11  *                                                       *
  12  *           CHIRP-Z TRANSFORM ALGORITHM                 *
  13  *                                                       *
  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    *
  18  *     THE Z-PLANE.                                      *
  19  *                                                       *
  20  *     REFERENCE:   RABINER, SCHAFER, AND RADER,         *
  21  *                IEEE TRANSACTIONS ON                   *
  22  *                AUDIO AND ELECTROACOUSTICS, VOL.       *
  23  *                AU-17, NO.2                            *
  24  *                                                       *
  25  *     VERSION 0.0 : DAVID A. ROBERTS, 7/87              *
  26  *                                                       *
  27  *********************************************************
  28
  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
  34
  35  OUTPUT - ODATA: ARRAY CONTAINING FOURIER SAMPLES
  36
  37  NOTE: ALL COMPLEX QUANTITIES ARE STORED IN (real,imag) FORM.
  38  */
  39
  40 //      parameter (pi = 3.141592654)
  41 #include <cstdio>
  42 #include <cmath>
  43
  44 #include "blkdta.h"
  45 #include "snfft.h"
  46 #include "consts.h"
  47
  48 #include "lino.h"
  49
  50 void lino_class::czt(REAL* x, REAL* odata, REAL dold, REAL dnew, INTEGER nold, INTEGER nnew)
  51 {
  52         // bug 195 - correct array sizes - swr - 1/17/06
  53         REAL rl1, im1;
  54         REAL pdd;
  55
  56         REAL phase;
  57
  58         INTEGER n3[3];
  59         INTEGER L;
  60         INTEGER offset;
  61
  62         INTEGER i;
  63         INTEGER ii;
  64         INTEGER k;
  65
  66         offset = 2 * nold;
  67         pdd = Consts.pi * dold * dnew;
  68         L = 4;
  69
  70         while (!(L >= (3 * nold + 1)))
  71         {
  72                 L = L * 2;
  73         }
  74
  75         n3[0] = L;
  76         n3[1] = 1;
  77         n3[2] = 1;
  78
  79         // FORM Y(N) SEQUENCE - L VALUES
  80
  81         i = 0;
  82         ii = offset;
  83
  84         // bug 195 - correct array sizes - swr - 1/17/06
  85         REAL rl[nold + 2];
  86         REAL im[nold + 2];
  87         REAL y[2 * L];
  88
  89         for (k = 0; k <= nold; k++)
  90         {
  91                 phase = -pdd * k * 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];
  96                 i += 2;
  97                 ii += 2;
  98         }
  99
 100         for (k = 1; k <= 2 * (L - 2 * nold - 1); k++)
 101         {
 102                 y[i] = 0.0;
 103                 i++;
 104         }
 105
 106         ii = 0;
 107         for (k = -nold; k <= -1; k++)
 108         {
 109                 rl1 = rl[-k];
 110                 im1 = im[-k];
 111                 y[i] = rl1 * x[ii] - im1 * x[ii + 1];
 112                 y[i + 1] = im1 * x[ii] + rl1 * x[ii + 1];
 113                 i += 2;
 114                 ii += 2;
 115         }
 116
 117         // COMPUTE L-POINT FFT OF Y(N)
 118
 119         snfft(y, n3, 1);
 120
 121         // FORM V(N) SEQUENCE - L VALUES
 122
 123         REAL v[2 * L];   // bug 195 - correct array sizes - swr - 1/17/06
 124
 125         i = 0;
 126         for (k = 0; k <= L / 2; k++)
 127         {
 128                 phase = pdd * k * k;
 129                 v[i] = (REAL) cos(phase);
 130                 v[i + 1] = (REAL) sin(phase);
 131                 i += 2;
 132         }
 133
 134         ii = i - 4;
 135         for (k = -(L / 2) + 1; k <= -1; k++)
 136         {
 137                 v[i] = v[ii];
 138                 v[i + 1] = v[ii + 1];
 139                 ii -= 2;
 140                 i += 2;
 141         }
 142
 143         // COMPUTE L-POINT FFT OF V(N)
 144
 145         snfft(v, n3, 1);
 146
 147         // MULTIPLY THE ARRAYS Y AND V, POINT-BY-POINT; CALL THE RESULT G(N)
 148
 149         REAL g[2 * L];   // bug 195 - correct array sizes - swr - 1/17/06
 150         i = 0;
 151         for (k = 0; k < L; k++)
 152         {
 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];
 155                 i += 2;
 156         }
 157
 158         // COMPUTE L-POINT INVERSE FFT OF G(N)
 159
 160         snfft(g, n3, -1);
 161
 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)
 164
 165         i = 2 * nnew;
 166         ii = 0;
 167         for (k = 0; k <= nnew; k++)
 168         {
 169                 rl1 = dold * rl[k];
 170                 im1 = dold * im[k];
 171                 odata[i] = g[ii] * rl1 - g[ii + 1] * im1;
 172                 odata[i + 1] = g[ii + 1] * rl1 + g[ii] * im1;
 173                 i += 2;
 174                 ii += 2;
 175         }
 176
 177         i = 0;
 178         ii = 2 * L - 2 * nnew;
 179         for (k = -nnew; k <= -1; k++)
 180         {
 181                 rl1 = dold * rl[-k];
 182                 im1 = dold * im[-k];
 183                 odata[i] = g[ii] * rl1 - g[ii + 1] * im1;
 184                 odata[i + 1] = g[ii + 1] * rl1 + g[ii] * im1;
 185                 i += 2;
 186                 ii += 2;
 187         }
 188
 189         return;
 190 }