spectrum.c

/* spectrum.c */
/*
This file is part of a program that implements a Software-Defined Radio.

Copyright (C) 2004, 2005, 2006 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;

    // copy starting from there in circular fashion,
    // applying window as we go
    if(!sb->polyphase)
    {
        for (i = 0; i < sb->size; i++)
        {
            CXBdata(sb->timebuf, i) =
                Cscl(CXBdata(sb->accum, j), sb->window[i]);
            j = (++j & sb->mask);
        }
    } else
    {
        int k;
        for (i = 0; i < sb->size; i++)
        {
            CXBreal(sb->timebuf, i) = CXBreal(sb->accum, j) * sb->window[i];
            CXBimag(sb->timebuf, i) = CXBimag(sb->accum, j) * sb->window[i];
            for (k = 1; k < 8; k++)
            {
                int accumidx = (j + k * sb->size) & sb->mask;
                int winidx = i + k * sb->size;
                CXBreal(sb->timebuf, i) +=
                    CXBreal(sb->accum, accumidx) * sb->window[winidx];
                CXBimag(sb->timebuf, i) +=
                    CXBimag(sb->accum, accumidx) * sb->window[winidx];
            }
            j = (++j & sb->mask);
        }

    }
    sb->label = label;
}

void
snap_scope(SpecBlock *sb, int label)
{
    int i, j;

    // where most recent signal started
    j = sb->fill;

    // copy starting from there in circular fashion
    for (i = 0; i < sb->size; i++)
    {
        CXBdata(sb->timebuf, i) = CXBdata(sb->accum, j);
        j = (++j & sb->mask);
    }

    sb->label = label;
}

// snapshot -> frequency domain
void
compute_spectrum(SpecBlock *sb)
{
    int i, j, half = sb->size / 2;

    // assume timebuf has windowed current snapshot

    fftwf_execute(sb->plan);

    if(sb->scale == SPEC_MAG)
    {
        for (i = 0, j = half; i < half; i++, j++)
        {
            sb->output[i] = (float)Cmag(CXBdata(sb->freqbuf, j));
            sb->output[j] = (float)Cmag(CXBdata(sb->freqbuf, i));
        }
    } else
    {               // SPEC_PWR
        for (i = 0, j = half; i < half; i++, j++)
        {
            sb->output[i] =
                (float)(10.0 *
                        log10(Csqrmag(CXBdata(sb->freqbuf, j)) + 1e-60));
            sb->output[j] =
                (float)(10.0 *
                        log10(Csqrmag(CXBdata(sb->freqbuf, i)) + 1e-60));
        }
    }
}

void
init_spectrum(SpecBlock *sb)
{
    sb->fill = 0;
    sb->accum = newCXB(sb->size * 16, 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 * 16, "spectrum window");
    makewindow(BLACKMANHARRIS_WINDOW, sb->size, sb->window);
    sb->mask = sb->size - 1;
    sb->polyphase = FALSE;
    sb->output =
        (float *)safealloc(sb->size, sizeof(float), "spectrum output");
    sb->plan =
        fftwf_plan_dft_1d(sb->size, (fftwf_complex *)CXBbase(sb->timebuf),
                          (fftwf_complex *)CXBbase(sb->freqbuf),
                          FFTW_FORWARD, sb->planbits);
}

void
reinit_spectrum(SpecBlock *sb)
{
    size_t polysize = 1;
    sb->fill = 0;
    if(sb->polyphase) polysize = 8;
    memset((char *)CXBbase(sb->accum), 0,
           polysize * 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);
        fftwf_destroy_plan(sb->plan);
    }
}