jDttSP/spectrum.c

   1 /* spectrum.c */
   2 /*
   3 This file is part of a program that implements a Software-Defined Radio.
   4
   5 Copyright (C) 2004 by Frank Brickle, AB2KT and Bob McGwier, N4HY
   6
   7 This program is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or
  10 (at your option) any later version.
  11
  12 This program is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with this program; if not, write to the Free Software
  19 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  20
  21 The authors can be reached by email at
  22
  23 ab2kt@arrl.net
  24 or
  25 rwmcgwier@comcast.net
  26
  27 or by paper mail at
  28
  29 The DTTS Microwave Society
  30 6 Kathleen Place
  31 Bridgewater, NJ 08807
  32 */
  33
  34 #include <spectrum.h>
  35
  36 // snapshot of current signal
  37 void
  38 snap_spectrum(SpecBlock *sb, int label) {
  39   int i, j;
  40
  41   // where most recent signal started
  42   j = sb->fill;
  43
  44   // copy starting from there in circular fashion,
  45   // applying window as we go
  46   for (i = 0; i < sb->size; i++) {
  47     COMPLEX z = CXBdata(sb->accum, j);
  48     sb->oscope[i] = (float) z.re;
  49     CXBdata(sb->timebuf, i) = Cscl(z, sb->window[i]);
  50     j = ++j % sb->size;
  51   }
  52
  53   sb->label = label;
  54 }
  55
  56 // snapshot -> frequency domain
  57 void
  58 compute_spectrum(SpecBlock *sb) {
  59   int i, half = sb->size / 2;
  60
  61   // assume timebuf has windowed current snapshot
  62
  63   fftw_one(sb->plan,
  64            (fftw_complex *) CXBbase(sb->timebuf),
  65            (fftw_complex *) CXBbase(sb->freqbuf));
  66
  67   if (sb->scale == SPEC_MAG)
  68
  69     for (i = 0; i < sb->size; i++)
  70       sb->output[i] =
  71         (float) Cmag(CXBdata(sb->freqbuf, (i + half) % sb->size));
  72
  73   else  // SPEC_PWR
  74
  75     for (i = 0; i < sb->size; i++)
  76       sb->output[i] =
  77         (float) (10.0 * log10(Csqrmag(CXBdata(sb->freqbuf, (i + half) % sb->size)) + 1e-60));
  78 }
  79
  80 void
  81 init_spectrum(SpecBlock *sb) {
  82   sb->fill = sb->size - sb->buflen;
  83   sb->accum = newCXB(sb->size, 0, "spectrum accum");
  84   sb->timebuf = newCXB(sb->size, 0, "spectrum timebuf");
  85   sb->freqbuf = newCXB(sb->size, 0, "spectrum freqbuf");
  86   sb->window = newvec_REAL(sb->size, "spectrum window");
  87   sb->output = (float *) safealloc(sb->size, sizeof(float), "spectrum output");
  88   sb->oscope = (float *) safealloc(sb->size, sizeof(float), "spectrum oscope");
  89   sb->plan = fftw_create_plan(sb->size, FFTW_FORWARD, sb->planbits);
  90 }
  91
  92 void
  93 reinit_spectrum(SpecBlock *sb) {
  94   sb->fill = sb->size - sb->buflen;
  95   memset((char *) CXBbase(sb->accum), 0, sb->size * sizeof(REAL));
  96   memset((char *) sb->output, 0, sb->size * sizeof(float));
  97   memset((char *) sb->oscope, 0, sb->size * sizeof(float));
  98 }
  99
 100 void
 101 finish_spectrum(SpecBlock *sb) {
 102   if (sb) {
 103     delCXB(sb->accum);
 104     delCXB(sb->timebuf);
 105     delCXB(sb->freqbuf);
 106     delvec_REAL(sb->window);
 107     safefree((char *) sb->output);
 108     safefree((char *) sb->oscope);
 109     fftw_destroy_plan(sb->plan);
 110   }
 111 }