Blog › ICP guides

MATLAB developer on retainer: Simulink, Signal Processing, and scientific computing on monthly retainer

September 25, 2026 · ~22 min read

An aerospace systems company had a flight control surface actuator model developed in Simulink over seven years. The model was the authoritative source for the embedded C code running on the actuator control unit, generated via Embedded Coder from a Simulink model and deployed after hardware-in-the-loop testing on an iron bird test rig. The team upgraded the development environment from MATLAB R2021b to R2023b during a tool qualification renewal. After regenerating the C code and running the existing HIL test suite, 14 of 400 test vectors failed — the generated code produced different numerical output for specific input sequences that exercised the DSP subsystem at saturation boundaries. The Simulink model had not changed. Only the MATLAB version had changed.

A MATLAB developer on Simulink retainer traced the problem in four hours. MathWorks changed the saturation behavior for certain fixed-point configurations in the Saturation block between R2022a and R2023b. The model's DSP subsystem used a model-wide Fixed-Point Toolbox data type override setting — Configuration Parameters → Simulation Target → Data Type Override: Double for simulation, with Scale to output type enabled for code generation. In R2023b, the interaction between the model-wide override and the block-level Saturate On Integer Overflow property changed for blocks whose output data type was derived from context rather than explicitly specified. Four Saturation blocks in the DSP subsystem had their output data type set to Inherit: Same as input rather than an explicit fixed-point type — in R2023b, these blocks received a different saturation semantic during code generation when the model-wide override was active. Fixing the issue required explicitly setting the output data type on each of the four blocks to the specific fixed-point type they were intended to use (fixdt(1, 16, 12) for signed 16-bit with 12 fractional bits), removing reliance on the model-wide override for these blocks. After regeneration, all 400 test vectors passed. The retainer work produced four block property changes in the Simulink model. A MATLAB developer on monthly retainer does this work continuously: auditing Simulink models against MathWorks release notes before tool upgrades introduce numerical behavior changes in generated code, profiling MATLAB M-file pipelines before vectorization regressions degrade runtime from 1.1 seconds to 8.4 seconds, and reviewing Control System Toolbox discretization method selections before integrator pole placement errors cause controller wind-up in production hardware.

MATLAB language fundamentals, numerical computing, and performance design

MATLAB's matrix-first data model treats all variables as arrays: a scalar is a 1×1 array; a row vector is 1×N; a column vector is N×1; a matrix is N×M. Operations applied to arrays are broadcast element-wise: A .* B is element-wise multiplication; A * B is matrix multiplication. Implicit expansion (R2016b+) allows operations between arrays of compatible sizes without explicit repmat: A + b where A is N×M and b is 1×M broadcasts b across all rows of A without allocating an intermediate matrix. Logical indexing provides concise in-place filtering: x(x < 0) = 0 sets all negative elements to zero; y = x(logical_mask) extracts elements where the mask is true. The colon operator generates sequences: 1:0.1:10 generates a row vector from 1 to 10 in steps of 0.1; A(:, 3) selects the third column of a matrix; A(end-2:end, :) selects the last three rows.

Cell arrays ({}) store heterogeneous data where each element can be any MATLAB type: c = {'string', 42, [1 2 3], struct('a', 1)}. Access uses curly braces: c{2} returns the number 42 as a double; c(2) returns a 1×1 cell array containing 42. Structs provide named field access with dot notation: s.name = 'signal'; s.fs = 8000; s.data = randn(1000,1). Struct arrays extend this to N structures with identical fields, vectorized over the struct array: signals(3).data accesses the data field of the third element. Function handles enable functional programming patterns: f = @(x) x.^2 + 2*x creates a lambda; arrayfun(f, [1 2 3 4]) applies it element-wise; cellfun(@numel, c) applies numel to each cell element. The nargin/nargout pattern enables functions with optional arguments: if nargin < 2, window = 'hann'; end provides a default for the second argument. inputParser provides named parameter parsing: p = inputParser; p.addRequired('signal'); p.addParameter('window', 'hann'); p.addParameter('overlap', 0.5); p.parse(sig, 'window', 'hamming') — p.Results.window returns the parsed value.

