diff --git a/src/frontend/Makefile.am b/src/frontend/Makefile.am index 07a77e90e..46a47f5f6 100644 --- a/src/frontend/Makefile.am +++ b/src/frontend/Makefile.am @@ -59,6 +59,8 @@ libfte_la_SOURCES = \ com_qpac.h \ com_qpnoise.c \ com_qpnoise.h \ + com_hbosc.c \ + com_hbosc.h \ com_hb.c \ com_hb.h \ com_checkpoint.c \ diff --git a/src/frontend/com_commands.h b/src/frontend/com_commands.h index 36db9e5f0..bb0f3cfb3 100644 --- a/src/frontend/com_commands.h +++ b/src/frontend/com_commands.h @@ -10,6 +10,8 @@ void com_optimize(wordlist *wl); /* Enhancement-130 */ void com_qpss(wordlist *wl); /* Enhancement-133 */ void com_qpac(wordlist *wl); /* Enhancement-137 */ void com_qpnoise(wordlist *wl); /* Enhancement-138 */ +void com_hbosc(wordlist *wl); /* Enhancement-140 */ +void com_phasenoise(wordlist *wl); /* Enhancement-140 */ void com_hb(wordlist *wl); /* Enhancement-134 */ void com_savestate(wordlist *wl); /* Enhancement-131 */ void com_loadstate(wordlist *wl); /* Enhancement-131 */ diff --git a/src/frontend/com_hbosc.c b/src/frontend/com_hbosc.c new file mode 100644 index 000000000..68d95add5 --- /dev/null +++ b/src/frontend/com_hbosc.c @@ -0,0 +1,173 @@ +/********** +Enhancement-140: autonomous harmonic balance for oscillators + phase noise. + + hbosc [fguess] [tstab] -- autonomous HB: find the oscillator's + steady state (harmonics + frequency) + phasenoise [points] -- the phase-noise spectrum L(df) + +`hbosc` runs a short transient (from the deck's .ic) to seed the limit cycle, estimates +the oscillation frequency and amplitude, and hands them to HBOSCanalyze() which refines +(V, w0) by a bordered Newton and retains the operating point. `phasenoise` then extracts +the perturbation projection vector (PPV) and folds the device noise to L(df). The engines +are in spicelib/analysis/dcpss.c. +**********/ + +#include "ngspice/ngspice.h" +#include "ngspice/cpdefs.h" +#include "ngspice/cktdefs.h" +#include "ngspice/ftedefs.h" +#include "ngspice/dvec.h" +#include "ngspice/fteext.h" +#include "ngspice/wordlist.h" +#include "ngspice/cpextern.h" + +#include "circuits.h" +#include "com_hbosc.h" + +static void hbosc_run_cmd(const char *cmdstr) +{ + wordlist *wl = cp_lexer((char *) cmdstr); + int i; + if (!wl || !wl->wl_word) { if (wl) wl_free(wl); return; } + for (i = 0; cp_coms[i].co_comname; i++) + if (strcasecmp(cp_coms[i].co_comname, wl->wl_word) == 0) + break; + if (cp_coms[i].co_comname && cp_coms[i].co_func) + cp_coms[i].co_func(wl->wl_next); + wl_free(wl); +} + +static double hboscnum(const char *w) +{ + char *s = (char *) w; + double v = 0.0; + if (ft_numparse(&s, FALSE, &v) < 0) + v = atof(w); + return v; +} + +static int hbosc_node(CKTcircuit *ckt, const char *name) +{ + int numNames = 0, i, num = 0; + IFuid *nameList = NULL; + if (CKTnames(ckt, &numNames, &nameList) != OK || !nameList) + return 0; + for (i = 0; i < numNames; i++) + if (nameList[i] && strcmp((const char *) nameList[i], name) == 0) { num = i + 1; break; } + tfree(nameList); + return num; +} + +void +com_hbosc(wordlist *wl) +{ + CKTcircuit *ckt; + const char *oscname; + double fguess = 0.0, tstab = 0.0, ampseed, f0est; + int oscNode, K, verbose, err; + char cmd[128]; + struct pnode *pn; + struct dvec *v, *sc; + double *tt, *vv; + int n, i0, i, ncross; + double vmax, tlast, tprev; + + if (!ft_curckt || !ft_curckt->ci_ckt) { + fprintf(cp_err, "Error: hbosc: there is no circuit loaded.\n"); + return; + } + ckt = ft_curckt->ci_ckt; + if (!wl || !wl->wl_next) { + fprintf(cp_err, "Usage: hbosc [fguess] [tstab] (the deck needs a .ic to start the oscillation)\n"); + return; + } + oscname = wl->wl_word; + K = (int) hboscnum(wl->wl_next->wl_word); + if (wl->wl_next->wl_next) { + fguess = hboscnum(wl->wl_next->wl_next->wl_word); + if (wl->wl_next->wl_next->wl_next) + tstab = hboscnum(wl->wl_next->wl_next->wl_next->wl_word); + } + if (K < 1) { fprintf(cp_err, "Error: hbosc: need K >= 1.\n"); return; } + if (fguess <= 0.0) fguess = 1e6; /* a default if none given */ + if (tstab <= 0.0) tstab = 300.0 / fguess; /* ~300 periods to settle */ + + /* build up the limit cycle with a transient (uic honours the deck's .ic) */ + (void) snprintf(cmd, sizeof cmd, "tran %.6g %.6g uic", 1.0/(fguess*40.0), tstab); + hbosc_run_cmd(cmd); + + /* fetch v(oscnode) + its time scale */ + (void) snprintf(cmd, sizeof cmd, "v(%s)", oscname); + pn = ft_getpnames_from_string(cmd, TRUE); + if (!pn) { fprintf(cp_err, "Error: hbosc: cannot read v(%s).\n", oscname); return; } + v = ft_evaluate(pn); + sc = (v && v->v_scale) ? v->v_scale : vec_get("time"); + if (!v || !isreal(v) || v->v_length < 8 || !sc || sc->v_length < v->v_length) { + fprintf(cp_err, "Error: hbosc: no usable transient for v(%s).\n", oscname); + if (pn && !pn->pn_value && v) vec_free(v); + if (pn) free_pnode(pn); + return; + } + tt = sc->v_realdata; vv = v->v_realdata; n = v->v_length; + + /* estimate the oscillation: amplitude = max|v| over the last third; frequency from + * the mean spacing of upward zero crossings there. */ + i0 = (2*n)/3; + vmax = 0.0; + for (i = i0; i < n; i++) if (fabs(vv[i]) > vmax) vmax = fabs(vv[i]); + ncross = 0; tlast = tprev = 0.0; + for (i = i0 + 1; i < n; i++) + if (vv[i-1] <= 0.0 && vv[i] > 0.0) { /* upward crossing */ + double tc = tt[i-1] + (tt[i]-tt[i-1]) * (-vv[i-1])/(vv[i]-vv[i-1]); + if (ncross == 0) tprev = tc; + tlast = tc; ncross++; + } + f0est = (ncross > 1 && tlast > tprev) ? (double)(ncross-1)/(tlast-tprev) : fguess; + + ampseed = vmax; + oscNode = hbosc_node(ckt, oscname); + if (pn && !pn->pn_value && v) vec_free(v); + if (pn) free_pnode(pn); + if (oscNode <= 0) { fprintf(cp_err, "Error: hbosc: unknown oscillator node '%s'.\n", oscname); return; } + if (vmax < 1e-9) { + fprintf(cp_err, "Error: hbosc: no oscillation detected (add a `.ic` to start it).\n"); + return; + } + + verbose = cp_getvar("hbosc_verbose", CP_BOOL, NULL, 0); + err = HBOSCanalyze(ckt, oscNode, K, 0, f0est, ampseed, 60, 1e-11, verbose ? 1 : 0); + if (err != OK) + fprintf(cp_err, "hbosc: autonomous harmonic balance did not complete (error %d).\n", err); +} + +void +com_phasenoise(wordlist *wl) +{ + CKTcircuit *ckt; + double fstart, fstop; + int npts = 21, verbose, err; + + if (!ft_curckt || !ft_curckt->ci_ckt) { + fprintf(cp_err, "Error: phasenoise: there is no circuit loaded.\n"); + return; + } + ckt = ft_curckt->ci_ckt; + if (!wl || !wl->wl_next) { + fprintf(cp_err, "Usage: phasenoise [points] (run `hbosc` first)\n"); + return; + } + fstart = hboscnum(wl->wl_word); + fstop = hboscnum(wl->wl_next->wl_word); + if (wl->wl_next->wl_next) + npts = (int) hboscnum(wl->wl_next->wl_next->wl_word); + if (fstart <= 0.0 || fstop < fstart) { + fprintf(cp_err, "Error: phasenoise: need 0 < fstart <= fstop.\n"); + return; + } + if (npts < 1) npts = 1; + + verbose = cp_getvar("phasenoise_verbose", CP_BOOL, NULL, 0); + err = PhaseNoiseAnalyze(ckt, fstart, fstop, npts, verbose ? 1 : 0); + if (err != OK) + fprintf(cp_err, "phasenoise: did not complete (error %d).\n", err); +} diff --git a/src/frontend/com_hbosc.h b/src/frontend/com_hbosc.h new file mode 100644 index 000000000..d54558f38 --- /dev/null +++ b/src/frontend/com_hbosc.h @@ -0,0 +1,8 @@ +#ifndef ngspice_COM_HBOSC_H +#define ngspice_COM_HBOSC_H + +/* Enhancement-140: autonomous HB for oscillators + phase noise. */ +void com_hbosc(wordlist *wl); +void com_phasenoise(wordlist *wl); + +#endif diff --git a/src/frontend/commands.c b/src/frontend/commands.c index 46fb9e197..565fa7892 100644 --- a/src/frontend/commands.c +++ b/src/frontend/commands.c @@ -444,6 +444,14 @@ struct comm spcp_coms[] = { { 040, 040, 040, 040 }, E_DEFHMASK, 2, LOTS, NULL, "output_node f_in [cyclo] : two-tone QPnoise -- output/input noise density at f_in, folding device noise over all sidebands around the `qpss ... hb` operating point (`cyclo` = cyclostationary PSD)." }, + { "hbosc", com_hbosc, TRUE, FALSE, /* Enhancement-140 */ + { 040, 040, 040, 040 }, E_DEFHMASK, 2, LOTS, + NULL, + "oscnode K [fguess] [tstab] : autonomous harmonic balance -- an oscillator's steady state (harmonics + oscillation frequency), seeded from a transient." }, + { "phasenoise", com_phasenoise, TRUE, FALSE, /* Enhancement-140 */ + { 040, 040, 040, 040 }, E_DEFHMASK, 2, LOTS, + NULL, + "fstart fstop [points] : oscillator phase-noise spectrum L(df) via the PPV, around the `hbosc` operating point." }, { "hb", com_hb, TRUE, FALSE, /* Enhancement-134 */ { 040, 040, 040, 040 }, E_DEFHMASK, 2, LOTS, NULL, diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index a0249e3d2..e865fcfec 100644 --- a/src/include/ngspice/cktdefs.h +++ b/src/include/ngspice/cktdefs.h @@ -479,6 +479,8 @@ extern int HBanalyze(CKTcircuit *, double f0, int K, int P, int maxiter, double extern int QPSShb(CKTcircuit *, double f1, double f2, int K1, int K2, int P1, int P2, int maxiter, double tol, int verbose); /* E-136 */ extern int QPACanalyze(CKTcircuit *, double f_in, int verbose); /* E-137 */ extern int QPnoiseAnalyze(CKTcircuit *, int outNode, double f_in, int cyclo, int verbose); /* E-138 / -139 */ +extern int HBOSCanalyze(CKTcircuit *, int oscNode, int K, int P, double f0seed, double ampseed, int maxiter, double tol, int verbose); /* E-140 */ +extern int PhaseNoiseAnalyze(CKTcircuit *, double fstart, double fstop, int npts, int verbose); /* E-140 */ #endif diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c index 0df416eab..953177eab 100644 --- a/src/spicelib/analysis/dcpss.c +++ b/src/spicelib/analysis/dcpss.c @@ -2425,6 +2425,276 @@ QPnoiseAnalyze(CKTcircuit *ckt, int outNode, double f_in, int cyclo, int verbose } +/* ====================================================================== + * Enhancement-140: oscillator phase noise. + * + * (1) HBOSCanalyze -- AUTONOMOUS harmonic balance. An oscillator has no + * driving source, so the HB residual F(V) = I_R + [dq/dt] = 0 is solved + * for the harmonics V *and* the unknown oscillation frequency w0, with a + * phase gauge (the solution's phase is free). The Jacobian dF/dV is the + * conversion matrix H -- SINGULAR, its right null space the phase mode + * u_k = jk V_k (shifting the oscillator's phase is a symmetry). Newton on + * the bordered system [ H dF/dw0 ; u*^T 0 ] [dV; dw0] = [-F; 0] + * (nonsingular by bordering) refines (V, w0) from a transient seed. + * (2) PhaseNoiseAnalyze -- the phase-noise spectrum L(df). The PPV (Demir's + * perturbation projection vector) is the LEFT null vector of H, normalized + * so = 1. The phase diffusion constant c is the device + * noise folded through the PPV (same machinery as pnoise, transfer -> PPV), + * and L(df) = 10 log10( f0^2 c / df^2 ) -- the classic 1/df^2 (-20 dB/dec) + * oscillator phase noise, saturating into a Lorentzian near the carrier. + * ====================================================================== */ + +static struct pac_harm osc_hd; /* retained oscillator conversion data */ +static double *osc_Vr = NULL, *osc_Vi = NULL; /* [Ntot] retained oscillator spectrum */ +static double osc_f0 = 0.0; +static int osc_node = 0, osc_valid = 0; + +int +HBOSCanalyze(CKTcircuit *ckt, int oscNode, int K, int Pin, double f0seed, + double ampseed, int maxiter, double tol, int verbose) +{ + int N = SMPmatSize(ckt->CKTmatrix); + int P = Pin > 0 ? Pin : ((8*K < 32) ? 32 : 8*K); + int Ntot = (2*K+1)*N, Naug = Ntot + 1; + int iter, s, i, k, c, rc = E_ITERLIM, have_hd = 0; + double f0 = f0seed, w0 = 2.0*M_PI*f0, fnorm = 0.0; + double *Vr, *Vi, *vsamp, *IRr, *IRi, *Kr, *Ki, *Jr, *Ji, *Fr, *Fi, *Ar, *Ai, *br, *bi; + struct pac_harm hd; + + if (N <= 0 || K < 1 || oscNode <= 0 || oscNode > N) { + fprintf(stderr, "hbosc: bad size or oscillator node.\n"); return E_PARMVAL; + } + if (Naug > 1600) { fprintf(stderr, "hbosc: system too large (reduce K).\n"); return E_PARMVAL; } + + Vr = TMALLOC(double, Ntot); Vi = TMALLOC(double, Ntot); + IRr = TMALLOC(double, Ntot); IRi = TMALLOC(double, Ntot); + Kr = TMALLOC(double, Ntot); Ki = TMALLOC(double, Ntot); + Fr = TMALLOC(double, Ntot); Fi = TMALLOC(double, Ntot); + Jr = TMALLOC(double, (size_t)Ntot*(size_t)Ntot); Ji = TMALLOC(double, (size_t)Ntot*(size_t)Ntot); + Ar = TMALLOC(double, (size_t)Naug*(size_t)Naug); Ai = TMALLOC(double, (size_t)Naug*(size_t)Naug); + br = TMALLOC(double, Naug); bi = TMALLOC(double, Naug); + vsamp = TMALLOC(double, (size_t)N*(size_t)P); + + /* seed: fundamental at the osc node (amplitude A -> |V_1| = A/2, real, phase 0) */ + for (i = 0; i < Ntot; i++) { Vr[i] = 0.0; Vi[i] = 0.0; } + Vr[(size_t)(1+K)*(size_t)N + (size_t)(oscNode-1)] = ampseed * 0.5; + Vr[(size_t)(K-1)*(size_t)N + (size_t)(oscNode-1)] = ampseed * 0.5; + + for (iter = 0; iter < maxiter; iter++) { + w0 = 2.0*M_PI*f0; + for (s = 0; s < P; s++) + for (i = 0; i < N; i++) { + double v = 0.0; + for (k = -K; k <= K; k++) { + double ang = 2.0*M_PI*k*s/(double)P; + v += Vr[(size_t)(k+K)*(size_t)N+(size_t)i]*cos(ang) + - Vi[(size_t)(k+K)*(size_t)N+(size_t)i]*sin(ang); + } + vsamp[(size_t)s*(size_t)N+(size_t)i] = v; + } + if (hb_extract(ckt, vsamp, N, P, K, &hd, IRr, IRi)) { + fprintf(stderr, "hbosc: device extraction failed.\n"); rc = E_PARMVAL; break; + } + have_hd = 1; + pac_build_matrix(&hd, f0, 0.0, Jr, Ji); + { /* I_C = (J - Jg) V */ + struct pac_harm hg = hd; + double *Jgr = TMALLOC(double, (size_t)Ntot*(size_t)Ntot); + double *Jgi = TMALLOC(double, (size_t)Ntot*(size_t)Ntot); + hg.Cmr = TMALLOC(double, (size_t)hd.nnz*(size_t)(hd.H+1)); + hg.Cmi = TMALLOC(double, (size_t)hd.nnz*(size_t)(hd.H+1)); + pac_build_matrix(&hg, f0, 0.0, Jgr, Jgi); + FREE(hg.Cmr); FREE(hg.Cmi); + for (i = 0; i < Ntot; i++) { + double cr = 0, ci = 0; + for (c = 0; c < Ntot; c++) { + double ar = Jr[(size_t)i*(size_t)Ntot+(size_t)c] - Jgr[(size_t)i*(size_t)Ntot+(size_t)c]; + double ai = Ji[(size_t)i*(size_t)Ntot+(size_t)c] - Jgi[(size_t)i*(size_t)Ntot+(size_t)c]; + cr += ar*Vr[c] - ai*Vi[c]; ci += ar*Vi[c] + ai*Vr[c]; + } + Kr[i] = cr; Ki[i] = ci; + } + FREE(Jgr); FREE(Jgi); + } + fnorm = 0.0; + for (i = 0; i < Ntot; i++) { + Fr[i] = IRr[i] + Kr[i]; Fi[i] = IRi[i] + Ki[i]; + fnorm += Fr[i]*Fr[i] + Fi[i]*Fi[i]; + } + fnorm = sqrt(fnorm); + if (verbose) + fprintf(stderr, "hbosc iter %2d: |F| = %.6e f0 = %.8g Hz\n", iter, fnorm, f0); + if (isnan(fnorm) || fnorm > 1e300) { pac_free_harmonics(&hd); have_hd = 0; break; } + if (fnorm < tol) { rc = OK; break; } /* converged: keep hd for retention */ + + /* bordered Newton: [ H d ; u*^T 0 ] [dV; dw0] = [-F; 0] */ + memset(Ar, 0, (size_t)Naug*(size_t)Naug*sizeof(double)); + memset(Ai, 0, (size_t)Naug*(size_t)Naug*sizeof(double)); + for (i = 0; i < Ntot; i++) + for (c = 0; c < Ntot; c++) { + Ar[(size_t)i*(size_t)Naug+(size_t)c] = Jr[(size_t)i*(size_t)Ntot+(size_t)c]; + Ai[(size_t)i*(size_t)Naug+(size_t)c] = Ji[(size_t)i*(size_t)Ntot+(size_t)c]; + } + for (i = 0; i < Ntot; i++) { /* d = dF/dw0 = I_C / w0 */ + Ar[(size_t)i*(size_t)Naug+(size_t)Ntot] = Kr[i]/w0; + Ai[(size_t)i*(size_t)Naug+(size_t)Ntot] = Ki[i]/w0; + } + for (i = 0; i < Ntot; i++) { /* row = conj(u), u_k = jk V_k */ + int kk = (i/N) - K; + Ar[(size_t)Ntot*(size_t)Naug+(size_t)i] = -(double)kk * Vi[i]; /* Re conj(u) */ + Ai[(size_t)Ntot*(size_t)Naug+(size_t)i] = -(double)kk * Vr[i]; /* Im conj(u) */ + } + for (i = 0; i < Ntot; i++) { br[i] = -Fr[i]; bi[i] = -Fi[i]; } + br[Ntot] = 0.0; bi[Ntot] = 0.0; + pac_free_harmonics(&hd); have_hd = 0; + if (pss_csolve(Naug, Ar, Ai, br, bi)) { + fprintf(stderr, "hbosc: singular augmented system.\n"); rc = E_SINGULAR; break; + } + for (i = 0; i < Ntot; i++) { Vr[i] += br[i]; Vi[i] += bi[i]; } + f0 += br[Ntot] / (2.0*M_PI); /* dw0 = Re(x[Ntot]) */ + } + + if (rc == OK && have_hd) { + int numNames, error; + IFuid *nameList = NULL; + /* retain the operating point for phasenoise */ + if (osc_valid) { pac_free_harmonics(&osc_hd); FREE(osc_Vr); FREE(osc_Vi); } + osc_hd = hd; /* transfer ownership of hd's arrays */ + osc_Vr = TMALLOC(double, Ntot); osc_Vi = TMALLOC(double, Ntot); + memcpy(osc_Vr, Vr, (size_t)Ntot*sizeof(double)); + memcpy(osc_Vi, Vi, (size_t)Ntot*sizeof(double)); + osc_f0 = f0; osc_node = oscNode; osc_valid = 1; + + error = CKTnames(ckt, &numNames, &nameList); + fprintf(stdout, "\nHBOSC: autonomous oscillator steady state\n" + " oscillation frequency f0 = %.9g Hz (converged, |F| = %.3e)\n" + " node harmonic frequency [Hz] |V| phase [deg]\n", + f0, fnorm); + for (i = 0; i < N; i++) { + const char *nm = (!error && i < numNames) ? (const char *) nameList[i] : "?"; + for (k = 0; k <= K; k++) { + double sc = (k == 0) ? 1.0 : 2.0; + double vr = sc*Vr[(size_t)(k+K)*(size_t)N+(size_t)i]; + double vi = sc*Vi[(size_t)(k+K)*(size_t)N+(size_t)i]; + fprintf(stdout, " %-8s %6d %16.6e %14.6e %10.3f\n", + nm, k, k*f0, hypot(vr, vi), (k==0)?0.0:atan2(vi, vr)*180.0/M_PI); + } + } + if (nameList) tfree(nameList); + } else { + if (have_hd) pac_free_harmonics(&hd); + if (rc != OK) + fprintf(stderr, "hbosc: did not converge to an oscillation (try a better fguess/tstab).\n"); + } + + FREE(Vr); FREE(Vi); FREE(IRr); FREE(IRi); FREE(Kr); FREE(Ki); + FREE(Fr); FREE(Fi); FREE(Jr); FREE(Ji); FREE(Ar); FREE(Ai); FREE(br); FREE(bi); FREE(vsamp); + return rc; +} + +int +PhaseNoiseAnalyze(CKTcircuit *ckt, double fstart, double fstop, int npts, int verbose) +{ + struct pac_harm *hd = &osc_hd; + int N, M, Ntot, i, j, k, rc = OK, onode, ni, mi, ei; + double f0 = osc_f0; + double *Psr, *Psi, Pcar, freq, mult; + NOISEAN nj; Ndata data; JOB *oldJob; + + if (!osc_valid) { + fprintf(stderr, "phasenoise: no oscillator operating point -- run `hbosc` first.\n"); + return E_NOTFOUND; + } + N = hd->N; M = hd->M; Ntot = hd->Ntot; onode = osc_node; + Psr = TMALLOC(double, Ntot); Psi = TMALLOC(double, Ntot); + + /* carrier power at the oscillator node: single-sided amplitude A = 2|V_1|, mean + * square A^2/2 = 2|V_1|^2 (the reference the sideband noise is measured against). */ + { + size_t c1 = (size_t)(M+1)*(size_t)N + (size_t)(onode-1); + Pcar = 2.0 * (osc_Vr[c1]*osc_Vr[c1] + osc_Vi[c1]*osc_Vi[c1]); + if (Pcar <= 0.0) Pcar = 1.0; + } + + /* bias the devices at the oscillator op-point (phase-0 sample) for their noise PSD */ + for (j = 1; j <= N; j++) { + double v = 0.0; + for (k = -M; k <= M; k++) + v += osc_Vr[(size_t)(k+M)*(size_t)N+(size_t)(j-1)]; + ckt->CKTrhsOld[j] = v; + } + ckt->CKTrhsOld[0] = 0.0; + ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG; + CKTload(ckt); + memset(&nj, 0, sizeof(nj)); + nj.NstartFreq = fstart; nj.NstopFreq = fstop; nj.NnumSteps = npts; + nj.NstpType = 0; nj.NStpsSm = 0; nj.JOBname = "phasenoise"; + memset(&data, 0, sizeof(data)); data.prtSummary = FALSE; + oldJob = ckt->CKTcurJob; ckt->CKTcurJob = (JOB *) &nj; + for (i = 0; i < DEVmaxnum; i++) + if (DEVices[i] && DEVices[i]->DEVnoise && ckt->CKThead[i]) { + double dummy = 0.0; + DEVices[i]->DEVnoise(N_DENS, N_OPEN, ckt->CKThead[i], ckt, &data, &dummy); + } + + fprintf(stdout, "\nPhaseNoise: oscillator phase noise (f0 = %.9g Hz, carrier power %.4e)\n" + " offset [Hz] L(df) [dBc/Hz]\n", f0, Pcar); + mult = (npts > 1) ? pow(fstop/fstart, 1.0/(double)(npts-1)) : 1.0; + + for (freq = fstart, i = 0; i < npts; i++, freq *= mult) { + /* adjoint of H at OFFSET f_in = df, unit at the CARRIER sideband (m=1) of the + * osc node: Psi is the transimpedance from a noise injection at every (node, + * sideband) to the carrier. As df -> 0, H(df) -> the singular limit-cycle matrix, + * so Psi (through the phase mode) blows up as 1/df -> Sv ~ 1/df^2 = phase noise. */ + double *Ar = TMALLOC(double, (size_t)Ntot*(size_t)Ntot); + double *Ai = TMALLOC(double, (size_t)Ntot*(size_t)Ntot); + double Sv = 0.0; + memset(Ar, 0, (size_t)Ntot*(size_t)Ntot*sizeof(double)); + memset(Ai, 0, (size_t)Ntot*(size_t)Ntot*sizeof(double)); + for (ni = 0; ni <= 2*M; ni++) + for (mi = 0; mi <= 2*M; mi++) { + int dm = (ni - M) - (mi - M); + double omega = 2.0*M_PI*(freq + (double)(mi - M)*f0); + for (ei = 0; ei < hd->nnz; ei++) { + size_t hi = (size_t)ei*(size_t)(hd->H+1) + (size_t)abs(dm); + double gr = hd->Gmr[hi], gi = hd->Gmi[hi], cr = hd->Cmr[hi], ci = hd->Cmi[hi]; + double er, eii; size_t row, col; + if (dm < 0) { gi = -gi; ci = -ci; } + er = gr - omega*ci; eii = gi + omega*cr; + row = (size_t)ni*(size_t)N + (size_t)(hd->rr[ei]-1); + col = (size_t)mi*(size_t)N + (size_t)(hd->cc[ei]-1); + Ar[col*(size_t)Ntot+row] += er; /* transpose: [col][row] = H^T */ + Ai[col*(size_t)Ntot+row] += eii; + } + } + memset(Psr, 0, (size_t)Ntot*sizeof(double)); + memset(Psi, 0, (size_t)Ntot*sizeof(double)); + Psr[(size_t)(M+1)*(size_t)N + (size_t)(onode-1)] = 1.0; /* carrier sideband m=1 */ + if (pss_csolve(Ntot, Ar, Ai, Psr, Psi) == 0) { + data.freq = freq; data.delFreq = 0.0; data.prtSummary = FALSE; + for (k = -M; k <= M; k++) { + double dens = 0.0; + size_t blk = (size_t)(k+M)*(size_t)N; + for (j = 1; j <= N; j++) { ckt->CKTrhs[j] = Psr[blk+(size_t)(j-1)]; ckt->CKTirhs[j] = Psi[blk+(size_t)(j-1)]; } + ckt->CKTrhs[0] = 0.0; ckt->CKTirhs[0] = 0.0; + for (ei = 0; ei < DEVmaxnum; ei++) + if (DEVices[ei] && DEVices[ei]->DEVnoise && ckt->CKThead[ei]) + DEVices[ei]->DEVnoise(N_DENS, N_CALC, ckt->CKThead[ei], ckt, &data, &dens); + Sv += dens; + } + } + FREE(Ar); FREE(Ai); + fprintf(stdout, " %14.6e %12.4f\n", freq, 10.0*log10(Sv / Pcar)); + if (npts == 1) break; + } + ckt->CKTcurJob = oldJob; + (void) verbose; + + FREE(Psr); FREE(Psi); + return rc; +} + + int DCpss(CKTcircuit *ckt, int restart) /* forced restart flag */ diff --git a/visualc/vngspice.vcxproj b/visualc/vngspice.vcxproj index 36d4afae3..31a9f2402 100644 --- a/visualc/vngspice.vcxproj +++ b/visualc/vngspice.vcxproj @@ -887,6 +887,7 @@ + @@ -1507,6 +1508,7 @@ +