#  spice3f4/3f5 patch level 2
#  Adds spec command for spectrum analysis using DFT
#  July 11, 1994
#  Anthony Parker (tonyp@mpce.mq.edu.au)
#  Macquarie University, Sydney Australia 2109 
#
diff -ruN src.orig/include/fteext.h src/include/fteext.h
--- src.orig/include/fteext.h	Tue Jul 27 14:51:18 1993
+++ src/include/fteext.h	Sun May 18 13:25:52 2003
@@ -208,6 +208,10 @@
 
 extern void com_fourier();
 
+/* spec.c */
+
+extern void com_spec();
+
 /* ginterface.c */
 
 extern bool gi_init();
diff -ruN src.orig/lib/fte/makedefs src/lib/fte/makedefs
--- src.orig/lib/fte/makedefs	Mon Mar 22 20:56:10 1993
+++ src/lib/fte/makedefs	Sun May 18 13:25:44 2003
@@ -11,7 +11,8 @@
 		  parse.c plot5.c plotcurv.c points.c postcoms.c postsc.c \
 		  rawfile.c resource.c runcoms.c shyu.c signal.c spcmdtab.c \
 		  spiceif.c subckt.c types.c vectors.c where.c x10.c x11.c \
-		  gens.c newcoms.c dimens.c xgraph.c runcoms2.c breakp2.c
+		  gens.c newcoms.c dimens.c xgraph.c runcoms2.c breakp2.c \
+                  spec.c
 
 COBJS		= agraf.o arg.o aspice.o breakp.o circuits.o clip.o cmath1.o \
 		  cmath2.o cmath3.o cmath4.o compose.o cpitf.o debugcom.o \
@@ -22,7 +23,8 @@
 		  parse.o plot5.o plotcurv.o points.o postcoms.o postsc.o \
 		  rawfile.o resource.o runcoms.o shyu.o signal.o spcmdtab.o \
 		  spiceif.o subckt.o types.o vectors.o where.o x10.o x11.o \
-		  gens.o newcoms.o dimens.o xgraph.o runcoms2.o breakp2.o
+		  gens.o newcoms.o dimens.o xgraph.o runcoms2.o breakp2.o \
+                  spec.o
 
 LIBRARY		= fte
 LIB_TARGET	= $(OBJLIB_DIR)/$(LIBRARY).a
@@ -92,3 +94,4 @@
 xgraph.o: xgraph.c
 runcoms2.o: runcoms2.c
 breakp2.o: breakp2.c
+spec.o: spec.c
diff -ruN src.orig/lib/fte/msc51.bat src/lib/fte/msc51.bat
--- src.orig/lib/fte/msc51.bat	Fri Jul 30 04:26:44 1993
+++ src/lib/fte/msc51.bat	Sun May 18 13:25:44 2003
@@ -61,4 +61,5 @@
 cl /I..\..\include /c xgraph.c >> ..\..\msc.out
 cl /I..\..\include /c runcoms2.c >> ..\..\msc.out
 cl /I..\..\include /c breakp2.c >> ..\..\msc.out
+cl /I..\..\include /c spec.c >> ..\..\msc.out
 lib ..\fte.lib @response.lib >> ..\..\msc.out
diff -ruN src.orig/lib/fte/response.lib src/lib/fte/response.lib
--- src.orig/lib/fte/response.lib	Fri Jul 30 04:26:45 1993
+++ src/lib/fte/response.lib	Sun May 18 13:25:44 2003
@@ -59,4 +59,5 @@
 +dimens.obj&
 +xgraph.obj&
 +runcoms2.obj&
++spec.obj&
 +breakp2.obj;
diff -ruN src.orig/lib/fte/spcmdtab.c src/lib/fte/spcmdtab.c
--- src.orig/lib/fte/spcmdtab.c	Thu Jul 29 20:49:33 1993
+++ src/lib/fte/spcmdtab.c	Sun May 18 13:25:44 2003
@@ -162,6 +162,10 @@
       { 0, 040000, 040000, 040000 }, E_DEFHMASK, 1, LOTS,
       (int (*)()) NULL,
       "fund_freq vector ... : Do a fourier analysis of some data." } ,
