diff --git a/src/include/ngspice/pssdefs.h b/src/include/ngspice/pssdefs.h
index f63d6f95c..3fe2fe83f 100644
--- a/src/include/ngspice/pssdefs.h
+++ b/src/include/ngspice/pssdefs.h
@@ -58,6 +58,9 @@ 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 */
+    int      PSSpnCyclo;     /* Enhancement-126: 1 = cyclostationary noise (per-sample
+                              * bias + time-domain transfer, averaged over the period);
+                              * 0 = stationary (noise PSD at one operating point) */
 
     /* 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
@@ -87,6 +90,7 @@ enum {
     PNOISE_INSRC,
     PXF_DO,
     PXF_OUT,
+    PNOISE_CYCLO,
 };
 
 #endif
diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c
index 717e9e066..a1576e41e 100644
--- a/src/spicelib/analysis/dcpss.c
+++ b/src/spicelib/analysis/dcpss.c
@@ -769,10 +769,97 @@ pnoise_sweep(CKTcircuit *ckt, PSSan *job)
     linstep = (np > 1) ? (fstop - fstart) / (np - 1) : 0.0;
 
     fprintf(stderr, "PNOISE sweep: %s from %.6g to %.6g Hz around f0 = %.6g Hz; "
-                    "output node %d; folding %d sidebands\n",
+                    "output node %d; folding %d sidebands%s\n",
             (stepType == 1) ? "dec" : (stepType == 2) ? "oct" : "lin",
-            fstart, fstop, f0, outNode, 2*M + 1);
+            fstart, fstop, f0, outNode, 2*M + 1,
+            job->PSSpnCyclo ? "; cyclostationary" : "");
 
+    if (job->PSSpnCyclo) {
+        /* Enhancement-126: cyclostationary noise. The device noise PSD S(t) varies
+         * along the PSS period, and its harmonics couple sidebands. Using the
+         * identity onoise = (1/P) Σ_s S(t_s)·|ΔA_s|², where A_s(j) = Σ_k Ψ_k(j)·
+         * exp(j·2π·k·s/P) is the inverse-DFT of the sideband adjoint transfers, this
+         * is computed by evaluating each device's noise at every sample's bias
+         * (CKTload per sample) and folding through the time-domain transfer, then
+         * averaging over the period. Reduces to the stationary case (and hence
+         * .noise) when S(t) is constant, by Parseval. */
+        long   P = job->PSSopPoints, s;
+        int    Nf = 0, fi, c;
+        double *freqs, *onz, *Pr_all, *Pi_all;
+
+        for (freq = fstart; freq <= fstop * (1.0 + 1e-9); ) {   /* count points */
+            Nf++;
+            if (stepType == 0) { if (np <= 1) break; freq += linstep; }
+            else               { freq *= mult; }
+        }
+        freqs  = TMALLOC(double, Nf);
+        onz    = TMALLOC(double, Nf);
+        Pr_all = TMALLOC(double, (size_t)Nf * (size_t)hd.Ntot);
+        Pi_all = TMALLOC(double, (size_t)Nf * (size_t)hd.Ntot);
+        c = 0;
+        for (freq = fstart; freq <= fstop * (1.0 + 1e-9); ) {   /* fill + adjoints */
+            freqs[c] = freq;  onz[c] = 0.0;
+            if (pac_solve_adjoint(&hd, f0, freq, outNode, Psr, Psi) != 0) {
+                memset(Psr, 0, (size_t)hd.Ntot * sizeof(double));
+                memset(Psi, 0, (size_t)hd.Ntot * sizeof(double));
+            }
+            memcpy(Pr_all + (size_t)c * (size_t)hd.Ntot, Psr, (size_t)hd.Ntot * sizeof(double));
+            memcpy(Pi_all + (size_t)c * (size_t)hd.Ntot, Psi, (size_t)hd.Ntot * sizeof(double));
+            c++;
+            if (stepType == 0) { if (np <= 1) break; freq += linstep; }
+            else               { freq *= mult; }
+        }
+
+        for (s = 0; s < P; s++) {   /* evaluate device noise at each sample's bias */
+            double ang0 = 2.0 * M_PI * (double)s / (double)P;
+            for (i = 1; i <= N; i++)
+                ckt->CKTrhsOld[i] = job->PSSopVoltages[(i - 1) + s * N];
+            ckt->CKTrhsOld[0] = 0.0;
+            if (ns > 0)
+                memcpy(ckt->CKTstate0, job->PSSopStates + (size_t)s * (size_t)ns,
+                       (size_t)ns * sizeof(double));
+            ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG;
+            CKTload(ckt);
+            for (fi = 0; fi < Nf; fi++) {
+                double dens = 0.0;
+                double *pr = Pr_all + (size_t)fi * (size_t)hd.Ntot;
+                double *pi = Pi_all + (size_t)fi * (size_t)hd.Ntot;
+                for (j = 1; j <= N; j++) {   /* A_s(j) = IDFT_k Ψ_k(j) */
+                    double ar = 0.0, ai = 0.0;
+                    for (k = -M; k <= M; k++) {
+                        size_t idx = (size_t)(k + M) * (size_t)N + (size_t)(j - 1);
+                        double cs = cos((double)k * ang0), sn = sin((double)k * ang0);
+                        ar += pr[idx] * cs - pi[idx] * sn;
+                        ai += pr[idx] * sn + pi[idx] * cs;
+                    }
+                    ckt->CKTrhs[j] = ar;  ckt->CKTirhs[j] = ai;
+                }
+                ckt->CKTrhs[0] = 0.0;  ckt->CKTirhs[0] = 0.0;
+                data.freq = freqs[fi];  data.delFreq = 0.0;  data.prtSummary = FALSE;
+                for (i = 0; i < DEVmaxnum; i++)
+                    if (DEVices[i] && DEVices[i]->DEVnoise && ckt->CKThead[i])
+                        DEVices[i]->DEVnoise(N_DENS, N_CALC, ckt->CKThead[i],
+                                             ckt, &data, &dens);
+                onz[fi] += dens;
+            }
+        }
+
+        for (fi = 0; fi < Nf; fi++) {   /* period-average, gain, output */
+            double onoise = onz[fi] / (double)P, gain2 = 1.0, gsi;
+            IFvalue refVal, valData;
+            double out[2];
+            if (hd.has_src && pac_solve_at(&hd, f0, freqs[fi], outNode, 1, Xr, Xi) == 0) {
+                size_t oidx = (size_t)M * (size_t)N + (size_t)(outNode - 1);
+                gain2 = Xr[oidx] * Xr[oidx] + Xi[oidx] * Xi[oidx];
+            }
+            gsi = 1.0 / MAX(gain2, N_MINGAIN);
+            out[0] = onoise;  out[1] = onoise * gsi;
+            refVal.rValue = freqs[fi];
+            valData.v.numValue = 2;  valData.v.vec.rVec = out;
+            SPfrontEnd->OUTpData(plot, &refVal, &valData);
+        }
+        FREE(freqs); FREE(onz); FREE(Pr_all); FREE(Pi_all);
+    } else
     for (freq = fstart; freq <= fstop * (1.0 + 1e-9); ) {
         double onoise = 0.0, gain2 = 1.0, gsi;
 
diff --git a/src/spicelib/analysis/psssetp.c b/src/spicelib/analysis/psssetp.c
index d285bfc0b..b8785daaf 100644
--- a/src/spicelib/analysis/psssetp.c
+++ b/src/spicelib/analysis/psssetp.c
@@ -86,6 +86,11 @@ PSSsetParm(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value)
         job->PxOutNode = value->nValue;
         break;
 
+    /* Enhancement-126: cyclostationary-noise flag */
+    case PNOISE_CYCLO:
+        job->PSSpnCyclo = value->iValue;
+        break;
+
     default:
         return(E_BADPARM);
     }
