Multitaper PSD — Control Systems/Spectral Measurements

Control_Systems/Spectral_Measurements/Multitaper_PSD · 1 input / 1 output port(s) at insert · exports to Python, MATLAB, Java, Rust, C, C++, VHDL, Verilog, SystemVerilog, PLC Structured Text

Description#

The block's own DESCRIPTION_HTML, rendered verbatim — the same text the config dialog's info panel and the library navigator show. Fix a wrong sentence in the block's .cpp (R-D9), never here.

Multitaper PSD

Control Systems / Spectral Measurements

The one-sided power spectral density of the signal by Thomson's multitaper method – MATLAB's pmtm – over a running window of the last N samples. The window is multiplied by each of K orthogonal Slepian tapers e_k (the discrete prolate spheroidal sequences of time-halfbandwidth product NW), each product's periodogram is taken, and the K eigenspectra are combined:

S_k(f) = |Σ_n e_k[n]·x[n]·e^−j2πfn/N|², P(f) = c(f)·S(f) ÷ fs, with c = 1 at DC and Nyquist and 2 between.

Where Welch PSD lowers the variance by averaging segments, this block averages tapers over the SAME samples, so a short record keeps its full resolution. It is pmtm(x, E, V, N, fs, method) with the tapers of this window, and the defaults are pmtm's: NW = 4 and 2·NW − 1 = 7 tapers, weighted adaptively.

Ports

u – the sampled signal. Scalar: one channel and its own window.
p – the density, [N/2+1, 1], one-sided with DC first and Nyquist last. This is the vector every other block in this family takes on its own p port, so it wires straight into Band Power, Occupied Bandwidth, Spectral Entropy and the rest.

Parameters

Record Length – N, the window and transform length. An even whole number from 4 to 64; even because a one-sided spectrum needs a Nyquist bin. There is no zero-padding: pmtm's own default FFT length, max(256, 2^nextpow2(N)), would interpolate the same spectrum onto 256 points.
Time-Halfbandwidth Product – NW, the tapers' concentration: each estimate is smoothed over a band of 2·NW bins. At least 0.75 (pmtm's own floor for tapers it is given) and below N/2. Defaults to 4.
Number of Tapers – K, how many of the most concentrated Slepian sequences are used, a whole number from 2 to N. pmtm's default is round(2·NW) − 1 – it designs 2·NW and drops the last – and this block's default of 7 is exactly that at NW = 4. Tapers beyond 2·NW leak, which is why pmtm stops short of them.
Weighting – how the eigenspectra are combined:
- Adaptive – pmtm's default. Weights b_k²λ_k with b_k = S ÷ (S·λ_k + σ²(1 − λ_k)), iterated from (S₁+S₂)/2 until the change falls below pmtm's own tolerance. λ_k is each taper's concentration and σ² the window's mean square.
- Eigenvalue – Σλ_kS_k ÷ K. Divided by K, not by Σλ_k, exactly as pmtm does.
- Unity – ΣS_k ÷ K, the plain average.
Bin Spacing – Δf, the frequency step between output bins, in hertz. The sampling rate follows as fs = N·Δf and is what makes the answer a density rather than a power. It is the same parameter, under the same name, that Welch PSD, Band Power and the rest of this family read.
Sampling Time (s) – zero or less inherits the solver's rate; a positive value runs the block at that period.

Code export

All ten targets: Python, MATLAB, Java, Rust, C, C++, VHDL, Verilog, SystemVerilog and PLC Structured Text. The tapers, their concentrations, the transform's twiddles and the one-sided scale are computed once at export time and carried as numbers; the eigenspectra and the weighting run as loops in the core, which is what the adaptive weighting needs: its weights depend on the data and it iterates to a tolerance. All five settings are structural, so re-export after changing any of them.

⚠ The three HDL targets run that loop in simulation-only real arithmetic, quantizing only at the two port boundaries: the adaptive weighting divides by the data, which a Q16.16 datapath does not carry. They are not offered as synthesizable, whatever the weighting.

The adaptive loop stops on pmtm's own rule, and additionally after 500 iterations, which pmtm does not have; the most measured was 54.

Simulink bridge