MEX files extend MATLAB with compiled C/C++ code for performance-critical operations that cannot be efficiently vectorized: the entry point is void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) where nlhs/nrhs are the number of left-hand/right-hand side arguments. mxGetDoubles(prhs[0]) returns a double* pointer to the first input array's data; mxGetM(prhs[0])/mxGetN(prhs[0]) return the number of rows and columns; mxCreateDoubleMatrix(m, n, mxREAL) allocates an output matrix; plhs[0] = outputMatrix sets the first output. MEX compilation: mex -R2018a myfunction.c compiles with the interleaved complex API (R2018a+). MATLAB's performance profiler — profile on; myFunction(data); profile viewer — shows cumulative time and call count for every M-file function and line, identifying hot spots that warrant MEX conversion or vectorization. The timeit function provides statistically stable timing by running the function multiple times and returning the median: timeit(@() myFunction(data)) avoids the JIT compilation artifact that affects the first tic/toc measurement.

Simulink model design, Embedded Coder configuration, and fixed-point arithmetic

Simulink models organize signal processing and control algorithms as block diagrams with explicit data flow. The model hierarchy uses subsystems to encapsulate functional modules — a Subsystem block contains internal blocks connected by signal lines; an Atomic Subsystem has its own execution context; a Model Reference references a separate .slx file, enabling incremental code generation where only modified referenced models are regenerated. Configuration set parameters determine simulation and code generation behavior: the Solver pane selects fixed-step vs variable-step integration — for embedded code generation, a fixed-step discrete solver with explicit step size matching the target hardware sample rate is required; variable-step solvers cannot generate embedded code. The Hardware Implementation pane specifies the target processor word sizes and byte ordering — correct settings are essential for generated integer arithmetic to match hardware behavior.

Simulink Coder generates portable C code; Embedded Coder generates optimized C/C++ for specific embedded targets with additional configuration for memory sections, function naming, and target-specific integer types. The code generation report — generated automatically with slbuild('model') and slreport — provides traceability between Simulink blocks and generated code lines, essential for certification workflows. Fixed-point design in Simulink uses the Fixed-Point Toolbox to specify data types explicitly: fixdt(1, 16, 12) specifies a signed (1) 16-bit word with 12 fractional bits (scaling factor 2^-12), representing values from -8 to +8 with resolution 2^-12 ≈ 0.000244. Block-level overflow behavior is set per-block in the block parameters dialog: Saturate On Integer Overflow checked means the block clamps output to the representable range on overflow; unchecked means wraparound overflow. The model-wide data type override (Configuration Parameters → Data Type Override) can override all fixed-point types to double for simulation, but the interaction between the override and the block-level saturation property changed in R2023b — blocks with inherited output types received different saturation semantics. The robust pattern is explicit data type specification on every Saturation block output rather than relying on type inheritance.

Stateflow adds hierarchical state machine design to Simulink. A Stateflow chart is a Simulink block containing states, transitions, and actions: entry: actions execute when entering a state; during: actions execute on every time step while in the state; exit: actions execute when leaving the state. Transitions between states have guards ([condition]), condition actions ({actions before transition}), and transition actions (/{actions after transition}). Hierarchical state machines use superstate containment — a superstate can contain multiple substates; a transition to the superstate boundary enters the default substate. Junctions allow complex branching logic within transitions: a history junction remembers the last active substate. MATLAB Function blocks embed M-code directly in the Simulink diagram with explicit input/output port declarations via function out = computeSignal(in1, in2) — the code is analyzed by Embedded Coder for code generation compatibility, requiring fixed-type variables (in1 = cast(in1, 'double')) and no dynamic allocation. The Simulink Model Advisor checks model configuration against best practices for simulation accuracy, code generation correctness, and certification guidelines — running it before each release migration identifies configuration warnings before they become generated code failures.

Signal Processing Toolbox, Control System Toolbox, and algorithm design

