+/* 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);
+ }
+}