None (Support::None). pmtm is a Signal Processing Toolbox function and that toolbox ships no Simulink library. DSP System Toolbox's dspspect3/Spectrum Estimator offers only the Welch and filter-bank methods, and no library block has a multitaper method. The bridge reports this block rather than dropping it silently, and it therefore has no parity testbench; code export verification still covers it across all ten languages.

Notes

Stateful, and discrete by nature (setDiscreteOnlyBlock(true)): one shift register of N samples, zero at the start of the run, so the first N−1 outputs are the spectrum of a record that is partly zeros.
⚠ pmtm drops the last taper even when the tapers are handed to it. pmtm(x, E, V) with K tapers uses K−1 unless 'DropLastTaper' is false; this block uses all K it is asked for. To compare with a pmtm call that passes an NW, set K to round(2·NW) − 1.
An all-zero window gives an all-zero density under every weighting – the adaptive loop's tolerance is zero and it does not start.

Code facts#

Fact	Value
registered type	`Control_Systems/Spectral_Measurements/Multitaper_PSD`
family	`Control_Systems/Spectral_Measurements`
solver environment class	`ICoreBlock_0_Control_Systems_1_Spectral_Measurements_2_Multitaper_PSD`
source	`src/ICoreBlocks/ICoreBlockLibrary/Blocks/Control_Systems/Spectral_Measurements/Multitaper_PSD/ICoreBlock_0_Control_Systems_1_Spectral_Measurements_2_Multitaper_PSD.cpp`
header	`src/ICoreBlocks/ICoreBlockLibrary/Blocks/Control_Systems/Spectral_Measurements/Multitaper_PSD/ICoreBlock_0_Control_Systems_1_Spectral_Measurements_2_Multitaper_PSD.h`
default size on canvas	140 × 76 px
ports at insert	1 in, 1 out
code generators implemented	`Python`, `MATLAB`, `Java`, `Rust`, `C`, `C++`, `VHDL`, `Verilog`, `SystemVerilog`, `PLC Structured Text`

Ports#

#	Direction	Signal type	Description label
1	in	`ICoreDouble`	u
2	out	`ICoreDouble`	p

Ports the constructor creates. A block whose port list changes with its configuration adds or removes ports at load time; the count above is the one a freshly inserted block has.

Configuration variables#

Config variable	Default	Simulink parameter
`Record Length`	`32`	not crossed
`Time-Halfbandwidth Product`	`4`	not crossed
`Number of Tapers`	`7`	not crossed
`Weighting`	`Adaptive%~%Eigenvalue%~%Unity~~Adaptive`	not crossed
`Bin Spacing`	`1`	not crossed

Every block also carries Sampling Time (s) from ICoreBlockSolverEnvironment: zero or less inherits the solver's rate, a positive value runs the block at that period.

Simulink bridge#


support	`Support::None`
Simulink path	—
port-count rule	`PortsParam::None`
`SampleTime` parameter	yes
deliberately not crossed	`Record Length`, `Time-Halfbandwidth Product`, `Number of Tapers`, `Weighting`, `Bin Spacing`

Caveat (shown to the user): pmtm is a Signal Processing Toolbox function and that toolbox ships no Simulink library. DSP System Toolbox's dspspect3/Spectrum Estimator offers only the Welch and filter-bank methods (its Method parameter, measured on R2026a), and a search of 119 block libraries found no multitaper, Thomson or Slepian block, so there is no path a diagram could name

Catalog contract: src/ICoreBlocks/ICoreCoder/ICoreCommandSystem/SimulinkBridge/ICoreSimulinkBlockCatalog.h

Description vs code#

The lists agree. check_block_descriptions.py finds no disagreement between the description's Ports, Parameters, Code export and Simulink bridge lists and the code's.

The verdict above is tools/docs/check_block_descriptions.py (P7.1), which compares LISTS. It cannot read a sentence: "stateless" on a block with a state, an initial-value semantic the recursion does not implement, a "not synthesizable" caveat the HDL banner contradicts. That is the agent audit (P7.3) on BLOCK_DESCRIPTION_AUDIT.md, and this tool's green is not a substitute for one.

The top comment of the block's .cpp — the maths, the realization and the export strategy, addressed to whoever changes it. It must not contradict the description above (P7.5).

Multitaper PSD -- the one-sided power spectral density of a stream, MATLAB's pmtm Over the last N samples x[0..N-1] (oldest first), with K Slepian tapers e_k of time-halfbandwidth product NW and their concentrations lambda_k:

