From ab0cc782c41a88fc9fa681943192552fa3f612f1 Mon Sep 17 00:00:00 2001 From: Meisam Bahadori Date: Sun, 26 Jul 2026 11:56:25 +0200 Subject: [PATCH] pa-132 --- src/include/ngspice/pssdefs.h | 9 ++ src/spicelib/analysis/dcpss.c | 243 ++++++++++++++++++++++++++++++-- src/spicelib/analysis/psssetp.c | 8 +- src/spicelib/parser/inp2dot.c | 87 ++++++++++++ 4 files changed, 333 insertions(+), 14 deletions(-) diff --git a/src/include/ngspice/pssdefs.h b/src/include/ngspice/pssdefs.h index 3fe2fe83f..92c44f3e7 100644 --- a/src/include/ngspice/pssdefs.h +++ b/src/include/ngspice/pssdefs.h @@ -68,6 +68,14 @@ typedef struct { * sweep fields + PACmaxSideband. */ int PSSdoPXF; /* 1 if this job runs a PXF sweep */ CKTnode *PxOutNode; /* PXF output node (reference = ground) */ + + /* Enhancement-132: periodic S-parameters (.psp). Runs PSS, then for each RF + * port injects a unit current in sideband 0 through the conversion matrix and + * reads the per-sideband port waves, building the periodic scattering matrix + * S^(k) = B^(k)*A^-1 vs input frequency. Reuses the PAC sweep fields + * (PACfStart/Stop/points/stepType/maxSideband). For a time-invariant circuit + * the sideband-0 result reduces to the ordinary .sp S-matrix. */ + int PSSdoPSP; /* 1 if this job runs a periodic S-parameter sweep */ } PSSan; enum { @@ -91,6 +99,7 @@ enum { PXF_DO, PXF_OUT, PNOISE_CYCLO, + PSP_DOPSP, }; #endif diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c index a1576e41e..8aa98fd4d 100644 --- a/src/spicelib/analysis/dcpss.c +++ b/src/spicelib/analysis/dcpss.c @@ -15,6 +15,11 @@ #include "ngspice/noisedef.h" /* Enhancement-124: NOISEAN/Ndata + CKTnoise for pnoise */ #include "ngspice/sperror.h" #include "ngspice/fteext.h" +#ifdef RFSPICE +#include "vsrc/vsrcdefs.h" /* Enhancement-132: RF port fields (z0, ki, branch) */ +#include "../maths/dense/dense.h" /* Enhancement-132: complex S = B*A^-1 */ +#include "../maths/dense/denseinlines.h" +#endif #ifdef XSPICE /* gtri - add - wbk - Add headers */ @@ -400,22 +405,15 @@ pac_free_harmonics(struct pac_harm *hd) FREE(hd->B0r); FREE(hd->B0i); } -/* Assemble the conversion matrix H_{nm} = G_{n-m} + j*omega_m*C_{n-m} at input - * frequency f_in and solve for a stimulus injected in the 0-th sideband. When - * `use_src` is set and the netlist supplied an `AC` source, the captured source - * RHS B_0 is the stimulus; otherwise a unit current is injected at node `inode`. - * The solution X (all sidebands, length Ntot) is written to Xr/Xi, which the caller - * allocates. Returns 0 on success, 1 if the matrix is singular. */ -static int -pac_solve_at(struct pac_harm *hd, double f0, double f_in, int inode, int use_src, - double *Xr, double *Xi) +/* Enhancement-121: assemble the dense (2M+1)N complex conversion matrix H_{nm} = + * G_{n-m} + j*omega_m*C_{n-m} at input frequency f_in into (Ar,Ai). Shared by the + * PAC (pac_solve_at) and PSP (psp_solve_port) solves. */ +static void +pac_build_matrix(struct pac_harm *hd, double f0, double f_in, double *Ar, double *Ai) { int N = hd->N, M = hd->M, H = hd->H, nnz = hd->nnz, Ntot = hd->Ntot; - int ni, mi, n, mm, e, rc, j; - double *Ar, *Ai; + int ni, mi, n, mm, e; - Ar = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); - Ai = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); 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++) { @@ -443,6 +441,23 @@ pac_solve_at(struct pac_harm *hd, double f0, double f_in, int inode, int use_src } } } +} + +/* Solve the conversion matrix at input frequency f_in for a stimulus injected in the + * 0-th sideband. When `use_src` is set and the netlist supplied an `AC` source, the + * captured source RHS B_0 is the stimulus; otherwise a unit current is injected at + * node `inode`. The solution X (all sidebands, length Ntot) is written to Xr/Xi, + * which the caller allocates. Returns 0 on success, 1 if the matrix is singular. */ +static int +pac_solve_at(struct pac_harm *hd, double f0, double f_in, int inode, int use_src, + double *Xr, double *Xi) +{ + int N = hd->N, M = hd->M, Ntot = hd->Ntot, rc, j; + double *Ar, *Ai; + + Ar = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + Ai = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + pac_build_matrix(hd, f0, f_in, Ar, Ai); /* stimulus in the 0-th sideband: netlist AC source RHS, or a unit current */ memset(Xr, 0, (size_t)Ntot * sizeof(double)); @@ -460,6 +475,33 @@ pac_solve_at(struct pac_harm *hd, double f0, double f_in, int inode, int use_src return rc; } +#ifdef RFSPICE +/* Enhancement-132: solve the conversion matrix with one RF port DRIVEN (V=1) in the + * 0-th sideband -- the per-port excitation for periodic S-parameters. The port is a + * unit voltage source: set its branch-equation RHS to 1 (exactly as .sp's + * VSRCspupdate does), so V(pos)-V(neg)=1 at that port and 0 at the others (which + * stay z0-terminated by the g0 shunt already in the Jacobian). Same conversion + * matrix as pac_solve_at, different RHS. */ +static int +psp_solve_port(struct pac_harm *hd, double f0, double f_in, int branch, + double *Xr, double *Xi) +{ + int N = hd->N, M = hd->M, Ntot = hd->Ntot, rc; + double *Ar, *Ai; + + Ar = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + Ai = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + pac_build_matrix(hd, f0, f_in, Ar, Ai); + + memset(Xr, 0, (size_t)Ntot * sizeof(double)); + memset(Xi, 0, (size_t)Ntot * sizeof(double)); + Xr[(size_t)M * (size_t)N + (size_t)(branch - 1)] = 1.0; + rc = pss_csolve(Ntot, Ar, Ai, Xr, Xi); + FREE(Ar); FREE(Ai); + return rc; +} +#endif + /* pick sideband count M from the requested harmonics; conversion matrix is * (2M+1)*N, so cap M so the dense solve stays small. Returns 0 if unusable. */ static int @@ -1017,6 +1059,174 @@ pxf_sweep(CKTcircuit *ckt, PSSan *job) } +#ifdef RFSPICE +/* Enhancement-132: periodic S-parameters (.psp). After PSS, excite each RF port in + * turn with a unit current in the 0-th sideband through the conversion matrix, read + * the per-sideband port waves a/b (Kurosawa power waves, matching .sp's convention), + * and form the periodic scattering matrix S^(k) = B^(k) * A^-1 at each swept input + * frequency. S = B*A^-1 is invariant to the excitation basis, so unit-current + * injection yields the same S as .sp's power-wave excitation; for a time-invariant + * circuit the conversion matrix is block-diagonal and the sideband-0 block reduces + * exactly to the ordinary .sp S-matrix. */ +static void +psp_sweep(CKTcircuit *ckt, PSSan *job) +{ + int N = job->PSSopMsize, np = ckt->CKTportCount; + int M, Ksb, nsb, stepType = job->PACstepType, npts = job->PACpoints; + int i, j, s, p, numNames, error; + double f0 = job->PSSopFreq, fstart = job->PACfStart, fstop = job->PACfStop; + double freq, mult, linstep; + struct pac_harm hd; + double *Xr, *Xi; + IFuid freqUid, *outNames = NULL; + runDesc *plot = NULL; + CMat *Amat = NULL, *Ainv = NULL, **Bmat = NULL; + + if (np < 1) { + fprintf(stderr, "PSP: no RF ports (add `portnum`/`z0` to voltage sources) " + "-- sweep skipped\n"); + return; + } + if (f0 <= 0.0 || fstart <= 0.0 || fstop < fstart || npts < 1) + return; + + /* every port's nodes + branch must live inside the conversion-matrix unknowns */ + for (p = 0; p < np; p++) { + VSRCinstance *pr = (VSRCinstance *) ckt->CKTrfPorts[p]; + if (pr->VSRCbranch > N || pr->VSRCposNode > N || pr->VSRCnegNode > N) { + fprintf(stderr, "PSP: port %d outside the PSS unknown set -- skipped\n", + pr->VSRCportNum); + return; + } + } + + M = pac_choose_M(ckt, job); + if (M < 1) { fprintf(stderr, "PSP: conversion matrix too large -- skipped\n"); return; } + if (pac_extract_harmonics(ckt, job, M, &hd)) + return; + + Ksb = job->PACmaxSideband; + if (Ksb < 0) Ksb = 0; + if (Ksb > M) Ksb = M; + nsb = 2 * Ksb + 1; + + /* output vectors: S__ per sideband (dest = measured port, src = + * excited). Sideband 0 uses the plain name so it lines up with an ordinary + * .sp run; conversion sidebands get an _usb/_lsb suffix. */ + numNames = nsb * np * np; + outNames = TMALLOC(IFuid, numNames); + { + int idx = 0; + for (s = 0; s < nsb; s++) { + int k = s - Ksb; + for (j = 1; j <= np; j++) + for (i = 1; i <= np; i++) { + char nm[64]; + if (k == 0) + snprintf(nm, sizeof(nm), "S_%d_%d", j, i); + else + snprintf(nm, sizeof(nm), "S_%d_%d_%csb%d", + j, i, (k > 0) ? 'u' : 'l', abs(k)); + SPfrontEnd->IFnewUid(ckt, &outNames[idx++], NULL, nm, UID_OTHER, NULL); + } + } + } + SPfrontEnd->IFnewUid(ckt, &freqUid, NULL, "frequency", UID_OTHER, NULL); + error = SPfrontEnd->OUTpBeginPlot(ckt, ckt->CKTcurJob, "PSP Analysis", + freqUid, IF_REAL, numNames, outNames, + IF_COMPLEX, &plot); + tfree(outNames); + if (error) { pac_free_harmonics(&hd); return; } + if (stepType != 0) + SPfrontEnd->OUTattributes(plot, NULL, OUT_SCALE_LOG, NULL); + + Xr = TMALLOC(double, hd.Ntot); + Xi = TMALLOC(double, hd.Ntot); + Amat = newcmat(np, np, 0.0, 0.0); + Bmat = TMALLOC(CMat *, nsb); + for (s = 0; s < nsb; s++) + Bmat[s] = newcmat(np, np, 0.0, 0.0); + + mult = (stepType == 1) ? pow(10.0, 1.0 / npts) : + (stepType == 2) ? pow(2.0, 1.0 / npts) : 0.0; + linstep = (npts > 1) ? (fstop - fstart) / (npts - 1) : 0.0; + + fprintf(stderr, "PSP sweep: %s from %.6g to %.6g Hz around f0 = %.6g Hz; " + "%d port%s, %d sideband%s\n", + (stepType == 1) ? "dec" : (stepType == 2) ? "oct" : "lin", + fstart, fstop, f0, np, (np == 1) ? "" : "s", nsb, (nsb == 1) ? "" : "s"); + + for (freq = fstart; freq <= fstop * (1.0 + 1e-9); ) { + int ok = 1; + /* excite each port -> one column of A (sb0) and of each B^(k) */ + for (i = 1; i <= np && ok; i++) { + VSRCinstance *pi = (VSRCinstance *) ckt->CKTrfPorts[i - 1]; + if (psp_solve_port(&hd, f0, freq, pi->VSRCbranch, Xr, Xi)) { + ok = 0; + break; + } + for (j = 1; j <= np; j++) { + VSRCinstance *pj = (VSRCinstance *) ckt->CKTrfPorts[j - 1]; + double ki = pj->VSRCki, z0 = pj->VSRCportZ0; + int pn = pj->VSRCposNode, nn = pj->VSRCnegNode, br = pj->VSRCbranch; + for (s = 0; s < nsb; s++) { + int k = s - Ksb; + size_t blk = (size_t)(k + M) * (size_t)N; + double Vr = 0.0, Vi = 0.0, Ir, Ii; + cplx a, b; + if (pn > 0) { Vr += Xr[blk + (size_t)(pn - 1)]; Vi += Xi[blk + (size_t)(pn - 1)]; } + if (nn > 0) { Vr -= Xr[blk + (size_t)(nn - 1)]; Vi -= Xi[blk + (size_t)(nn - 1)]; } + Ir = -Xr[blk + (size_t)(br - 1)]; + Ii = -Xi[blk + (size_t)(br - 1)]; + b.re = ki * (Vr - z0 * Ir); + b.im = ki * (Vi - z0 * Ii); + setc(Bmat[s], j - 1, i - 1, b); + if (k == 0) { + a.re = ki * (Vr + z0 * Ir); + a.im = ki * (Vi + z0 * Ii); + setc(Amat, j - 1, i - 1, a); + } + } + } + } + if (ok && (Ainv = cinverse(Amat)) != NULL) { + IFvalue freqData, valueData; + IFcomplex *data = TMALLOC(IFcomplex, numNames); + int idx = 0; + freqData.rValue = freq; + valueData.v.numValue = numNames; + valueData.v.vec.cVec = data; + for (s = 0; s < nsb; s++) { + CMat *S = cmultiply(Bmat[s], Ainv); /* S^(k) = B^(k) * A^-1 */ + for (j = 0; j < np; j++) + for (i = 0; i < np; i++) { + cplx sij = getcplx(S, j, i); + data[idx].real = sij.re; + data[idx].imag = sij.im; + idx++; + } + freecmat(S); + } + SPfrontEnd->OUTpData(plot, &freqData, &valueData); + FREE(data); + freecmat(Ainv); + Ainv = NULL; + } + if (stepType == 0) { if (npts <= 1) break; freq += linstep; } + else { freq *= mult; } + } + + SPfrontEnd->OUTendPlot(plot); + for (s = 0; s < nsb; s++) + freecmat(Bmat[s]); + FREE(Bmat); + freecmat(Amat); + FREE(Xr); FREE(Xi); + pac_free_harmonics(&hd); +} +#endif + + int DCpss(CKTcircuit *ckt, int restart) /* forced restart flag */ @@ -2050,6 +2260,13 @@ shootingexit: * report the input->output transfer at each sideband. */ if (job->PSSdoPXF) pxf_sweep (ckt, job) ; + +#ifdef RFSPICE + /* Enhancement-132: for a .psp card, inject at each RF port through + * the conversion matrix and emit the periodic S-parameters. */ + if (job->PSSdoPSP) + psp_sweep (ckt, job) ; +#endif } /****************************/ diff --git a/src/spicelib/analysis/psssetp.c b/src/spicelib/analysis/psssetp.c index b8785daaf..7ea773af2 100644 --- a/src/spicelib/analysis/psssetp.c +++ b/src/spicelib/analysis/psssetp.c @@ -91,6 +91,11 @@ PSSsetParm(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value) job->PSSpnCyclo = value->iValue; break; + /* Enhancement-132: periodic S-parameters flag */ + case PSP_DOPSP: + job->PSSdoPSP = value->iValue; + break; + default: return(E_BADPARM); } @@ -118,7 +123,8 @@ static IFparm PSSparms[] = { { "pnoise_insrc", PNOISE_INSRC, IF_SET|IF_STRING, "pnoise input source (for the input-referred spectrum)" }, { "pxf", PXF_DO, IF_SET|IF_INTEGER, "run a periodic transfer-function sweep after PSS" }, { "pxf_out", PXF_OUT, IF_SET|IF_STRING, "pxf output node" }, - { "pnoise_cyclo", PNOISE_CYCLO, IF_SET|IF_INTEGER, "pnoise cyclostationary mode (per-sample bias, period average)" } + { "pnoise_cyclo", PNOISE_CYCLO, IF_SET|IF_INTEGER, "pnoise cyclostationary mode (per-sample bias, period average)" }, + { "psp", PSP_DOPSP, IF_SET|IF_INTEGER, "run a periodic S-parameter sweep after PSS (needs RF ports)" } }; SPICEanalysis PSSinfo = { diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index 21220c8af..a61b7bd7b 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -786,6 +786,89 @@ dot_pac(char *line, void *ckt, INPtables *tab, struct card *current, return (0); } +/* Enhancement-132: Periodic S-parameters (PSP). Runs PSS then, for each RF port, + * injects through the harmonic conversion matrix and forms the periodic scattering + * matrix vs input frequency. Same PSS params + sweep tail as .pac; the excitation is + * the netlist's RF ports (portnum/z0), so there is no output node. */ +static int +dot_psp(char *line, void *ckt, INPtables *tab, struct card *current, + void *task, void *gnode, JOB *foo) +{ + int error; /* error code temporary */ + IFvalue ptemp; /* a value structure to package resistance into */ + IFvalue *parm; /* a pointer to a value struct for function returns */ + char *nname; /* the oscNode name */ + CKTnode *nnode; /* the oscNode node */ + int which; /* which analysis we are performing */ + char *steptype; /* psp sweep type: dec/oct/lin */ + + NG_IGNORE(gnode); + NG_IGNORE(current); + + /* .psp Fguess StabTime OscNode Points Harmonics SC_iter Steady_coeff + * NumPts Fstart Fstop [maxsideband] */ + which = ft_find_analysis("PSS"); + if (which == -1) { + LITERR("Periodic S-parameter (PSP) analysis unsupported.\n"); + return (0); + } + IFC(newAnalysis, (ckt, which, "Periodic S-parameter Analysis", &foo, task)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* Fguess */ + GCA(INPapName, (ckt, which, foo, "fguess", parm)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* StabTime */ + GCA(INPapName, (ckt, which, foo, "stabtime", parm)); + + INPgetNetTok(&line, &nname, 0); + INPtermInsert(ckt, &nname, tab, &nnode); + ptemp.nValue = nnode; + GCA(INPapName, (ckt, which, foo, "oscnode", &ptemp)); /* OscNode given as string */ + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* PSS points */ + GCA(INPapName, (ckt, which, foo, "points", parm)); + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* PSS harmonics */ + GCA(INPapName, (ckt, which, foo, "harmonics", parm)); + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* SC iterations */ + GCA(INPapName, (ckt, which, foo, "sc_iter", parm)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* Steady coefficient */ + GCA(INPapName, (ckt, which, foo, "steady_coeff", parm)); + + /* PSP sweep tail: NumPts Fstart Fstop */ + INPgetTok(&line, &steptype, 1); + ptemp.iValue = (strcmp(steptype, "dec") == 0) ? 1 : + (strcmp(steptype, "oct") == 0) ? 2 : 0; /* default LIN */ + tfree(steptype); + GCA(INPapName, (ckt, which, foo, "pac_step", &ptemp)); + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* number of points */ + GCA(INPapName, (ckt, which, foo, "pac_points", parm)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* fstart */ + GCA(INPapName, (ckt, which, foo, "pac_fstart", parm)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* fstop */ + GCA(INPapName, (ckt, which, foo, "pac_fstop", parm)); + + { /* optional trailing maxsideband: output conversion sidebands each side */ + char *p = line; + while (*p == ' ' || *p == '\t') + p++; + if (*p) { + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); + GCA(INPapName, (ckt, which, foo, "pac_maxsb", parm)); + } + } + + ptemp.iValue = 1; /* enable the PSP sweep */ + GCA(INPapName, (ckt, which, foo, "psp", &ptemp)); + + return (0); +} + /* Enhancement-124: Periodic noise (PNOISE). Runs PSS then folds each device's * noise through the conversion-matrix adjoint over all sidebands. Reuses the PSS * analysis (like .pac) with the pnoise output node, input source, and sweep set. */ @@ -1164,6 +1247,10 @@ INP2dot(CKTcircuit *ckt, INPtables *tab, struct card *current, TSKtask *task, CK } else if ((strcmp(token, ".pxf") == 0)) { rtn = dot_pxf(line, ckt, tab, current, task, gnode, foo); goto quit; + /* Enhancement-132: Periodic S-parameters */ + } else if ((strcmp(token, ".psp") == 0)) { + rtn = dot_psp(line, ckt, tab, current, task, gnode, foo); + goto quit; #endif #ifdef RFSPICE }