Yule Walker Method — Control Systems/Spectral Measurements

Yule-Walker Method

Control Systems / Spectral Measurements

Estimates the power spectral density of a signal frame parametrically: it fits an order-p autoregressive model by the Yule-Walker (autocorrelation) method – exactly as the Yule-Walker AR Estimator block does – and evaluates the spectrum that model implies,

P_xx[k] = G·Ts / |A(e^j2πk/nfft)|², k = 0 … nfft−1,

where A is the prediction-error filter, G the model's error variance and Ts the time between samples of the series. This is MATLAB's pyulear with 'twosided' and a sampling frequency of 1/Ts.

It forms the frame's biased autocorrelation r[m] = (1/N)·Σx[n]x[n+m] at lags 0 to p and solves the Yule-Walker equations for it with the Levinson-Durbin recursion; the error variance is the recursion's final prediction error.

Ports

• u – the signal frame x, an [N,1] column; N must be greater than the order p.
• Pxx – the power spectral density, [nfft,1]: all nfft bins, two-sided, bin k at frequency k/(nfft·Ts). Units of power per hertz when Ts is in seconds.

Parameters

• Inherit Estimation Order – on or off (default off, as Simulink's). On, the order is p = N−1, the largest the frame allows, and Estimation Order is ignored; off, it is Estimation Order.
• Estimation Order – p, a whole number of 1 or more. Defaults to 6, as Simulink's does.
• FFT Length – nfft, the number of frequency bins, a whole number of 1 or more. Defaults to 256. It need not be a power of two. An nfft smaller than p+1 samples the spectrum of A wrapped modulo nfft rather than truncated – measured against the Simulink block, which does the same.
• Series Sample Time (s) – Ts, the time between samples of the series the frame was cut from, a positive scalar. Defaults to 1. It only scales the answer: halving Ts halves every bin.
• 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 frame length and the order are baked into the core at export time, and so are the FFT length, Ts and the nfft-point cosine and sine tables.

The three HDL targets are simulation-only: they carry the arithmetic in real and quantize only at the port boundary. The solve divides by a running error power that shrinks at every order, and the spectrum divides by |A|², which a Q16.16 datapath does not survive. The cores simulate correctly and are not offered as synthesizable.

Simulink bridge

Import and export, mapped to dspspect3/Yule-Walker Method – the DSP System Toolbox block. "Inherit Estimation Order" to inheritOrd (on/off, one to one), "Estimation Order" to ord, "FFT Length" to fftsize and "Series Sample Time (s)" to Ts.

inheritFFT = off and inheritTs = off are always emitted. With inheritFFT on the Simulink block picks its own length, and with inheritTs on it takes Ts from the frame period divided by N – measured: a 16-sample frame arriving every second scales every bin by exactly 1/16. This block runs one frame per step and has no frame period to divide, so Ts is always the dialog's value. Importing a model whose block inherits its Ts keeps the block's dialog Ts, so its spectrum can come out scaled differently from Simulink's, by exactly that ratio.

"Sampling Time (s)" does not cross. dspspect3/Yule-Walker Method defines no SampleTime parameter at all – verified against the R2026a block dialog – and set_param on a parameter a block does not define is a hard error that aborts the whole generated script.

Notes

• Stateless and algebraic: the whole solve runs on this sample's frame. Nothing is carried between steps.
• An all-zero frame is answered, not divided by. The recursion's division is guarded, so zeros give A = [1 0 … 0], K = 0 and G = 0 – measured on the Simulink block, and the same rule Autocorrelation LPC's zero-input handling applies. The other three estimators answer NaN there.
• Two-sided, every bin. Bins 1 … nfft−1 mirror about nfft/2 for a real frame; nothing is folded or doubled. For the one-sided density pyulear gives by default, double bins 1 to nfft/2−1.
• The three HDL targets' output port is Q16.16, whose largest value is 32767. The spectrum is G·Ts/|A|² and nothing keeps A's zeros off the unit circle – a sharp enough peak overflows that port on those three targets while the other seven carry it. Scale Ts down to keep every bin in range; it only scales the answer.
• Reach for the Yule-Walker AR Estimator block for the model itself rather than its spectrum.