luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

pa-124

This commit is contained in:
Meisam Bahadori 2026-07-26 11:11:15 +02:00 committed by Holger Vogt
parent da174d9950
commit 951a6d4a14
4 changed files with 307 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/include/ngspice/pssdefs.h
View File

@ -51,6 +51,13 @@ typedef struct {
int PACpoints; /* points per decade/octave, or total for linear */
int PACstepType; /* 0 = LIN, 1 = DEC, 2 = OCT */
int PACmaxSideband; /* number of output conversion sidebands each side (0 = sb0 only) */
/* Enhancement-124: periodic noise (.pnoise). Runs PSS, then folds each device's
* noise through the conversion-matrix adjoint over all sidebands to get the
* output noise spectrum. Reuses the PAC sweep fields for the frequency sweep. */
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 */
} PSSan;
enum {
@ -68,6 +75,9 @@ enum {
PAC_POINTS,
PAC_STEPTYPE,
PAC_MAXSB,
PNOISE_DO,
PNOISE_OUT,
PNOISE_INSRC,
};
#endif

204
src/spicelib/analysis/dcpss.c
View File

@ -12,6 +12,7 @@
#include "ngspice/devdefs.h" /* Enhancement-120: DEVbindCSCComplex for the KLU AC stamp */
#include "ngspice/smpdefs.h" /* Enhancement-120: SMPfindElt to read the Jacobian */
#include "ngspice/spmatrix.h" /* Enhancement-120: spSetComplex (Sparse complex mode) */
#include "ngspice/noisedef.h" /* Enhancement-124: NOISEAN/Ndata + CKTnoise for pnoise */
#include "ngspice/sperror.h"
#include "ngspice/fteext.h"
@ -628,6 +629,204 @@ pac_sweep(CKTcircuit *ckt, PSSan *job)
}
/* Enhancement-124: solve the ADJOINT conversion system Hᵀ Psi = e_{out,0}. Psi
* then holds, for every (node j, sideband k), the transfer from a unit injection
* at (j,k) to the output at sideband 0 -- the conversion transimpedance the noise
* folding needs. Assembles Hᵀ (the transpose of the pac_solve_at matrix) and puts a
* unit at the output node in the 0-th sideband. Returns 0 on success, 1 if singular. */
static int
pac_solve_adjoint(struct pac_harm *hd, double f0, double f_in, int outNode,
double *Psr, double *Psi)
{
int N = hd->N, M = hd->M, H = hd->H, nnz = hd->nnz, Ntot = hd->Ntot;
int ni, mi, n, mm, e, rc;
double *Ar, *Ai;
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++) {
n = ni - M;
for (mi = 0; mi <= 2*M; mi++) {
int dm;
double omega;
mm = mi - M;
dm = n - mm;
omega = 2.0 * M_PI * (f_in + (double)mm * f0);
for (e = 0; e < nnz; e++) {
double gr, gi, cr, ci;
size_t hi = (size_t)e * (size_t)(H + 1) + (size_t)abs(dm);
gr = hd->Gmr[hi]; gi = hd->Gmi[hi];
cr = hd->Cmr[hi]; ci = hd->Cmi[hi];
if (dm < 0) { gi = -gi; ci = -ci; }
{
double er = gr - omega * ci;
double ei = gi + omega * cr;
size_t row = (size_t)ni * (size_t)N + (size_t)(hd->rr[e] - 1);
size_t col = (size_t)mi * (size_t)N + (size_t)(hd->cc[e] - 1);
Ar[col * (size_t)Ntot + row] += er; /* transpose: [col][row] */
Ai[col * (size_t)Ntot + row] += ei;
}
}
}
}
memset(Psr, 0, (size_t)Ntot * sizeof(double));
memset(Psi, 0, (size_t)Ntot * sizeof(double));
Psr[(size_t)M * (size_t)N + (size_t)(outNode - 1)] = 1.0;
rc = pss_csolve(Ntot, Ar, Ai, Psr, Psi);
FREE(Ar); FREE(Ai);
return rc;
}
/* Enhancement-124: periodic noise (.pnoise). Runs off the retained operating point:
* folds every device's noise through the conversion-matrix adjoint over all
* sidebands to get the output noise spectrum. The device noise routines
* (DEVnoise/NevalSrc, OSDI load_noise) compute S*|dTransimp|^2 reading the
* transimpedance from CKTrhs/CKTirhs -- so loading the sideband-k adjoint into
* CKTrhs and summing over k = -M..M folds the noise exactly. A local NOISEAN job
* gives those routines their expected context. For a linear (block-diagonal)
* circuit only sideband 0 contributes, so the result reduces to ordinary .noise. */
static void
pnoise_sweep(CKTcircuit *ckt, PSSan *job)
{
int N = job->PSSopMsize, ns = job->PSSopNumStates;
int outNode = job->PnOutNode ? job->PnOutNode->number : 0;
int M, i, j, k, np = job->PACpoints, stepType = job->PACstepType, error;
double f0 = job->PSSopFreq, fstart = job->PACfStart, fstop = job->PACfStop;
double freq, mult, linstep;
struct pac_harm hd;
double *Psr, *Psi, *Xr, *Xi;
NOISEAN nj;
Ndata data;
JOB *oldJob;
IFuid freqUid, nlist[2];
runDesc *plot = NULL;
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, "PNOISE: conversion matrix too large -- sweep skipped\n");
return;
}
if (pac_extract_harmonics(ckt, job, M, &hd))
return;
/* set the device bias to the (sample-0) operating point so each noise PSD
* (conductances, dc currents) is evaluated at the periodic operating point. */
for (i = 1; i <= N; i++)
ckt->CKTrhsOld[i] = job->PSSopVoltages[(i - 1)];
ckt->CKTrhsOld[0] = 0.0;
if (ns > 0)
memcpy(ckt->CKTstate0, job->PSSopStates, (size_t)ns * sizeof(double));
ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG;
CKTload(ckt);
/* output plot (onoise/inoise spectrum vs frequency), opened while CKTcurJob is
* still the persistent PSS job. */
SPfrontEnd->IFnewUid(ckt, &freqUid, NULL, "frequency", UID_OTHER, NULL);
SPfrontEnd->IFnewUid(ckt, &nlist[0], NULL, "onoise_spectrum", UID_OTHER, NULL);
SPfrontEnd->IFnewUid(ckt, &nlist[1], NULL, "inoise_spectrum", UID_OTHER, NULL);
error = SPfrontEnd->OUTpBeginPlot(ckt, ckt->CKTcurJob, "PNoise Analysis",
freqUid, IF_REAL, 2, nlist, IF_REAL, &plot);
if (error) { pac_free_harmonics(&hd); return; }
if (stepType != 0)
SPfrontEnd->OUTattributes(plot, NULL, OUT_SCALE_LOG, NULL);
/* a minimal NOISEAN context for the device noise routines (they cast
* CKTcurJob to NOISEAN* and read NStpsSm / NstartFreq). */
memset(&nj, 0, sizeof(nj));
nj.output = job->PnOutNode;
nj.outputRef = job->PnOutNode;
nj.input = job->PnInSrc;
nj.NstartFreq = fstart;
nj.NstopFreq = fstop;
nj.NnumSteps = np;
nj.NstpType = stepType;
nj.NStpsSm = 0; /* no per-device summary vectors */
nj.JOBname = "pnoise";
memset(&data, 0, sizeof(data));
data.prtSummary = FALSE; /* keep the routines from writing outpVector */
oldJob = ckt->CKTcurJob;
ckt->CKTcurJob = (JOB *) &nj;
/* let each device set up its noise state (a no-op naming pass with NStpsSm=0) */
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);
}
Psr = TMALLOC(double, hd.Ntot); Psi = TMALLOC(double, hd.Ntot);
Xr = TMALLOC(double, hd.Ntot); Xi = 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, "PNOISE sweep: %s from %.6g to %.6g Hz around f0 = %.6g Hz; "
"output node %d; folding %d sidebands\n",
(stepType == 1) ? "dec" : (stepType == 2) ? "oct" : "lin",
fstart, fstop, f0, outNode, 2*M + 1);
for (freq = fstart; freq <= fstop * (1.0 + 1e-9); ) {
double onoise = 0.0, gain2 = 1.0, gsi;
data.freq = freq;
data.delFreq = 0.0; /* density only -- we do not integrate here */
data.prtSummary = FALSE;
/* transfer from every (node, sideband) to the output at sideband 0 */
if (pac_solve_adjoint(&hd, f0, freq, outNode, Psr, Psi) == 0) {
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 (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);
onoise += dens; /* sum device noise over sidebands */
}
}
/* input-referred: divide by the source->output conversion gain (sideband 0) */
if (hd.has_src && pac_solve_at(&hd, f0, freq, 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);
{
IFvalue refVal, valData;
double out[2];
out[0] = onoise; /* output noise density (V^2/Hz) */
out[1] = onoise * gsi; /* input-referred density */
refVal.rValue = freq;
valData.v.numValue = 2;
valData.v.vec.rVec = out;
SPfrontEnd->OUTpData(plot, &refVal, &valData);
}
if (stepType == 0) { if (np <= 1) break; freq += linstep; }
else { freq *= mult; }
}
SPfrontEnd->OUTendPlot(plot);
ckt->CKTcurJob = oldJob;
FREE(Psr); FREE(Psi); FREE(Xr); FREE(Xi);
pac_free_harmonics(&hd);
}
int
DCpss(CKTcircuit *ckt,
int restart) /* forced restart flag */
@ -1651,6 +1850,11 @@ shootingexit:
* emit the sideband-0 node responses as a complex plot. */
if (job->PSSdoPAC)
pac_sweep (ckt, job) ;
/* Enhancement-124: for a .pnoise card, fold each device's noise
* through the conversion matrix to get the output noise spectrum. */
if (job->PSSdoPnoise)
pnoise_sweep (ckt, job) ;
}
/****************************/