S_k(f) = | SUM_n e_k[n] x[n] exp(-j 2 pi f n / N) |^2 the K eigenspectra S(f) = SUM_k w_k(f) S_k(f) / SUM_k w_k(f) adaptive = SUM_k lambda_k S_k(f) / K eigenvalue = SUM_k S_k(f) / K unity P(f) = c(f) S(f) / fs, c = 1 at DC and Nyquist, 2 between one-sided

⚠ EVERY CONVENTION HERE WAS READ OUT OF R2026a's pmtm.m AND dpss.m AND THEN MEASURED, by a prototype of exactly this formulation against pmtm over N = 8..32, NW = 1.5..4, all three weightings: the largest disagreement on any bin is 8e-15 relative. What that pinned:

THE EIGENVALUE WEIGHTING DIVIDES BY K, NOT BY SUM(lambda): pmtm writes Sk*V/k.
pmtm DROPS THE LAST TAPER EVEN WHEN THE TAPERS ARE HANDED TO IT: pmtm(x, E, V) with K

tapers uses K-1 unless 'DropLastTaper' is false. Its default, pmtm(x, NW), designs round(2 NW) tapers and drops one, so this block's "Number of Tapers" default of 7 at NW = 4 IS pmtm's default -- measured identical to the explicit form at relative 0.

THE ADAPTIVE LOOP stops when SUM |S - S_prev| / nfft over ALL nfft bins, both halves,

falls to 0.0005 * (x'x/N) / nfft, starting from S = (S_1 + S_2)/2 and S_prev = 0, with weights b_k^2 lambda_k, b_k = S / (S lambda_k + (x'x/N)(1 - lambda_k)). The block runs the test over both halves in MATLAB's order, reading the negative-frequency bins off their positive twins, so it stops on the same iteration pmtm does. pmtm has no iteration limit; this block stops at 500, and the most measured was 54 (3000 windows per configuration, including start-up windows that are mostly zeros).

THE CONCENTRATION lambda_k is dpss's own: the taper's autocorrelation r against

[2W, 4W sinc(2Wm)], clamped to [0, 1]. The tapers are the eigenvectors of the same tridiagonal matrix dpss uses, and agree with dpss's to 4.5e-15 (signs chosen as dpss chooses them, although a spectrum cannot see a sign).

The FFT length is N, so the answer is pmtm(x, E, V, N, fs, method) with the tapers of this window -- pmtm's own default FFT length, max(256, 2^nextpow2(N)), would zero-pad to 256.

⚠ THE ADAPTIVE WEIGHTS DEPEND ON THE DATA AND ITERATE TO A TOLERANCE, which no fixed set of dot products expresses, so every target runs the same statement list: the tapers, the twiddles and the scale are computed once at export time and carried as numbers, and the loops run in the core. The three HDL targets run it in simulation-only real arithmetic.

Sample results#

The same rig also ran:

Stimulus	What it is	Output range
`impulse`	Impulse: one sample of 1 at k = 5, 0 elsewhere (Repeating Sequence Stair)	0 … 0.001091
`ramp`	Ramp: slope 1 from t = 0	0 … 2.727
`sine`	Sine Wave: amplitude 1, 2 rad/s, no phase, no bias	0 … 0.06968
`table`	Repeating Sequence Stair: [-2 -1 -0.5 0 0.5 1 2 3], one entry per sample	3.86e-8 … 0.08063

Plotted: step — Step: 0 -> 1 at t = 1 s

Category dynamic · sample time 0.1 · 60 steps · commit 93133d604 · produced by docsSample --out <folder> --blocks Gain_Scheduled_Lead_Lag Controller_1D Controller_Blend_1D Controller_2D Controller_3D Observer_Form_1D Self_Conditioned_1D Line_Of_Sight_Access Orbit_Propagator_Kepler Attitude_Dynamics Attitude_Profile_Nadir_Pointing Attitude_Profile_Geographic_Pointing Multitaper_PSD Cross_Power_Spectral_Density Transfer_Function_Estimate Envelope_Spectrum Compose_String Scan_String --steps 60 · data docs/generated/samples/Control_Systems__Spectral_Measurements__Multitaper_PSD.json · the SVG is generated from those numbers by tools/docs/plot_svg.py, so it is a run and not a drawing (R-D10).