Major changes. Added metering and power spectrum, other fixes. Rearranged headers...

[dttsp.git] / jDttSP / spectrum.c

diff --git a/jDttSP/spectrum.c b/jDttSP/spectrum.c
new file mode 100644 (file)
index 0000000..441f92e
--- /dev/null
+++ b/jDttSP/spectrum.c
@@ -0,0 +1,107 @@
+/* spectrum.c */
+/*
+This file is part of a program that implements a Software-Defined Radio.
+
+Copyright (C) 2004 by Frank Brickle, AB2KT and Bob McGwier, N4HY
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+The authors can be reached by email at
+
+ab2kt@arrl.net
+or
+rwmcgwier@comcast.net
+
+or by paper mail at
+
+The DTTS Microwave Society
+6 Kathleen Place
+Bridgewater, NJ 08807
+*/  
+
+#include <spectrum.h>
+
+// snapshot of current signal
+void
+snap_spectrum(SpecBlock *sb, int label) {
+  int i, j;
+
+  // where most recent signal started
+  j = (sb->fill - sb->buflen + sb->size) % sb->size;
+
+  // 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]);
+    j = ++j % sb->size;
+  }
+  
+  sb->label = label;
+}
+
+// snapshot -> frequency domain
+void
+compute_spectrum(SpecBlock *sb) {
+  int i, half = sb->size / 2;
+
+  // assume timebuf has windowed current snapshot
+
+  fftw_one(sb->plan,
+          (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);
+  }
+}
+
+void
+init_spectrum(SpecBlock *sb) {
+  sb->fill = sb->size - sb->buflen;
+  sb->accum = newCXB(sb->size, 0, "spectrum accum");
+  sb->timebuf = newCXB(sb->size, 0, "spectrum timebuf");
+  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->plan = fftw_create_plan(sb->size, FFTW_FORWARD, sb->planbits);
+}
+
+void
+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));
+}
+
+void
+finish_spectrum(SpecBlock *sb) {
+  if (sb) {
+    delCXB(sb->accum);
+    delCXB(sb->timebuf);
+    delCXB(sb->freqbuf);
+    delvec_REAL(sb->window);
+    safefree((char *) sb->output);
+    fftw_destroy_plan(sb->plan);
+  }
+}