Fourier Fit — Control Systems/Curve Fitting

Fourier Fit

Control Systems / Curve Fitting

Fits a truncated Fourier series to the last W samples and emits its coefficients:

ŷ(x) = a₀ + Σ_m=1..n [ a_m·cos(m·ω·x) + b_m·sin(m·ω·x) ]

This is MATLAB's fourier<n> model, and the output carries its coefficients in that model's own order – a₀, a₁, b₁, …, a_n, b_n. The newest sample sits at x = 0 and the sample i steps back at x = −i·h, the same abscissa Polynomial Fit uses.

Ports

• u – the sampled signal being fitted. Scalar: one channel and its own window – see Notes.
• c – the coefficient column, [2n+1, 1], in the order above. Its height follows Number of Harmonics alone.

Parameters

• Window Length – W, how many samples the fit sees. A whole number from 3 to 101, and it must be at least 2n+1: that many coefficients cannot be determined by fewer points. Bounded above because the dot products are unrolled at export.
• Number of Harmonics – n, a whole number from 1 to 8, matching MATLAB's fourier1…fourier8. Each harmonic adds a cosine and a sine, so the output grows by two.
• Fundamental (rad per unit) – ω, the angular frequency of the first harmonic, in radians per unit of x. A single positive number. ⚠ This is a parameter here and an answer in MATLAB – see Notes, because it is the difference between a linear fit and a nonlinear one.
• Abscissa Step – h, the spacing between consecutive samples on the fitted axis. A single positive number. Only the product ω·h reaches the arithmetic, so setting h to the sampling period and ω in radians per second is the same fit as leaving h at 1 and giving ω in radians per sample.
• 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 (2n+1)×W coefficient bank is structural and is inlined into the arithmetic at export time rather than exposed as a tunable parameter: it follows from all four settings, and changing any of them changes how many multiplies and how many outputs the core contains. Re-export after changing them. No sine or cosine is evaluated in any generated core – they were evaluated once, at export time, to build the bank.

The three HDL targets are genuine synthesizable Q16.16: a shift register and one fixed multiply-accumulate per coefficient, with the products accumulated at full width and shifted back once per row.

Simulink bridge

None (Support::None). Fitting a Fourier series to a running window is a Curve Fitting Toolbox function (fit with a fourier<n> model), and that toolbox ships no Simulink library at all, so there is no path a diagram could name. 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)): the window advances once per sample.
• ⚠ Its output is NOT a polynomial, and Evaluate Fit cannot read it. Evaluate Fit, Fit Derivative and Fit Integral all speak coefficients in descending powers; these are amplitudes against a basis of cosines and sines. Wiring this block into any of the three produces a number rather than an error, so the mistake is silent – reconstruct the series with Trigonometry blocks and a Sum instead.
• ⚠ The fundamental is a parameter here and MATLAB solves for it. That is the one real departure from fit(x, y, 'fourier2'), and it is what makes the block possible: with ω unknown the fit is nonlinear and needs an iterative solver every sample; with ω given it is an ordinary linear least squares whose design matrix is fixed by the configuration. If the frequency is what you are looking for, this is the wrong block – reach for Machine Learning / Feature Engineering / FFT Magnitude or Transforms / Goertzel first, then set the result here.
• One solve, then arithmetic. The small normal-equation system is solved once when the configuration loads; every sample afterwards is 2n+1 fixed dot products.
• ⚠ Not every setting is a fit, and the block refuses rather than guessing. The basis degenerates on a sample grid: ω·h at an exact multiple of π makes a sine column vanish identically, and a harmonic above the grid's Nyquist rate is indistinguishable from a lower one. Either makes the normal matrix singular, which is detected when the configuration loads and reported with a reason – rather than emitting the arbitrary answer a pseudo-inverse would have produced.
• Verified against MATLAB, at two conditionings, because the difference is the caveat. Both compare against R2026a's own least-squares solve of the same design matrix (a QR solve) over a 13-sample window. At W = 13, n = 3, h = 0.3, ω = 1.7 – where ω·h = 0.51 rad per sample and cond(A'A) = 2.21 – this block agrees to 4.3e−16, which is the last bit. Move to a poorly spread basis (W = 11, n = 2, h = 0.3, ω = 0.8, cond(A'A) = 3103) and the agreement loosens to 2.0e−13. Nothing is wrong in either case: solving the normal equations squares the condition number where a QR solve does not, and the normal equations are what let the whole fit collapse into one inlined bank of numbers instead of putting a factorization in ten generated cores. Keep n·ω·h well under π and the loss never appears.
• The window is zero-prefilled, and the zeros count. The first W−1 coefficient vectors of a run are a startup transient, matching Moving Median, Detrend, Savitzky-Golay Filter and Polynomial Fit.
• Scalar only. One channel and its own history; wire one block per channel.
• No state space. Linear in the input, but through a fixed filter bank over W past samples rather than an A/B/C/D pair, so model reduction correctly declines to merge it.