Major update

[dttsp.git] / jDttSP / spectrum.c

diff --git a/jDttSP/spectrum.c b/jDttSP/spectrum.c
index 441f92e16a7f58f40428bfd563f246111537ed3f..e547567b7ade4b561f4ddc081c6cdd4cab4c661e 100644 (file)
--- a/jDttSP/spectrum.c
+++ b/jDttSP/spectrum.c
@@ -39,12 +39,14 @@ snap_spectrum(SpecBlock *sb, int label) {
   int i, j;
 
   // where most recent signal started
-  j = (sb->fill - sb->buflen + sb->size) % sb->size;
+  j = sb->fill;
 
   // copy starting from there in circular fashion,
   // applying window as we go
   for (i = 0; i < sb->size; i++) {
-    CXBdata(sb->timebuf, i) = Cscl(CXBdata(sb->timebuf, j), sb->window[i]);
+    COMPLEX z = CXBdata(sb->accum, j);
+    sb->oscope[i] = (float) z.re;
+    CXBdata(sb->timebuf, i) = Cscl(z, sb->window[i]);
     j = ++j % sb->size;
   }
   
@@ -62,18 +64,17 @@ compute_spectrum(SpecBlock *sb) {
           (fftw_complex *) CXBbase(sb->timebuf),
           (fftw_complex *) CXBbase(sb->freqbuf));
 
-  if (sb->scale == SPEC_MAG) {
-    for (i = 0; i < half; i++)
-      sb->output[i + half] = Cmag(CXBdata(sb->freqbuf, i));
-    for (; i < sb->size; i++)
-      sb->output[i] = Cmag(CXBdata(sb->freqbuf, i));
-
-  } else { // SPEC_PWR
-    for (i = 0; i < half; i++)
-      sb->output[i + half] = 10.0 * log10(Csqrmag(CXBdata(sb->freqbuf, i)) + 1e-60);
-    for (; i < sb->size; i++)
-      sb->output[i] = 10.0 * log10(Csqrmag(CXBdata(sb->freqbuf, i)) + 1e-60);
-  }
+  if (sb->scale == SPEC_MAG)
+
+    for (i = 0; i < sb->size; i++)
+      sb->output[i] =
+       (float) Cmag(CXBdata(sb->freqbuf, (i + half) % sb->size));
+
+  else // SPEC_PWR
+
+    for (i = 0; i < sb->size; i++)
+      sb->output[i] =
+       (float) (10.0 * log10(Csqrmag(CXBdata(sb->freqbuf, (i + half) % sb->size)) + 1e-60));
 }
 
 void
@@ -84,6 +85,7 @@ init_spectrum(SpecBlock *sb) {
   sb->freqbuf = newCXB(sb->size, 0, "spectrum freqbuf");
   sb->window = newvec_REAL(sb->size, "spectrum window");
   sb->output = (float *) safealloc(sb->size, sizeof(float), "spectrum output");
+  sb->oscope = (float *) safealloc(sb->size, sizeof(float), "spectrum oscope");
   sb->plan = fftw_create_plan(sb->size, FFTW_FORWARD, sb->planbits);
 }
 
@@ -92,6 +94,7 @@ reinit_spectrum(SpecBlock *sb) {
   sb->fill = sb->size - sb->buflen;
   memset((char *) CXBbase(sb->accum), 0, sb->size * sizeof(REAL));
   memset((char *) sb->output, 0, sb->size * sizeof(float));
+  memset((char *) sb->oscope, 0, sb->size * sizeof(float));
 }
 
 void
@@ -102,6 +105,7 @@ finish_spectrum(SpecBlock *sb) {
     delCXB(sb->freqbuf);
     delvec_REAL(sb->window);
     safefree((char *) sb->output);
+    safefree((char *) sb->oscope);
     fftw_destroy_plan(sb->plan);
   }
 }