/**********
Copyright 1994 Macquarie University, Sydney Australia.  All rights reserved.
Author:   1994 Anthony E. Parker, Department of Electronics, Macquarie Uni.
**********/

/*
 * Code to do fourier transforms on data.
 */

#include "ngspice/ngspice.h"
#include "ngspice/ftedefs.h"
#include "ngspice/dvec.h"
#include "ngspice/sim.h"

#include "spec.h"
#include "parse.h"
#include "variable.h"
#include "ngspice/missing_math.h"
#include "../misc/misc_time.h"


void
com_spec(wordlist *wl)
{
    ngcomplex_t **fdvec = NULL;
    double  **tdvec = NULL;
    double  *freq, *win = NULL, *time, *dc = NULL;
    double  startf, stopf, stepf, span;
    int     fpts, i, j, k, tlen, ngood;
    bool    trace;
    char    *s;
    struct dvec  *f, *vlist, *lv = NULL, *vec;
    struct pnode *pn, *names = NULL;

    if (!plot_cur || !plot_cur->pl_scale) {
        fprintf(cp_err, "Error: no vectors loaded.\n");
        goto done;
    }

    if (!isreal(plot_cur->pl_scale) ||
        ((plot_cur->pl_scale)->v_type != SV_TIME)) {
        fprintf(cp_err, "Error: spec needs real time scale\n");
        goto done;
    }

    s = wl->wl_word;
    tlen = (plot_cur->pl_scale)->v_length;
    if ((freq = ft_numparse(&s, FALSE)) == NULL || (*freq < 0.0)) {
        fprintf(cp_err, "Error: bad start freq %s\n", wl->wl_word);
        goto done;
    }
    startf = *freq;

    wl = wl->wl_next;
    s = wl->wl_word;
    if ((freq = ft_numparse(&s, FALSE)) == NULL || (*freq <= startf)) {
        fprintf(cp_err, "Error: bad stop freq %s\n", wl->wl_word);
        goto done;
    }
    stopf = *freq;

    wl = wl->wl_next;
    s = wl->wl_word;
    if ((freq = ft_numparse(&s, FALSE)) == NULL || !(*freq <= (stopf-startf))) {
        fprintf(cp_err, "Error: bad step freq %s\n", wl->wl_word);
        goto done;
    }
    stepf = *freq;

    wl = wl->wl_next;
    time = (plot_cur->pl_scale)->v_realdata;
    span = time[tlen-1] - time[0];
    if (stopf > 0.5*tlen/span) {
        fprintf(cp_err,
                "Error: nyquist limit exceeded, try stop freq less than %e Hz\n",
                tlen/2/span);
        goto done;
    }
    span = ((int)(span*stepf*1.000000000001))/stepf;
    if (span > 0) {
        startf = (int)(startf/stepf*1.000000000001) * stepf;
        fpts = (int)((stopf - startf)/stepf + 1.);
        if (stopf > startf + (fpts-1)*stepf)
            fpts++;
    } else {
        fprintf(cp_err, "Error: time span limits step freq to %1.1e Hz\n",
                1/(time[tlen-1] - time[0]));
        goto done;
    }
    win = TMALLOC(double, tlen);
    {
        char   window[BSIZE_SP];
        double maxt = time[tlen-1];
        if (!cp_getvar("specwindow", CP_STRING, window, sizeof(window)))
            strcpy(window, "hanning");
        if (eq(window, "none"))
            for (i = 0; i < tlen; i++)
                win[i] = 1;
        else if (eq(window, "rectangular"))
            for (i = 0; i < tlen; i++) {
                if (maxt-time[i] > span) {
                    win[i] = 0;
                } else {
                    win[i] = 1;
                }
            }
        else if (eq(window, "hanning") || eq(window, "cosine"))
            for (i = 0; i < tlen; i++) {
                if (maxt-time[i] > span) {
                    win[i] = 0;
                } else {
                    win[i] = 1 - cos(2*M_PI*(time[i]-maxt)/span);
                }
            }
        else if (eq(window, "hamming"))
            for (i = 0; i < tlen; i++) {
                if (maxt-time[i] > span) {
                    win[i] = 0;
                } else {
                    win[i] = 1 - 0.92/1.08*cos(2*M_PI*(time[i]-maxt)/span);
                }
            }
        else if (eq(window, "triangle") || eq(window, "bartlet"))
            for (i = 0; i < tlen; i++) {
                if (maxt-time[i] > span) {
                    win[i] = 0;
                } else {
                    win[i] = 2 - fabs(2+4*(time[i]-maxt)/span);
                }
            }
        else if (eq(window, "blackman")) {
            int order;
            if (!cp_getvar("specwindoworder", CP_NUM, &order, 0))
                order = 2;
            if (order < 2)      /* only order 2 supported here */
                order = 2;
            for (i = 0; i < tlen; i++) {
                if (maxt-time[i] > span) {
                    win[i] = 0;
                } else {
                    win[i]  = 1;
                    win[i] -= 0.50/0.42*cos(2*M_PI*(time[i]-maxt)/span);
                    win[i] += 0.08/0.42*cos(4*M_PI*(time[i]-maxt)/span);
                }
            }
        } else if (eq(window, "gaussian")) {
            int order;
            double scale;
            if (!cp_getvar("specwindoworder", CP_NUM, &order, 0))
                order = 2;
            if (order < 2)
                order = 2;
            scale = pow(2*M_PI/order, 0.5)*(0.5-erfc(pow(order, 0.5)));
            for (i = 0; i < tlen; i++) {
                if (maxt-time[i] > span) {
                    win[i] = 0;
                } else {
                    win[i] = exp(-0.5*order*(1-2*(maxt-time[i])/span)
                                 *(1-2*(maxt-time[i])/span))/scale;
                }
            }
        } else {
            fprintf(cp_err, "Warning: unknown window type %s\n", window);
            goto done;
        }
    }

    names = ft_getpnames(wl, TRUE);
    vlist = NULL;
    ngood = 0;
    for (pn = names; pn; pn = pn->pn_next) {
        vec = ft_evaluate(pn);
        for (; vec; vec = vec->v_link2) {

            if (vec->v_length != tlen) {
                fprintf(cp_err, "Error: lengths don't match: %d, %d\n",
                        vec->v_length, tlen);
                continue;
            }

            if (!isreal(vec)) {
                fprintf(cp_err, "Error: %s isn't real!\n", vec->v_name);
                continue;
            }

            if (vec->v_type == SV_TIME) {
                continue;
            }

            if (!vlist)
                vlist = vec;
            else
                lv->v_link2 = vec;

            lv = vec;
            ngood++;
        }
    }

    if (!ngood)
        goto done;

    plot_cur = plot_alloc("spectrum");
    plot_cur->pl_next = plot_list;
    plot_list = plot_cur;
    plot_cur->pl_title = copy((plot_cur->pl_next)->pl_title);
    plot_cur->pl_name = copy("Spectrum");
    plot_cur->pl_date = copy(datestring());

    f = dvec_alloc(copy("frequency"),
                   SV_FREQUENCY,
                   VF_REAL | VF_PERMANENT | VF_PRINT,
                   fpts, NULL);
    vec_new(f);
    freq = f->v_realdata;

    tdvec = TMALLOC(double  *, ngood);
    fdvec = TMALLOC(ngcomplex_t *, ngood);
    for (i = 0, vec = vlist; i < ngood; i++) {
        tdvec[i] = vec->v_realdata;
        f = dvec_alloc(vec_basename(vec),
                       vec->v_type,
                       VF_COMPLEX | VF_PERMANENT,
                       fpts, NULL);
        vec_new(f);
        fdvec[i] = f->v_compdata;
        vec = vec->v_link2;
    }

    dc = TMALLOC(double, ngood);
    for (i = 0; i < ngood; i++)
        dc[i] = 0;

    for (k = 1; k < tlen; k++) {
        double amp = win[k]/(tlen-1);
        for (i = 0; i < ngood; i++) {
            dc[i] += tdvec[i][k]*amp;
        }
    }
    trace = cp_getvar("spectrace", CP_BOOL, NULL, 0);
    for (j = (startf == 0 ? 1 : 0); j < fpts; j++) {
        freq[j] = startf + j*stepf;
        if (trace)
            fprintf(cp_err, "spec: %e Hz: \r", freq[j]);
        for (i = 0; i < ngood; i++) {
            fdvec[i][j].cx_real = 0;
            fdvec[i][j].cx_imag = 0;
        }
        for (k = 1; k < tlen; k++) {
            double
                amp = 2*win[k]/(tlen-1),
                rad = 2*M_PI*time[k]*freq[j],
                cosa = amp*cos(rad),
                sina = amp*sin(rad);
            for (i = 0; i < ngood; i++) {
                double value = tdvec[i][k]-dc[i];
                fdvec[i][j].cx_real += value*cosa;
                fdvec[i][j].cx_imag += value*sina;
            }
        }
#ifdef HAS_PROGREP
        SetAnalyse("spec", (int)(j * 1000./ fpts));
#endif
    }

    if (startf == 0) {
        freq[0] = 0;
        for (i = 0; i < ngood; i++) {
            fdvec[i][0].cx_real = dc[i];
            fdvec[i][0].cx_imag = 0;
        }
    }

    if (trace)
        fprintf(cp_err, "                           \r");

#ifdef KEEPWINDOW
        f = dvec_alloc(copy("win"),
                       SV_NOTYPE,
                       VF_REAL | VF_PERMANENT,
                       tlen, win);
        win = NULL;
        vec_new(f);
#endif

done:
    tfree(dc);

    tfree(tdvec);
    tfree(fdvec);
    tfree(win);

    free_pnode(names);
}