The Signal Processing Toolbox provides the core functions for spectral analysis, filter design, and digital signal processing algorithm implementation. Fast Fourier Transform: X = fft(x, N) computes the N-point DFT; the output has N complex values, with the first N/2+1 values representing frequencies from DC to the Nyquist frequency. For spectral analysis, windowing the input before the FFT reduces spectral leakage: X = fft(x .* hann(length(x))', N) applies a Hann window. The spectrogram function computes the short-time Fourier transform: [S, F, T, P] = spectrogram(x, window, noverlap, nfft, fs) returns the complex spectrogram, frequency vector, time vector, and power spectral density. Filter design with designfilt provides a specification-based API: d = designfilt('lowpassiir', 'FilterOrder', 8, 'HalfPowerFrequency', 0.4, 'DesignMethod', 'butter') creates an 8th-order Butterworth lowpass filter; filter(d, x) applies it; fvtool(d) visualizes the frequency response. For high-order IIR filters, direct-form implementation accumulates numerical errors — sosfilt(d.Coefficients, x) applies the filter as a cascade of second-order sections, which is numerically stable for high-order designs where the direct [b,a] form would exhibit coefficient quantization errors.

The Control System Toolbox models dynamic systems as transfer functions, state-space representations, or zero-pole-gain models. G = tf([1], [1, 2, 1]) creates the transfer function 1/(s^2 + 2s + 1); G = ss(A, B, C, D) creates a state-space model with system matrices; G = zpk(zeros, poles, gain) creates a zero-pole-gain model. Analysis functions: bode(G) plots the Bode diagram (magnitude and phase vs frequency); step(G) plots the step response; pzmap(G) plots poles and zeros in the s-plane; margin(G) returns gain margin and phase margin for stability analysis. pidtune(G, 'PID') automatically synthesizes a PID controller for the plant G using internal loop shaping; the result is a pid controller object that can be analyzed with bode, step, or passed to feedback for closed-loop analysis. Discretization of continuous-time controllers for implementation on digital hardware uses c2d: Cd = c2d(C, Ts, 'tustin') converts the continuous controller C to a discrete-time controller at sample period Ts using the bilinear (Tustin) transformation, which preserves stability by mapping the left-half s-plane to the interior of the unit disk in the z-plane. The Euler forward method ('foh' or 'forward euler') maps the integrator pole at s=0 to z=1+Ts — outside the unit circle for any positive Ts — making the discrete integrator marginally unstable; Tustin maps the same pole to exactly z=1, preserving the continuous-time marginally stable integrator semantics.

Retainer work auditing Control System Toolbox designs focuses on discretization method correctness for the intended sample rate and controller structure. A PID controller with integral action discretized at a sample period that is large relative to the integrator time constant requires Tustin or ZOH discretization — not Euler forward, which introduces integrator instability at any sample period. Anti-windup design — limiting the integrator state when the actuator saturates — prevents the classical integrator wind-up problem where the controller drives the integrator to large values during saturation, causing overshoot when the saturation constraint releases. In the discrete domain, anti-windup is implemented by clamping the integrator state update: if u > u_max; I = I - Ki * e * Ts; end — the integrator increment is reversed when the saturated output indicates the actuator has already reached its limit, preventing further accumulation. Comparing the continuous-time controller design against the discrete implementation using bode(C) and bode(Cd) on the same figure verifies that the frequency response matches at frequencies well below the Nyquist rate before deploying to hardware.

Parallel Computing Toolbox and parfor optimization

MATLAB's Parallel Computing Toolbox enables multi-core and cluster parallelism through the parfor loop, spmd blocks, and GPU array operations. parfor replaces for for loops where iterations are independent — the MATLAB worker pool executes iterations across available workers: parfor i = 1:N; result(i) = computeSignal(data(:,i), params); end. The MATLAB code analyzer enforces parfor constraints at edit time: loop variables cannot have inter-iteration dependencies; loop body cannot use break, return, or continue; arrays indexed by the loop variable are sliced (each worker receives only its slice, not the full array); arrays read but not written are broadcast (sent to all workers in full at loop start). The broadcast variable cost is significant for large arrays — a 100MB parameter matrix broadcast to 8 workers requires 800MB of inter-process data transfer before the first iteration executes. The fix: extract invariant subsets of large parameters outside the loop, or convert the broadcast array to a parallel.pool.Constant which is sent to workers once at creation time and cached across multiple parfor loops: params_c = parallel.pool.Constant(params); parfor i = 1:N; result(i) = computeSignal(data(:,i), params_c.Value); end.

GPU arrays accelerate matrix-intensive computation by offloading arithmetic to NVIDIA GPU cores. A_gpu = gpuArray(A) transfers a MATLAB array to GPU memory; most MATLAB built-in functions operate transparently on gpuArrays without code changes — fft(A_gpu) runs cuFFT on the GPU; A_gpu * B_gpu uses cuBLAS GEMM. gather(result_gpu) transfers the result back to host memory. GPU acceleration is most effective for large matrix operations (GEMM, FFT) where the arithmetic intensity — flops per byte transferred — is high enough to amortize the PCIe transfer cost; for small matrices or algorithms with low arithmetic intensity, GPU overhead exceeds the computation cost. The gputimeit function provides GPU-accurate timing by synchronizing the GPU before and after measurement, avoiding the asynchronous execution issue where tic/toc measures only the time to submit work to the GPU, not the time for the GPU to complete it. Tall arrays extend MATLAB's normal array operations to datasets too large to fit in memory: t = tall(datastore('*.csv')) creates a tall array backed by a datastore; mean(t), sum(t)`, and most element-wise operations work on tall arrays via deferred evaluation — the computation executes when gather(result) is called, processing data in chunks that fit in memory.