@@ -112,7 +117,8 @@ static IFparm PSSparms[] = {
     { "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)" },
     { "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" }
+    { "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)" }
 };
 
 SPICEanalysis PSSinfo  = {
diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c
index e13d6bafb..21220c8af 100644
--- a/src/spicelib/parser/inp2dot.c
+++ b/src/spicelib/parser/inp2dot.c
@@ -854,6 +854,23 @@ dot_pnoise(char *line, void *ckt, INPtables *tab, struct card *current,
     parm = INPgetValue(ckt, &line, IF_REAL, tab);		/* fstop */
     GCA(INPapName, (ckt, which, foo, "pac_fstop", parm));
 
+    {   /* Enhancement-126: optional trailing "cyclo" keyword */
+        char *p = line;
+        while (*p == ' ' || *p == '\t')
+            p++;
+        if (*p) {
+            char *word;
+            INPgetTok(&line, &word, 1);
+            if (strcmp(word, "cyclo") == 0) {
+                ptemp.iValue = 1;
+                GCA(INPapName, (ckt, which, foo, "pnoise_cyclo", &ptemp));
+            } else {
+                fprintf(stderr, "Error: unknown parameter %s on .pnoise - ignored\n", word);
+            }
+            tfree(word);
+        }
+    }
+
     ptemp.iValue = 1;						/* enable the pnoise sweep */
     GCA(INPapName, (ckt, which, foo, "pnoise", &ptemp));