jDttSP/win/power_spectrum.c

   1 #include <common.h>
   2
   3 #ifdef _WINDOWS
   4
   5 PWS
   6 newPowerSpectrum(int size, SPECTRUMTYPE type) {
   7   PWS p = (PWS) safealloc(1, sizeof(powerspectrum), "New Power Spectrum");
   8
   9   p->buf.win = newRLB(size, NULL, "Power Spectrum Window Buffer");
  10   p->buf.res = newRLB(32768, NULL, "Save buffer");
  11   p->buf.wdat = newCXB(size, NULL, "Power Spectrum Windowed Signal buffer");
  12   p->buf.zdat = newCXB(size, NULL, "Power Spectrum FFT Results Buffer");
  13   p->fft.size = size;
  14   p->flav.spec = type;
  15
  16   makewindow(BLACKMANHARRIS_WINDOW, size, RLBbase(p->buf.win));
  17
  18 #ifdef notdef
  19   {
  20     REAL *tmp = (REAL *) alloca(size * sizeof(REAL));
  21     int n = size / 2, nb = n * sizeof(REAL);
  22     memcpy(&tmp[n], RLBbase(p->buf.win), nb);
  23     memcpy(tmp, &RLBdata(p->buf.win, n), nb);
  24   }
  25 #endif
  26
  27   p->fft.plan = fftw_create_plan(size, FFTW_FORWARD, uni.wisdom.bits);
  28
  29   return p;
  30 }
  31
  32 void
  33 delPowerSpectrum(PWS pws) {
  34   delCXB(pws->buf.wdat);
  35   delCXB(pws->buf.zdat);
  36   delRLB(pws->buf.win);
  37   delRLB(pws->buf.win2);
  38   fftw_destroy_plan(pws->fft.plan);
  39   safefree((char *) pws);
  40 }
  41
  42 DttSP_EXP BOOLEAN
  43 process_spectrum(float *results) {
  44   int i, j, n = CXBsize(uni.powsp.buf);
  45
  46   if (getChan_nowait(uni.powsp.chan.c,
  47                      (char *) CXBbase(uni.powsp.buf),
  48                      uni.powsp.fft.size * sizeof(COMPLEX))
  49       == TRUE) {
  50
  51     // Window the data to minimize "splatter" in the powerspectrum
  52     for (i = 0; i < uni.powsp.fft.size; i++)
  53       CXBdata(uni.powsp.gen->buf.wdat, i) = Cscl(CXBdata(uni.powsp.buf, i),
  54                                                  RLBdata(uni.powsp.gen->buf.win, i));
  55
  56     // Compute the FFT
  57     fftw_one(uni.powsp.gen->fft.plan,
  58              (fftw_complex *) CXBbase(uni.powsp.gen->buf.wdat),
  59              (fftw_complex *) CXBbase(uni.powsp.gen->buf.zdat));
  60     for (i = 0; i < n / 2; i++) {
  61       RLBdata(uni.powsp.gen->buf.res, 8 * i) =
  62         results[8 * i] =
  63         (float) (10.0 *
  64                  log10(1e-160 +
  65                        Csqrmag(CXBdata(uni.powsp.gen->buf.zdat,
  66                                        (i + uni.buflen)))));
  67       RLBdata(uni.powsp.gen->buf.res, 8 * (i + uni.buflen)) =
  68         results[8 * (i + uni.buflen)] =
  69         (float) (10.0 *
  70                   log10(1e-160 +
  71                         Csqrmag(CXBdata(uni.powsp.gen->buf.zdat, i))));
  72     }
  73
  74     // Interpolate for display on high resolution monitors
  75
  76     for (i = 0; i < 32768; i += 8)
  77       for (j = 1; j < 8; j++)
  78         RLBdata(uni.powsp.gen->buf.res, i + j) =
  79           results[i + j] = (float) 0.125 * (results[i] * (float) (8 - j)
  80                                             + results[i + 8] * (float) j);
  81     return TRUE;
  82   } else {
  83     for (i = 0; i < 32768; i++)
  84       results[i] = (float) RLBdata(uni.powsp.gen->buf.res, i);
  85     return FALSE;
  86   }
  87 }
  88
  89 DttSP_EXP BOOLEAN
  90 process_phase(float *results, int numpoints) {
  91   int i, j;
  92   if (getChan_nowait(uni.powsp.chan.c,
  93                      (char *) CXBbase(uni.powsp.buf),
  94                      uni.powsp.fft.size * sizeof(COMPLEX))
  95       == TRUE) {
  96     for (i = 0, j = 0; i < numpoints; i++, j += 2) {
  97       RLBdata(uni.powsp.gen->buf.res, j) =
  98         results[j] = (float) CXBreal(uni.powsp.buf, i);
  99       RLBdata(uni.powsp.gen->buf.res, j + 1) =
 100         results[j + 1] = (float) CXBimag(uni.powsp.buf, i);
 101     } return TRUE;
 102   } else {
 103     for (i = 0; i < 2 * numpoints; i++)
 104       results[i] = (float) RLBdata(uni.powsp.gen->buf.res, i);
 105     return FALSE;
 106   }
 107 }
 108
 109 // Send out Real Signal, reals only
 110
 111 DttSP_EXP BOOLEAN
 112 process_scope(float *results) {
 113   int i;
 114   if (getChan_nowait(uni.powsp.chan.c,
 115                      (char *) CXBbase(uni.powsp.buf),
 116                      uni.powsp.fft.size * sizeof(COMPLEX))
 117       == TRUE) {
 118     for (i = 0; i < uni.powsp.fft.size; i++) {
 119       results[i] = (float) CXBreal(uni.powsp.buf, i);
 120     } return TRUE;
 121   } else {
 122     for (i = 0; i < uni.powsp.fft.size; i++)
 123       results[i] = (float) CXBreal(uni.powsp.buf, i);
 124     return FALSE;
 125   }
 126 }
 127
 128 #endif