16
src/spicelib/analysis/psssetp.c
View File

@ -67,6 +67,17 @@ PSSsetParm(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value)
job->PACmaxSideband = value->iValue;
break;
/* Enhancement-124: pnoise parameters (.pnoise card) */
case PNOISE_DO:
job->PSSdoPnoise = value->iValue;
break;
case PNOISE_OUT:
job->PnOutNode = value->nValue;
break;
case PNOISE_INSRC:
job->PnInSrc = value->uValue;
break;
default:
return(E_BADPARM);
}
@ -88,7 +99,10 @@ static IFparm PSSparms[] = {
{ "pac_fstop", PAC_FSTOP, IF_SET|IF_REAL, "PAC input sweep stop frequency" },
{ "pac_points", PAC_POINTS, IF_SET|IF_INTEGER, "PAC points per decade/octave (or total for linear)" },
{ "pac_step", PAC_STEPTYPE, IF_SET|IF_INTEGER, "PAC sweep step type (0 lin, 1 dec, 2 oct)" },
{ "pac_maxsb", PAC_MAXSB, IF_SET|IF_INTEGER, "PAC output conversion sidebands each side (0 = sideband 0 only)" }
{ "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)" }
};
SPICEanalysis PSSinfo = {

78
src/spicelib/parser/inp2dot.c
View File

@ -785,6 +785,80 @@ dot_pac(char *line, void *ckt, INPtables *tab, struct card *current,
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. */
static int
dot_pnoise(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 */
char *sname; /* the input source name */
CKTnode *nnode; /* a node pointer */
int which; /* which analysis we are performing */
char *steptype; /* pnoise sweep type: dec/oct/lin */
NG_IGNORE(gnode);
NG_IGNORE(current);
/* .pnoise Fguess StabTime OscNode Points Harmonics SC_iter Steady_coeff
* OutNode InSrc <DEC|OCT|LIN> NumPts Fstart Fstop */
which = ft_find_analysis("PSS");
if (which == -1) {
LITERR("Periodic noise (PNOISE) analysis unsupported.\n");
return (0);
}
IFC(newAnalysis, (ckt, which, "Periodic Noise 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, "pnoise_out", &ptemp));
INPgetTok(&line, &sname, 1); /* InSrc */
INPinsert(&sname, tab);
ptemp.uValue = sname;
GCA(INPapName, (ckt, which, foo, "pnoise_insrc", &ptemp));
/* sweep tail: <DEC|OCT|LIN> NumPts Fstart Fstop */
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));
ptemp.iValue = 1; /* enable the pnoise sweep */
GCA(INPapName, (ckt, which, foo, "pnoise", &ptemp));
return (0);
}
#endif
@ -986,6 +1060,10 @@ INP2dot(CKTcircuit *ckt, INPtables *tab, struct card *current, TSKtask *task, CK
} else if ((strcmp(token, ".pac") == 0)) {
rtn = dot_pac(line, ckt, tab, current, task, gnode, foo);
goto quit;
/* Enhancement-124: Periodic noise */
} else if ((strcmp(token, ".pnoise") == 0)) {
rtn = dot_pnoise(line, ckt, tab, current, task, gnode, foo);
goto quit;
#endif
#ifdef RFSPICE
}