1 /***************************************************************************
2 JSPICE3 adaptation of Spice3e2 - Copyright (c) Stephen R. Whiteley 1992
3 Copyright 1990 Regents of the University of California. All rights reserved.
4 Authors: 1985 Wayne A. Christopher
5 1992 Stephen R. Whiteley
6 ****************************************************************************/
7
8 /*
9 * Code to do fourier transforms on data. Note that we do interpolation
10 * to get a uniform grid. Note that if polydegree is 0 then no interpolation
11 * is done.
12 */
13
14 #include "spice.h"
15 #include "ftedefs.h"
16 #include "ftegraph.h"
17 #include "plotext.h"
18 #include "sperror.h"
19 #include "const.h"
20 #include "util.h"
21
22 static int CKTfour();
23
24 #define DEF_FOURGRIDSIZE 200
25
26 /* CKTfour(ndata,numFreq,thd,Time,Value,FundFreq,Freq,Mag,Phase,nMag,nPhase)
27 * len 10 ? inp inp inp out out out out out
28 */
29
30 void
com_fourier(wl)31 com_fourier(wl)
32 wordlist *wl;
33 {
34 struct dvec *time, *vec;
35 struct dvlist *dl0, *dl;
36 struct pnode *names;
37 double *ff, fundfreq, *dp, *stuff;
38 int nfreqs, fourgridsize, polydegree;
39 double *freq, *mag, *phase, *nmag, *nphase; /* Outputs from CKTfour */
40 double thd, *timescale, *grid, d;
41 char *s;
42 int i, err, fw;
43 char xbuf[20];
44 int shift;
45 extern char *kw_nfreqs, *kw_polydegree, *kw_fourgridsize;
46
47 sprintf(xbuf, "%1.1e", 0.0);
48 shift = strlen(xbuf) - 7;
49 if (!plot_cur || !plot_cur->pl_scale) {
50 fprintf(cp_err, "Error: no vectors loaded.\n");
51 return;
52 }
53
54 if ((!cp_getvar(kw_nfreqs, VT_NUM, (char *) &nfreqs)) || (nfreqs < 1))
55 nfreqs = 10;
56 if ((!cp_getvar(kw_polydegree, VT_NUM, (char *) &polydegree)) ||
57 (polydegree < 0))
58 polydegree = 1;
59 if ((!cp_getvar(kw_fourgridsize, VT_NUM, (char *) &fourgridsize)) ||
60 (fourgridsize < 1))
61 fourgridsize = DEF_FOURGRIDSIZE;
62
63 time = plot_cur->pl_scale;
64 if (!isreal(time)) {
65 fprintf(cp_err, "Error: fourier needs real time scale\n");
66 return;
67 }
68 s = wl->wl_word;
69 if (!(ff = ft_numparse(&s, false)) || (*ff <= 0.0)) {
70 fprintf(cp_err, "Error: bad fund freq %s\n", wl->wl_word);
71 return;
72 }
73 fundfreq = *ff;
74
75 wl = wl->wl_next;
76 if ((names = ft_getpnames(wl, true)) == NULL)
77 return;
78 if ((dl0 = ft_dvlist(names)) == NULL)
79 return;
80
81 if (polydegree) {
82 /* Build the grid... */
83 dp = ft_minmax(time, true);
84
85 /* Now get the last fund freq... */
86 d = 1 / fundfreq; /* The wavelength... */
87 if (dp[1] - dp[0] < d) {
88 fprintf(cp_err,
89 "Error: wavelength longer than time span\n");
90 return;
91 }
92 else if (dp[1] - dp[0] > d) {
93 dp[0] = dp[1] - d;
94 }
95 }
96
97 freq = (double *) tmalloc(nfreqs * sizeof (double));
98 mag = (double *) tmalloc(nfreqs * sizeof (double));
99 phase = (double *) tmalloc(nfreqs * sizeof (double));
100 nmag = (double *) tmalloc(nfreqs * sizeof (double));
101 nphase = (double *) tmalloc(nfreqs * sizeof (double));
102 if (polydegree) {
103 grid = (double *) tmalloc(fourgridsize * sizeof (double));
104 stuff = (double *) tmalloc(fourgridsize * sizeof (double));
105
106 d = (dp[1] - dp[0]) / fourgridsize;
107 for (i = 0; i < fourgridsize; i++)
108 grid[i] = dp[0] + i * d;
109 timescale = grid;
110 }
111 else
112 timescale = time->v_realdata;
113
114
115 for (dl = dl0; dl; dl = dl->dl_next) {
116 vec = dl->dl_dvec;
117 if (vec->v_length != time->v_length) {
118 fprintf(cp_err,
119 "Error: lengths don't match: %d, %d\n",
120 vec->v_length, time->v_length);
121 continue;
122 }
123 if (!isreal(vec)) {
124 fprintf(cp_err, "Error: %s isn't real!\n",
125 vec->v_name);
126 continue;
127 }
128
129 if (polydegree) {
130
131 /* Now interpolate the stuff... */
132 if (!ft_interpolate(vec->v_realdata, stuff,
133 time->v_realdata, vec->v_length,
134 grid, fourgridsize, polydegree)) {
135 fprintf(cp_err,
136 "Error: can't interpolate\n");
137 break;
138 }
139 }
140 else {
141 fourgridsize = vec->v_length;
142 stuff = vec->v_realdata;
143 }
144
145 err = CKTfour(fourgridsize, nfreqs, &thd, timescale,
146 stuff, fundfreq, freq, mag, phase, nmag, nphase);
147 if (err != OK) {
148 ft_sperror(err, "fourier");
149 break;
150 }
151
152 out_printf("Fourier analysis for %s:\n", vec->v_name);
153 out_printf(
154 " No. Harmonics: %d, THD: %g %%, Gridsize: %d, Interpolation Degree: %d\n\n",
155 nfreqs, thd, fourgridsize, polydegree);
156 /* Each field will have width cp_numdgt + 6 (or 7
157 * with HP-UX) + 1 if there is a - sign.
158 */
159 fw = ((cp_numdgt > 0) ? cp_numdgt : 6) + 5 + shift;
160 out_printf( "Harmonic %-*s %-*s %-*s %-*s %-*s\n",
161 fw, "Frequency", fw, "Magnitude",
162 fw, "Phase", fw, "Norm. Mag",
163 fw, "Norm. Phase");
164 out_printf( "-------- %-*s %-*s %-*s %-*s %-*s\n",
165 fw, "---------", fw, "---------",
166 fw, "-----", fw, "---------",
167 fw, "-----------");
168 for (i = 0; i < nfreqs; i++) {
169 out_printf("%-4d %-*s ", i, fw, printnum(freq[i]));
170 out_printf("%-*s ", fw, printnum(mag[i]));
171 out_printf("%-*s ", fw, printnum(phase[i]));
172 out_printf("%-*s ", fw, printnum(nmag[i]));
173 out_printf("%-*s\n", fw, printnum(nphase[i]));
174 }
175 out_send("\n");
176 }
177
178 vec_dlfree(dl0);
179 txfree((char*)freq);
180 txfree((char*)mag);
181 txfree((char*)phase);
182 txfree((char*)nmag);
183 txfree((char*)nphase);
184 if (polydegree) {
185 txfree((char*)grid);
186 txfree((char*)stuff);
187 }
188 }
189
190
191 /*
192 * CKTfour() - perform fourier analysis of an output vector.
193 * Due to the construction of the program which places all the
194 * output data in the post-processor, the fourier analysis can not
195 * be done directly. This function allows the post processor to
196 * hand back vectors of time and data values to have the fourier analysis
197 * performed on them.
198 *
199 */
200
201 /*ARGSUSED*/
202 static int
CKTfour(ndata,numFreq,thd,Time,Value,FundFreq,Freq,Mag,Phase,nMag,nPhase)203 CKTfour(ndata,numFreq,thd,Time,Value,FundFreq,Freq,Mag,Phase,nMag,nPhase)
204 int ndata; /* number of entries in the Time and Value arrays */
205 int numFreq; /* number of harmonics to calculate */
206 double *thd; /* total harmonic distortion (percent) to be returned */
207 double *Time; /* times at which the voltage/current values were measured*/
208 double *Value; /* voltage or current vector whose transform is desired */
209 double FundFreq;/* the fundamental frequency of the analysis */
210 double *Freq; /* the frequency value of the various harmonics */
211 double *Mag; /* the Magnitude of the fourier transform */
212 double *Phase; /* the Phase of the fourier transform */
213 double *nMag; /* the normalized magnitude of the transform: nMag(fund)=1*/
214 double *nPhase; /* the normalized phase of the transform: Nphase(fund)=0 */
215 /* note we can consider these as a set of arrays: The sizes are:
216 * Time[ndata], Value[ndata]
217 * Freq[numFreq],Mag[numfreq],Phase[numfreq],nMag[numfreq],nPhase[numfreq]
218 * The arrays must all be allocated by the caller.
219 * The Time and Value array must be reasonably distributed over at
220 * least one full period of the fundamental Frequency for the
221 * fourier transform to be useful. The function will take the
222 * last period of the frequency as data for the transform.
223 */
224
225 {
226 /* we are assuming that the caller has provided exactly one period
227 * of the fundamental frequency.
228 */
229 int i;
230 int j;
231 double tmp;
232 /* clear output/computation arrays */
233
234 for(i=0;i<numFreq;i++) {
235 Mag[i]=0;
236 Phase[i]=0;
237 }
238 for(i=0;i<ndata;i++) {
239 for(j=0;j<numFreq;j++) {
240 Mag[j] += (Value[i]*sin(j*2.0*M_PI*i/((double) ndata)));
241 Phase[j] += (Value[i]*cos(j*2.0*M_PI*i/((double) ndata)));
242 }
243 }
244
245 Mag[0] = Phase[0]/ndata;
246 Phase[0]=nMag[0]=nPhase[0]=Freq[0]=0;
247 *thd = 0;
248 for(i=1;i<numFreq;i++) {
249 tmp = Mag[i]*2.0 /ndata;
250 Phase[i] *= 2.0/ndata;
251 Freq[i] = i * FundFreq;
252 Mag[i] = sqrt(tmp*tmp+Phase[i]*Phase[i]);
253 Phase[i] = atan2(Phase[i],tmp)*180.0/M_PI;
254 nMag[i] = Mag[i]/Mag[1];
255 nPhase[i] = Phase[i]-Phase[1];
256 if(i>1) *thd += nMag[i]*nMag[i];
257 }
258 *thd = 100*sqrt(*thd);
259 return(OK);
260 }
261