# 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;
+}
+
