diff --git a/src/include/ngspice/pssdefs.h b/src/include/ngspice/pssdefs.h index 7a33fd5db..f63d6f95c 100644 --- a/src/include/ngspice/pssdefs.h +++ b/src/include/ngspice/pssdefs.h @@ -58,6 +58,13 @@ typedef struct { int PSSdoPnoise; /* 1 if this job runs a pnoise sweep */ CKTnode *PnOutNode; /* pnoise output node (reference = ground) */ IFuid PnInSrc; /* input source name, for the input-referred spectrum */ + + /* Enhancement-125: periodic transfer function (.pxf). The ADJOINT counterpart of + * PAC: solve Hᵀ Ψ = e_{out,0} and dot Ψ with the netlist AC-source pattern to get + * the transfer from the input to a fixed output at each sideband. Reuses the PAC + * sweep fields + PACmaxSideband. */ + int PSSdoPXF; /* 1 if this job runs a PXF sweep */ + CKTnode *PxOutNode; /* PXF output node (reference = ground) */ } PSSan; enum { @@ -78,6 +85,8 @@ enum { PNOISE_DO, PNOISE_OUT, PNOISE_INSRC, + PXF_DO, + PXF_OUT, }; #endif diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c index ab1325584..717e9e066 100644 --- a/src/spicelib/analysis/dcpss.c +++ b/src/spicelib/analysis/dcpss.c @@ -827,6 +827,109 @@ pnoise_sweep(CKTcircuit *ckt, PSSan *job) } +/* Enhancement-125: periodic transfer function (.pxf). The adjoint counterpart of + * .pac: solve Hᵀ Ψ = e_{out,0} once per frequency and dot each sideband block of Ψ + * with the netlist AC-source pattern B_0 to get the transfer from the input to the + * fixed output at each sideband, xf_k = Σ_j Ψ_k(j)·B0(j). By the identity + * (H⁻¹B)_out = (H⁻ᵀe_out)ᵀB, the sideband-0 transfer equals the PAC response at the + * output exactly -- the reciprocity cross-check. Emits xf (sideband 0) plus + * xf_usb/xf_lsb conversion transfers as a complex plot vs frequency. */ +static void +pxf_sweep(CKTcircuit *ckt, PSSan *job) +{ + int N = job->PSSopMsize, outNode = job->PxOutNode ? job->PxOutNode->number : 0; + int M, Ksb, nsb, s, k, j, error, stepType = job->PACstepType, np = job->PACpoints; + double f0 = job->PSSopFreq, fstart = job->PACfStart, fstop = job->PACfStop; + double freq, mult, linstep; + struct pac_harm hd; + double *Psr, *Psi; + IFuid freqUid, *outNames = NULL; + runDesc *plot = NULL; + char nm[64]; + + if (outNode <= 0 || outNode > N || f0 <= 0.0 || fstart <= 0.0 || + fstop < fstart || np < 1) + return; + M = pac_choose_M(ckt, job); + if (M < 1) { + fprintf(stderr, "PXF: conversion matrix too large -- sweep skipped\n"); + return; + } + if (pac_extract_harmonics(ckt, job, M, &hd)) + return; + if (!hd.has_src) { + fprintf(stderr, "PXF: no netlist AC source found -- give the input source an " + "AC value; sweep skipped\n"); + pac_free_harmonics(&hd); + return; + } + + Ksb = job->PACmaxSideband; + if (Ksb < 0) Ksb = 0; + if (Ksb > M) Ksb = M; + nsb = 2 * Ksb + 1; + + /* one transfer vector per output sideband: xf (sb0), xf_usb, xf_lsb */ + outNames = TMALLOC(IFuid, nsb); + for (s = 0; s < nsb; s++) { + k = s - Ksb; + if (k == 0) + (void) snprintf(nm, sizeof(nm), "xf"); + else + (void) snprintf(nm, sizeof(nm), "xf_%csb%d", (k > 0) ? 'u' : 'l', abs(k)); + SPfrontEnd->IFnewUid(ckt, &outNames[s], NULL, nm, UID_OTHER, NULL); + } + + SPfrontEnd->IFnewUid(ckt, &freqUid, NULL, "frequency", UID_OTHER, NULL); + error = SPfrontEnd->OUTpBeginPlot(ckt, ckt->CKTcurJob, "PXF Analysis", + freqUid, IF_REAL, nsb, outNames, + IF_COMPLEX, &plot); + tfree(outNames); + if (error) { pac_free_harmonics(&hd); return; } + if (stepType != 0) + SPfrontEnd->OUTattributes(plot, NULL, OUT_SCALE_LOG, NULL); + + Psr = TMALLOC(double, hd.Ntot); + Psi = TMALLOC(double, hd.Ntot); + mult = (stepType == 1) ? pow(10.0, 1.0 / np) : + (stepType == 2) ? pow(2.0, 1.0 / np) : 0.0; + linstep = (np > 1) ? (fstop - fstart) / (np - 1) : 0.0; + + fprintf(stderr, "PXF sweep: %s from %.6g to %.6g Hz around f0 = %.6g Hz; " + "output node %d; %d sideband%s (adjoint)\n", + (stepType == 1) ? "dec" : (stepType == 2) ? "oct" : "lin", + fstart, fstop, f0, outNode, nsb, (nsb == 1) ? "" : "s"); + + for (freq = fstart; freq <= fstop * (1.0 + 1e-9); ) { + if (pac_solve_adjoint(&hd, f0, freq, outNode, Psr, Psi) == 0) { + IFvalue freqData, valData; + IFcomplex *data = TMALLOC(IFcomplex, nsb); + freqData.rValue = freq; + valData.v.numValue = nsb; + valData.v.vec.cVec = data; + for (s = 0; s < nsb; s++) { + size_t blk = (size_t)(s - Ksb + M) * (size_t)N; + double xr = 0.0, xi = 0.0; /* xf_k = sum_j Psi_k(j) * B0(j) */ + for (j = 0; j < N; j++) { + xr += Psr[blk + (size_t)j] * hd.B0r[j] - Psi[blk + (size_t)j] * hd.B0i[j]; + xi += Psr[blk + (size_t)j] * hd.B0i[j] + Psi[blk + (size_t)j] * hd.B0r[j]; + } + data[s].real = xr; + data[s].imag = xi; + } + SPfrontEnd->OUTpData(plot, &freqData, &valData); + FREE(data); + } + if (stepType == 0) { if (np <= 1) break; freq += linstep; } + else { freq *= mult; } + } + + SPfrontEnd->OUTendPlot(plot); + FREE(Psr); FREE(Psi); + pac_free_harmonics(&hd); +} + + int DCpss(CKTcircuit *ckt, int restart) /* forced restart flag */ @@ -1855,6 +1958,11 @@ shootingexit: * through the conversion matrix to get the output noise spectrum. */ if (job->PSSdoPnoise) pnoise_sweep (ckt, job) ; + + /* Enhancement-125: for a .pxf card, solve the conversion adjoint and + * report the input->output transfer at each sideband. */ + if (job->PSSdoPXF) + pxf_sweep (ckt, job) ; } /****************************/ diff --git a/src/spicelib/analysis/psssetp.c b/src/spicelib/analysis/psssetp.c index 13d4a420e..d285bfc0b 100644 --- a/src/spicelib/analysis/psssetp.c +++ b/src/spicelib/analysis/psssetp.c @@ -78,6 +78,14 @@ PSSsetParm(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value) job->PnInSrc = value->uValue; break; + /* Enhancement-125: pxf parameters (.pxf card) */ + case PXF_DO: + job->PSSdoPXF = value->iValue; + break; + case PXF_OUT: + job->PxOutNode = value->nValue; + break; + default: return(E_BADPARM); } @@ -102,7 +110,9 @@ static IFparm PSSparms[] = { { "pac_maxsb", PAC_MAXSB, IF_SET|IF_INTEGER, "PAC output conversion sidebands each side (0 = sideband 0 only)" }, { "pnoise", PNOISE_DO, IF_SET|IF_INTEGER, "run a periodic-noise sweep after PSS" }, { "pnoise_out", PNOISE_OUT, IF_SET|IF_STRING, "pnoise output node" }, - { "pnoise_insrc", PNOISE_INSRC, IF_SET|IF_STRING, "pnoise input source (for the input-referred spectrum)" } + { "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" } }; SPICEanalysis PSSinfo = { diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index 18f684684..e13d6bafb 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -859,6 +859,85 @@ dot_pnoise(char *line, void *ckt, INPtables *tab, struct card *current, return (0); } + +/* Enhancement-125: Periodic transfer function (PXF). The adjoint of PAC: runs PSS + * then solves Hᵀ Ψ = e_{out,0} and dots Ψ with the netlist AC-source pattern to get + * the transfer from the input to a fixed output at each sideband. Reuses the PSS + * analysis (like .pac). */ +static int +dot_pxf(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; /* a node name */ + CKTnode *nnode; /* a node pointer */ + int which; /* which analysis we are performing */ + char *steptype; /* pxf sweep type: dec/oct/lin */ + + NG_IGNORE(gnode); + NG_IGNORE(current); + + /* .pxf Fguess StabTime OscNode Points Harmonics SC_iter Steady_coeff + * OutNode NumPts Fstart Fstop [maxsideband] */ + which = ft_find_analysis("PSS"); + if (which == -1) { + LITERR("Periodic transfer-function (PXF) analysis unsupported.\n"); + return (0); + } + IFC(newAnalysis, (ckt, which, "Periodic Transfer Function 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); /* OscNode */ + INPtermInsert(ckt, &nname, tab, &nnode); + ptemp.nValue = nnode; + GCA(INPapName, (ckt, which, foo, "oscnode", &ptemp)); + 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)); + + INPgetNetTok(&line, &nname, 0); /* OutNode */ + INPtermInsert(ckt, &nname, tab, &nnode); + ptemp.nValue = nnode; + GCA(INPapName, (ckt, which, foo, "pxf_out", &ptemp)); + + /* sweep tail: NumPts Fstart Fstop [maxsideband] */ + INPgetTok(&line, &steptype, 1); + ptemp.iValue = (strcmp(steptype, "dec") == 0) ? 1 : + (strcmp(steptype, "oct") == 0) ? 2 : 0; + 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 */ + 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 pxf sweep */ + GCA(INPapName, (ckt, which, foo, "pxf", &ptemp)); + + return (0); +} #endif @@ -1064,6 +1143,10 @@ INP2dot(CKTcircuit *ckt, INPtables *tab, struct card *current, TSKtask *task, CK } else if ((strcmp(token, ".pnoise") == 0)) { rtn = dot_pnoise(line, ckt, tab, current, task, gnode, foo); goto quit; + /* Enhancement-125: Periodic transfer function */ + } else if ((strcmp(token, ".pxf") == 0)) { + rtn = dot_pxf(line, ckt, tab, current, task, gnode, foo); + goto quit; #endif #ifdef RFSPICE }