+    { "spec", com_spec, false, false, true,
+      { 0, 0, 0, 0 }, E_DEFHMASK, 4, LOTS,
+      (int (*)()) NULL,
+      "start_freq stop_freq step_freq vector ... : Create a frequency domain plot." } ,
     { "show", com_show, false, true, false,
       { 040, 040, 040, 040 }, E_DEFHMASK, 0, LOTS,
       (int (*)()) NULL,
diff -ruN src.orig/lib/fte/spec.c src/lib/fte/spec.c
--- src.orig/lib/fte/spec.c	Wed Dec 31 19:00:00 1969
+++ src/lib/fte/spec.c	Sun May 18 13:25:44 2003
@@ -0,0 +1,286 @@
+/**********
+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 "spice.h"
+#include "ftedefs.h"
+#include "ftedata.h"
+#include "util.h"
+#include "suffix.h"
+
+void
+com_spec(wl)
+    wordlist *wl;
+{
+    complex **fdvec;
+    double  **tdvec;
+    double  *freq, *win, *time, *dc;
+    double  startf, stopf, stepf, span;
+    int     fpts, i, j, k, tlen, ngood;
+    bool    trace;
+    char    *s;
+    struct dvec  *f, *vlist, *lv, *vec;
+    struct pnode *names, *first_name;
+
+    if (!plot_cur || !plot_cur->pl_scale) {
+        fprintf(cp_err, "Error: no vectors loaded.\n");
+        return;
+    }
+    if (!isreal(plot_cur->pl_scale) || 
+        ((plot_cur->pl_scale)->v_type != SV_TIME)) {
+        fprintf(cp_err, "Error: spec needs real time scale\n");
+        return;
+    }
+    s = wl->wl_word;
+    tlen = (plot_cur->pl_scale)->v_length;
+    if (!(freq = ft_numparse(&s, false)) || (*freq < 0.0)) {
+        fprintf(cp_err, "Error: bad start freq %s\n", wl->wl_word);
+        return;
+    }
+    startf = *freq;
+    wl = wl->wl_next;
+    s = wl->wl_word;
+    if (!(freq = ft_numparse(&s, false)) || (*freq <= startf)) {
+        fprintf(cp_err, "Error: bad stop freq %s\n", wl->wl_word);
+        return;
+    }
+    stopf = *freq;
+    wl = wl->wl_next;
+    s = wl->wl_word;
+    if (!(freq = ft_numparse(&s, false)) || !(*freq <= (stopf-startf))) {
+        fprintf(cp_err, "Error: bad step freq %s\n", wl->wl_word);
+        return;
+    }
+    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);
+        return;
+    }
+    span = ((int)(span*stepf*1.000000000001))/stepf;
+    if (span > 0) {
+        startf = (int)(startf/stepf*1.000000000001) * stepf;
+        fpts = (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]));
+        return;
+    }
+    win = (double *) tmalloc(tlen * sizeof (double));
+    {
+       char   window[BSIZE_SP];
+       double maxt = time[tlen-1];
+       if (!cp_getvar("specwindow", VT_STRING, 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", VT_NUM, &order)) order = 2;
+          if (order < 2) order = 2;  /* only order 2 supported here */
+          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;
+          extern double erfc();
+          if (!cp_getvar("specwindoworder", VT_NUM, &order)) 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);
+          tfree(win);
+          return;
+       }
+    }
+    
+    names = ft_getpnames(wl, true);
+    first_name = names;
+    vlist = NULL;
+    ngood = 0;
+    while (names) {
+        vec = ft_evaluate(names);
+        names = names->pn_next;
+        while (vec) {
+            if (vec->v_length != tlen) {
+                fprintf(cp_err, "Error: lengths don't match: %d, %d\n",
+                        vec->v_length, tlen);
+                vec = vec->v_link2;
+                continue;
+            }
+            if (!isreal(vec)) {
+                fprintf(cp_err, "Error: %s isn't real!\n", 
+                        vec->v_name);
+                vec = vec->v_link2;
+                continue;
+            }
+            if (vec->v_type == SV_TIME) {
+                vec = vec->v_link2;
+                continue;
+            }
+	        if (!vlist)
+	           vlist = vec;
+	        else
+	           lv->v_link2 = vec;
+            lv = vec;
+            vec = vec->v_link2;
+            ngood++;
+        }
+    }
+    free_pnode(first_name);
+    if (!ngood) {
+       tfree(win);
+       return;
+    }
+ 
+    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( ));
+
+    freq = (double *) tmalloc(fpts * sizeof(double));
+    f = alloc(struct dvec);
+    ZERO(f, struct dvec);
+    f->v_name = copy("frequency");
+    f->v_type = SV_FREQUENCY;
+    f->v_flags = (VF_REAL | VF_PERMANENT | VF_PRINT);
+    f->v_length = fpts;
+    f->v_realdata = freq;
+    vec_new(f);
+
+    tdvec = (double  **) tmalloc(ngood * sizeof(double  *));
+    fdvec = (complex **) tmalloc(ngood * sizeof(complex *));
+    for (i = 0, vec = vlist; i<ngood; i++) {
+       tdvec[i] = vec->v_realdata;
+       fdvec[i] = (complex *) tmalloc(fpts * sizeof(complex));
+       f = alloc(struct dvec);
+       ZERO(f, struct dvec);
+       f->v_name = vec_basename(vec);
+       f->v_type = vec->v_type;
+       f->v_flags = (VF_COMPLEX | VF_PERMANENT);
+       f->v_length = fpts;
+       f->v_compdata = fdvec[i];
+       vec_new(f);
+       vec = vec->v_link2;
+    }
+
+    dc = (double  *) tmalloc(ngood * sizeof(double));
+    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;
+       }
+    }
+    cp_getvar("spectrace", VT_BOOL, &trace);
+    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;
+          }
+       }
+    }
+    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");
+    tfree(dc);
+    tfree(tdvec);
+    tfree(fdvec);
+
+#ifdef KEEPWINDOW
+    f = alloc(struct dvec);
+    ZERO(f, struct dvec);
+    f->v_name = copy("win");
+    f->v_type = SV_NOTYPE;
+    f->v_flags = (VF_REAL | VF_PERMANENT);
+    f->v_length = tlen;
+    f->v_realdata = win;
+    vec_new(f);
+#else
+    tfree(win);
+#endif
+	
+    return;
+}
+