pa-125

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

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<k>/xf_lsb<k> 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<k>, xf_lsb<k> */
+    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) ;
             }
             /****************************/
 

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  = {

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 <DEC|OCT|LIN> 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: <DEC|OCT|LIN> 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
     }