Optimising Time-Frequency Distributions: A Surface Metrology Approach

Time-frequency signal processing offers a significant advantage over temporal or frequency-only methods, but representations require optimisation for a given signal. Standard practice includes choosing the appropriate time-frequency distribution and fine-tuning its parameters, usually via visual inspection and various measures—the most commonly used ones are based on the Rényi entropies or energy concentration by Stanković. However, a discrepancy between the observed representation quality and reported numerical value may arise when the filter kernel has greater adaptability. Herein, a performance measure derived from the Abbot–Firestone curve similar to the volume parameters in surface metrology is proposed as the objective function to be minimised by the proposed minimalistic differential evolution variant that is parameter-free and uses a population of five members. Tests were conducted on two synthetic signals of different frequency modulations and one real-life signal. The multiform tiltable exponential kernel was optimised according to the Rényi entropy, Stanković’s energy concentration and the proposed measure. The resulting distributions were mutually evaluated using the same measures and visual inspection. The optimiser demonstrated a reliable convergence for all considered measures and signals, while the proposed measure showed consistent alignment of reported numerical values and visual assessments.


Introduction
Real-life signals, whether natural or artificial, frequently require joint time-frequency analysis to better understand the internal signal structure due to such signals' amplitude and frequency temporal variability. Time-frequency representations result from various time-frequency distributions and are methods of describing temporal variations of signals' spectral content [1]. In general, there are linear, quadratic and higher-order classes of distributions with applications investigated in diverse fields of research, such as instantaneous frequency estimation [2][3][4][5][6], parameter determination [7][8][9], signal modelling [10][11][12], wireless resource management [13,14], target localisation [15][16][17], and biomedical signal abnormalities [18][19][20][21][22][23]. Time-frequency representations augment and extend the classical analysis by offering additional features that can help better discriminate and classify diverse phenomena, and this is paramount for emerging pattern recognition and machine learning applications [24]. However, the fundamental challenge in successfully applying time-frequency analysis and processing methods lies in the appropriate choice of the distribution and its filtering kernel parameters.
According to the literature, the quadratic class of distributions is of the highest interest since it has a signal composition displayed as distinct concentrations of energy along the time-frequency (t, f ) plane that resemble a landscape of ridges of varying heights. Crests of these ridges correspond to signal components' instantaneous frequency (IF) and directly suitable as a performance measure due to its noise sensitivity, a threshold variant of the 0 norm was successfully applied in automatic SPWVD optimisation in [36].
Besides the concentration, TFD performance is commonly measured as the distribution complexity using generalised Rényi entropy with different normalisations applied in practice [1,37]. In [38], Sucic et al. based their instantaneous number of signal components estimation algorithm on the short-term Rényi entropy. Saulig et al. later used it in spectrogram optimisation [39]. To define performance differences and selection criteria between various TFDs, Boashash and Sucic introduced the normalised instantaneous resolution (NIR) measure based on per-time slice measurements of classical spectral performance indicators [40]: components' instantaneous bandwidths, side-lobe to main-lobe and crossterm to main-lobe amplitude ratios of two spectrally closest peaks. Although in-depth, the method proved too complicated for TFD optimisation due to its closest peak-finding sensitivity and because it requires an already well-performing exploratory distribution to identify the cross-terms [25]. The Reinhold measure, which was subsequently proposed, attempts to overcome NIR issues but introduces two additional user-supplied parameters while limiting the analysis to only two components at any moment [41].
The optimality of a time-frequency distribution is usually defined relative to the ideal representation, whereas the WVD has been found to only be optimal for a mono-component linear frequency modulated (LFM) signal [1]. Exploiting this optimality, Al-Sa'd et al. [25] proposed decomposing the signal into LFM segments from which a piece-wise spline WVD is constructed and used as a reference for comparing the accuracy and resolution of other TFDs. Based on the proposed performance evaluation process, popular distributions such as CKD, SPWVD and EMBD have again been recommended for the general case of multi-component non-stationary and non-linear instantaneous frequency laws signals.
Rényi entropy and Stanković's ECM are among the more commonly used measures [24,25]. As a global performance indicator in the form of a single value, such measures are convenient for implementing TFD optimisation algorithms. When stated as the general optimisation problem with only the analysed signal and TFD kernel parameters as input arguments, such algorithms allow simple modular inclusion into machine learning frameworks. However, a quantitatively optimised filtering kernel may distort some or all of the information in the resulting TFD, making a visual evaluation by a signal specialist indispensable. In addition to numerically capturing TFD performance, a credible metric should also align with how the TFD is visually perceived to perform.
In this study, our objective was to define a global measure relative to the uniform time-frequency distribution and agreeable with the visual assessment using definitions from the surface metrology as inspiration for the TFD performance indicators.
A product's texture and surface properties are often critical to functional specifications, and material sciences and manufacturing have long recognised the need for controlling and appropriately measuring those properties. Many parameters to capture and quantify surface complexities of processed materials were previously used in national standards, which resulted in what was known as the "parameter rash" [42]. In order to address such a situation, the regulating bodies under the SURFStand project [43] put forward a set of parameters based on statistical, functional and morphological surface features that ultimately became accepted and are continuously developed as the ISO-25178 international standard for geometric product specification [44]. Of particular interest, we found definitions of aerial surface texture parameters that served as the basis for the proposed measure.

Materials and Methods
Essentially, a quadratic (t, f ) distribution, denoted ρ z (t, f ), is a height map and can be viewed as a surface texture that results from (t, f ) smoothing of the rough WVD in the process defined as [1]: where γ(t, f ) is the filter kernel convolved ( * * ) with the Wigner-Ville distribution W z (t, f ) of the complex signal z(t). The WVD, in turn, results from the Fourier transform (FT) over the time-lag variable τ of the signal instantaneous auto-correlation function (IAF), denoted K z (t, τ), as in Equation (2): where F is the forward FT, and z(t) is the analytic associate of the real signal s(t). However, QTFDs are usually designed in the Doppler-lag (ν, τ) domain and then transformed back to the (t, f ) domain as in Equation (3): where F −1 is the inverse FT, g(ν, τ) is the kernel in the Doppler-lag form and A z (ν, τ) is the signal ambiguity function (AF) constructed from the IAF, or, equivalently, from the WVD as in Equation (4) [1]: Aside from efficient filtering, the Doppler-lag domain makes kernel design easier since it provides some separation between the cross-terms and auto-terms. Filter kernels in the (ν, τ) form are commonly the low-pass type that preserves the signal AF around its origin and de-emphasises locations further away, usually occupied by cross-terms. No smoothing is performed if the kernel is the all-pass type, i.e., g(ν, τ) = 1, and the result is the WVD. Conversely, if the kernel is zero everywhere but the signal AF origin, it will flatten the resulting TFD by distributing the signal energy uniformly across the (t, f ) plane. Compared to the uniform, any other TFD of the same signal is more informative and can be considered better-performing, so we use it as the reference for the proposed metrics in this study.
In order to establish the metrics, we follow the approach from the surface metrology by first constructing the Abbot-Firestone curve of a discrete TFD. The Abbot curve is more frequently known as the bearing area or material ratio curve (MRC) because of its usage in definitions of functional surface parameters [45]. Essentially, the Abbot curve is a form of an empirical cumulative function that shows the ratio of a surface cross-sectional area at a considered height to the evaluation area.
The material ratio curve of a discrete (t, f ) distribution ρ z [n, k] of size N · M, with n and k being discrete-time and frequency index, respectively, can be formed as the finite sequence of non-increasing values: obtained from ρ z [n, k] as the centred samples [45] where m [n,k] is the average of ρ z [n, k] defined as [1] Figure 1 shows the resulting MRCs of the uniform and the WVD of an example signal consisting of two crossed LFM components. The horizontal axis represents the areal material ratio mr, while the vertical axis is the height h of the surface cross-section at the considered material ratio relative to the reference surface, which is the uniform TFD.  (25) and (26), respectively. The shaded area is the material volume Vm(mr(c 2 )), while the smaller, darker part is the material volume Vm(mr(c 1 )) corresponding to the significant auto-terms.
For a considered level c and the evaluation area equalling the size of ρ z [n, k], we 174 calculate mr(c) as in Equation (8) 175 which is essentially a threshold form of the ℓ 0 norm. The area between an MRC and 176 the intersecting line at c or at mr(c) represents a volume. For a given mr, similarly to the 177 volume parameters from the surface metrology [45], we define the area under the positive 178 section of an MRC as the material volume Vm(mr) which corresponds to positive energy 179 in ρ z [n, k]. As ρ z [n, k] can be considered a densely sampled surface, a straightforward 180 summation can be used to calculate the volume Vm(mr) as in Equation (9) 181 As evident in Figure 1, Vm(·) of the WVD at mr(0) has an equal-area counterpart 182 between level zero and the negative of the MRC since h i samples are centred. The shape of 183 this area and its deepest recess can be used to estimate the cross-term contribution to the 184 material volume due to their oscillatory nature. By comparing the shapes of the positive 185 and negative areas, it is clear how cross-terms make a substantial part of the material 186 The horizontal axis is the areal material ratio mr, and the vertical axis is the height h of the surface cross-section at the considered material ratio relative to the reference surface, i.e., the uniform TFD. The curves of the uniform (flat line at c 2 = 0) and Wigner-Ville distribution (S-shaped curve) of the example signal defined in Equation (24) consisting of two crossed LFM components as per Equations (25) and (26), respectively, are shown. The shaded area is the material volume Vm(mr(c 2 )), while the smaller, darker part is the material volume Vm(mr(c 1 )) corresponding to the significant auto-terms.
For a considered level c and the evaluation area equalling the size of ρ z [n, k], we calculate mr(c) as in Equation (8): which is essentially a threshold form of the 0 norm. The area between an MRC and the intersecting line at c or at mr(c) represents a volume. For a given mr, similar to the volume parameters from the surface metrology [45], we define the area under the positive section of an MRC as the material volume Vm(mr), which corresponds to positive energy in ρ z [n, k].
As ρ z [n, k] can be considered a densely sampled surface, a straightforward summation can be used to calculate the volume Vm(mr) as in Equation (9): As is evident in Figure 1, Vm(·) of the WVD at mr(0) has an equal-area counterpart between level zero and the negative of the MRC since h i samples are centred. The shape of this area and its deepest recess can be used to estimate the cross-term contribution to the material volume due to their oscillatory nature. By comparing the shapes of the positive and negative areas, it is clear how cross-terms make a substantial part of the material volume of the WVD and increase mr(0). In contrast, the uniform TFD has no material volume, but it takes the entire evaluation area N · M since its MRC is the line at level zero. These observations demonstrate that the MRC of the optimal TFD should be nearly flat in its negative part with the material ratio mr(0) smaller than that of the WVD.
Based on the above insight, the relevant indicators are then the material volume and ratios at levels c 1 = |min(h i )| and c 2 = 0, respectively, from which the proposed performance metrics follow as with mr(·) and Vm(·) calculated according to Equations (8) and (9), respectively. V + is the positive volume of ρ z [n, k] calculated as Since p 1 can be close to zero for nearly flat distributions or the WVD, p 3 is needed to discriminate higher concentrations further. However, p 1 and p 3 can promote cross-terms, so p 2 provides the counterbalance by only considering values greater or equal to c 1 , i.e., the significant auto-terms, as the principal concentration contributors. For the uniform TFD, p 1 , p 2 and p 3 all equal one, so it is uniquely ranked as the worst-performing TFD. On the other extreme, in the ideal TFD case, the p values are minimal, and it holds that p 2 = p 3 .
The so-defined p metrics evaluate the performance of ρ z [n, k] based on its values, but, in general, ρ z [n, k] may not strictly uphold the energy conservation property and, therefore, may not correspond well with the actual signal. Round-off errors, non-compact support kernels and aliasing in discrete realisations of (t, f ) distributions can violate this property. However, most of the signal's energy should be preserved within the evaluation area of size N · M. Otherwise, it can be argued that a (t, f ) distribution no longer represents the signal's energy distribution. In order to measure the relative amount of this violation, we introduced a discrepancy metric between the average across the whole matrix of ρ z [n, k] and the average power of the analysed signal z[n] as Smaller p values indicate a better-performing TFD, but all p metrics must be considered simultaneously and equally important, thus presenting a multi-objective optimisation problem with equally weighted individual objectives. In this sense, we formulated the proposed measure as that is, the average of the proposed p metrics. As a single-valued global indicator, the proposed measure PM can directly serve as the objective function in a TFD optimisation algorithm. The problem of selecting optimal kernel parameters in such an algorithm can generally be stated as arg min kernel parameters that is, find the kernel parameters for which some performance measure, such as Equation (15), (22) or (23), attains a minimal value globally. Within such a framework, there are three parts to an algorithm then: the kernel, the performance measure, and the core optimiser. As a solution to the posed problem, kernel parameters should be globally optimal; therefore, the core optimiser must also be of the global type. From the plethora of algorithms for global optimisation, we chose the differential evolution (DE) because of its proven effectiveness, algorithmic simplicity and the small number of control parameters. Differential evolution is an instance of the evolutionary class of global optimisers in evolutionary computing. It is a stochastic direct search, a derivative-free algorithm similar to other popular methods such as the genetic algorithm [46] and particle swarm optimisation [47]. The DE principle of operation is based on the iterative transformation of a population of candidate solutions by subjecting them to rounds of mutation, recombination, assessment and selection. New solutions are produced by perturbing the current population members with their weighted differences calculated between two or more randomly chosen individuals [48,49]. Combined with the selection operation, this has the effect of a selfscaling step when sampling the solution space providing population members with a local contour-matching property that allows fast identification of attractive basins that might contain the global optimum [50].
The surveyed literature shows that there are many DE variants and flavours of those variants [51,52]. Differences between the proposed improvements are commonly used in choosing and adapting the values of the control parameters and what search strategy or combination of strategies are applied and possibly augmented with other search algorithms [52,53]. No variant can be singled out as the best and certainly none are universal. When algorithmic performance is evaluated according to selected benchmark functions commonly used in competitions [52], the majority shows some advantage over the basic algorithm proposed in [49]. However, the original DE with fixed, appropriately set control parameter values is still competitive.
The successful application of DE strongly depends on the control parameters: population size N P , difference vector scaling factor F and the crossover rate C r [52]. However, if those could be removed while retaining the same search ability, a genuine black-box optimiser could be achieved. Here, we present a minimalistic DE variant as a straightforward black-box search algorithm to provide reliable results when optimising (t, f ) distributions. The pseudo-code listing in Algorithm 1 presents the general structure of the proposed optimiser.
The algorithm uses a single population X of N P members x i that are sampled in the g = 0 th generation within the bounds [x min , x max ] D using Latin hyper-cube sampling (LHS) [54] so that a pseudo-random initialisation is achieved. The LHS guarantees a diverse but search-space-representative population considering the population size and problem dimensionality D. A random initialisation does not produce optimal coverage and may result in crowds of similar solutions that will attract new solutions and quickly bring optimisation to a halt. The search strategy is built on the DE/current/1/ variant [48,55] as a mutation-only strategy described in Equation (17) as where y i is the generated trial solution, • is the Hadamard product, U (·) is a uniformly distributed random value, i ∈ {1, . . . , N P }, and ∆ max and ∆ min are maximum and minimum vector differences in the current generation X (g) , defined as and ∆ (g) The current target for replacement x i is chosen as the base vector for a new trial once in every generation, ensuring no selection bias is present [56]. Since no crossover is used, the algorithm is rotationally invariant; consequently, no C r parameter is needed. The difference vector range in each dimension D is uniformly sampled from the current extent of difference vectors. In addition, the direction is determined from the signum of uniformly sampled values in the (−1, 1) interval. This behaviour mimics the original algorithm but ensures diverse vectors are generated, and no mutation bias is present [56]. The process is equivalent to sampling from an infinite-size population in a bounded search space that can shrink or expand according to the population extent as the search progresses.
In line with observations about the generalised scale factor in [55], sampling from a dense population can produce difference vectors equivalent to other mutation operators of varying scale factors, and consequently, the parameter F is also not needed. The lower bound on the current range of difference vectors discourages exploitation proportionally to population spread in the search space. It provides simple diversity control to lower the chance of premature convergence by keeping the algorithm in an explorative state for longer. Since the proposed strategy can generate trial vectors with out-of-bound solutions, the mutation operator is repeated on a per-element basis until all constraints are satisfied and the trial can enter the objective function evaluation. The selection process is a local one-on-one competition with the current target x i . If better, the trial y i replaces its target in the same generation according to the asynchronous update model and the extent of the search space is updated, so it immediately affects the creation of the next trial vector. The only remaining control parameter that must be provided is the population size N P , which can be small. The population members' role is now to represent the attraction basin and serve as a base vector, while population diversity is provided by sampling difference vectors from the continuous interval.
In the present application of TFD optimisation, the population of N P = 5 members was experimentally determined as sufficient by observing the convergence rate, population diversity and the objective function variance. In addition, stopping conditions were set conservatively; the objective function variance value-to-reach was set to 1 × 10 −9 and the maximum number of generations to g max = 1 × 10 5 . The objective was to test whether reliable convergence of the proposed optimiser could be achieved on different signals using selected performance measures as the objective functions and whether those measures agree with the visual TFD assessment.
where ν 0 > 0 and τ 0 > 0 are the Doppler and lag scaling constants, and r is the tilt parameter usually bounded to the [−1, 1] interval. β and γ are coupled powers such that γ = 1/β, β ∈ [1, 2], and λ > 0 governs the slope of the pass band-to-stop band transition region. α > 0 ensures that TFD's time and frequency marginals are satisfied but since those are of lesser importance in signal processing [1], we set α = 0. The boundary constraints on the other MTEK parameters, as used in the tests, are given in Table 1.
Besides the proposed performance measure, optimisations were also performed using the volume normalised Rényi entropy with order α = 3 calculated as [1] and Stanković's ECM calculated according to Equation (23) as Equations (22) and (23) represent two of the more commonly used application-agnostic measures when objectively evaluating the performance of a (t, f ) distribution and optimising the representation. Commonly used measures are rooted in probability theory; hence (t, f ) distributions are viewed as two-dimensional probability density functions satisfying the unit sum constraint [1].
The Rényi entropy is a generalisation of Shannon's; the idea is that a smaller entropy value is associated with a better-performing distribution or a smaller number of signal components [1]. Typically, the lowest entropy order used is three since order two produces a zero entropy value for a normalised distribution. In contrast, order one is the Shannon entropy, which cannot be applied to distributions containing negative values. For oddordered entropy, oscillatory cross-terms in a QTFD cancel out, so different normalisations were introduced to capture their influence: prevalently, the signal energy-and volumenormalised versions. Given that the optimisation guided by the signal energy normalised the Rényi entropy results in the Wigner-Ville distribution as optimal, we used volume normalisation as in Equation (22) to capture the influence of cross-terms and distinguish better-performing distributions past the Wigner-Ville.
Stanković's ECM, as in Equation (23), is another prevalent measure based on lowerorder norms used in signal sparsity measurements [1]. Stanković proposed this measure as an alternative to kurtosis-type measures sensitive to varying powers of signal components. The idea is that better-concentrated distributions have shorter time and frequency support regions, and by measuring their extents, the energy concentration can be estimated, and better performance is associated with lower values of the measure. The lower-order form, as in Equation (23) with power term two, is commonly preferred since it is less sensitive to small TFD values.

Results
As test signals for the proposed algorithm, we used two synthetic and one real-life multi-component signal. For every considered signal and selected performance measure, 30 optimisation runs were conducted to assess the convergence of the optimiser and the quality of the produced TFDs. The algorithm was implemented in MATLAB 2020b, and signals were processed using Intel Haswell Core i7-4790K desktop processor with 16 GB of RAM. Most of the computational burden corresponds to calculating a timefrequency distribution via fast Fourier transform (FFT). Hence, the overall algorithm has O((NM) log (NM)) time complexity, whereas the selected measures are O(NM). However, differences in the reported average count of generations needed to reach the defined stopping condition are consequences of difficult local optima and plateaux in the objective function landscape, a function of a selected measure and kernel parameters. The results are presented in the sequel.

Linear FM Laws
In order to evaluate the algorithm's behaviour when there are only linear frequency modulations present, we considered a synthetic signal with two crossed linear FM components of unit amplitude and equal time support. The first component starts at 0.1 Hz and linearly increases with time to 0.4 Hz, while the second component starts at 0.4 Hz and decreases toward 0.1 Hz. The signal was obtained as a discrete sequence of N = 128 samples at a 1 Hz sampling rate. Its analytic form is given as where the components are The average results from 30 repeated runs are reported in Table 2, while Table 3 presents the worst and best results according to and evaluated by each considered measure. Figure 2 shows the material ratio curves of the TFDs generated according to each measure's best solution. The resulting TFDs are presented in Figure 3, along with the 0.1 per cent MTEK iso-contour overlaying the signal ambiguity function. The number of frequency bins used in generating TFDs is M = 128. The optimiser's behaviour shows consistency for each measure, as evident by comparing the average fitness in Table 2 to the respective measure's best value in Table 3, verifying that it can converge reliably.
Version June 17, 2023 submitted to Sensors 12 of 24 of the ECM result in Figure 3c. According to the RE3DV, the best PMs are second to 373 RE3DV performance. However, visual inspection of Figure 3a shows that one component 374 is missing in the RE3DV result, which is why PM labels ECMs and the WVD as better than 375 RE3DVs. In the considered example, the PM shows greater consistency than the other two 376 measures when the agreement among the perceived performances of results in Figure 3 377 and numerical indicators in Table 3 are considered.  (24). The curves correspond to the TFDs optimised according to the: RE3DV (grey), ECM (dashed grey) and PM (black). Thick grey is the WVD given for reference. The respective zero-level material ratio mr(0) related to the p 3 metric is marked on the horizontal axis. A short, steep transition, followed by an extended, nearly flat region, is characteristic of a better-performing distribution.  Table 3. In (a), component z 2 [n] (Equation (26)) is completely smeared, whereas z 1 [n] (Equation (25)) is highly concentrated. In (c), z 1 [n] and z 2 [n] are visible but lack concentration. In (e), both components are well-resolved and highly concentrated.  Table 3. In (a), component z 2 [n] (Equation (26)) is completely smeared, whereas z 1 [n] (Equation (25)) is highly concentrated. In (c), z 1 [n] and z 2 [n] are visible but lack concentration. In (e), both components are well-resolved and highly concentrated. In the parameter space, the ECM results show the smallest variance and the measure took the least amount of generations to reach the stopping conditions. On the other hand, the parameters according to the RE3DV exhibit more significant variance, which is not reflected in the objective function space. Different TFDs result from different parameter values but have similar fitness values, and there are also inconsistencies in the perceived TFD quality and reported numerical values. The RE3DV evaluates PM results as better than ECMs, but among ECMs, it labels the worst as better. According to the ECM, PMs are worse than ECMs, although a visual inspection shows a significantly lower concentration of the ECM result in Figure 3c. According to the RE3DV, the best PMs are second to the RE3DV performance. However, a visual inspection of Figure 3a shows that one component is missing in the RE3DV result, which is why the PM labels ECMs and the WVD as better than RE3DVs. In the considered example, the PM shows greater consistency than the other two measures when the agreement among the perceived performances of results in Figure 3 and numerical indicators in Table 3 are considered.

Mixed FM Laws
As, in practice, signals can have various frequency modulations, we considered a mix of two quadratic components and one linear FM component in the second test. All components are of unit amplitude and the same time support. The signal was obtained as a discrete sequence of N = 128 samples at a 1 Hz sampling rate. The analytic form is where the components are The average results from 30 repeated runs are reported in Table 4, and Table 5 shows the worst and best results according to and evaluated by each considered measure. Figure 4 shows the material ratio curves of the TFDs generated according to each measure's best solution. The resulting TFDs are presented in Figure 5, along with the 0.1 per cent MTEK iso-contour overlaying the signal ambiguity function. The number of frequency bins used in generating TFDs is M = 128. As in the case of the crossed LFMs, the optimiser shows reliable convergence for each measure, confirmed by how close the average fitness values in Table 4 are to the respective measure's best value in Table 5, which indicates it can deal with signals of mixed FM laws.    When the parameter space is considered in Table 5, the values obtained under each measure show low variance, reflected in the objective function space. Hence, no ambiguity is present-similar parameters produce similar results. However, as in the previous example, mutual evaluations show inconsistencies and are not aligned with the visual evaluation. As evident in Figure 5a, under RE3DV-guided optimisation, the MTEK emphasises linear FM segments resulting in a highly concentrated LFM component. However, two quadratic FMs are distorted, and cross-terms are amplified. Similarly, but to a lesser extent, this can also be observed in Figure 5c, which shows the resulting TFD of the ECM-guided optimisation. Here, the quadratic FMs have their shape preserved, but there is a significant loss in concentration toward the components' edge as if an amplitude modulation is present. On the other hand, the PM result in Figure 5e shows no such effects. All components are equally well concentrated, and their shapes are preserved from start to finish, corresponding to the ground truth. Nevertheless, RE3DV and ECM label PMs as second after their best results, whereas the PM values consistently reflect the perceived performance of TFDs in Figure 5.

Real-Life Signal
For a real-life test case, we used an echolocation signal emitted by the big brown bat (Eptesicus Fuscus) since it is commonly used as an example in the performance assessment of time-frequency processing methods. The signal has four nearly hyperbolic FM components with significant variations in the components' amplitude and duration. The signal length is 2.5 ms in N = 400 samples obtained at a 142.857 kHz sampling rate. The number of frequency bins used in generating TFDs is M = 512.
The average results from 30 repeated runs are shown in Table 6, while Table 7 shows the worst and best solutions according to and evaluated by each considered measure. Figure 6 shows the material ratio curves of the TFDs generated according to each measure's best solution. The resulting TFDs are presented in Figure 7, along with the 0.1 per cent MTEK iso-contour overlaying the signal ambiguity function. As in the previous examples, the average fitness values in Table 6 correspond well with the respective measure's best value in Table 7, confirming, on average, the consistent behaviour of the optimiser.  Figure 6. The material ratio curves for the echolocation signal as emitted by the Big brown bat. The curves correspond to the TFDs optimised according to the: RE3DV (grey), ECM (dashed grey) and PM (black). Thick grey is the WVD given for reference. Zero-level material ratio mr(0) related to the p 3 metric is only marked on the horizontal axis for the WVD. The near coalescence of the individual curves reflects the overall similarities of the resulting TFDs, but the far-right edge reveals significant cross-terms in the RE3DV case. Figure 6. The material ratio curves for the echolocation signal as emitted by the big brown bat. The curves correspond to the TFDs optimised according to the following: RE3DV (grey), ECM (dashed grey) and PM (black). Thick grey is the WVD given for reference. Zero-level material ratio mr(0) related to the p 3 metric is marked on the horizontal axis only for the WVD. The near coalescence of the individual curves reflects the overall similarities of the resulting TFDs, but the far-right edge reveals significant cross-terms in the RE3DV case.
The parameter space in the case of the RE3DV and ECM shows low variance, which is also reflected in the objective function space with no ambiguity in results. The worst of the PMs shows some parameter deviation, but this is appropriately differentiated in the objective function space, so there is also no ambiguity. The RE3DV and ECM label PM results as the worst or second after their best, respectively, but a closer visual inspection of TFDs in Figure 7 reveals a bias toward linear segments in the case of the RE3DV result in Figure 7a and to a lesser extent in the ECM result in Figure 7c. Consequently, subtleties in components' shapes are diminished, cross-terms appear near strong components, and weaker ones are smeared. In contrast, the PM result in Figure 7e shows equally wellconcentrated components regardless of amplitude and duration. Per the reported PM value, it is less distorted than the other two results. The worst PM result (not shown) has a four-quadrant symmetric MTEK topology similar to the one in Figure 7f that passes some cross-terms between the two weakest components, which the PM appropriately captures as lower performance.  Table 7. In (a), weaker and non-linear segments of the components are smeared, and there are cross-terms near the high-power segments. The components' shapes are better discernible in (c). However, there is some loss of concentration at the weaker segments along with the noticeable cross-term near the high-power segment of the lowest component. In (e), the components are equally concentrated in all segments regardless of varying components' amplitudes, and there are no observable cross-terms.  Table 7. In (a), weaker and non-linear segments of the components are smeared, and there are cross-terms near the high-power segments. The components' shapes are better discernible in (c). However, there is some loss of concentration at the weaker segments along with the noticeable cross-term near the high-power segment of the lowest component. In (e), the components are equally concentrated in all segments regardless of varying components' amplitudes, and there are no observable cross-terms.

Discussion
The proposed optimiser as a minimalistic differential evolution variant performed reliably in the considered examples regardless of the signal and performance measure as the objective function, signifying how algorithmic simplicity and a small population are not performance-limiting. Since it requires no additional control parameters, it represents a genuine black-box optimiser that, coupled with a reliable TFD performance measure and an adaptable kernel, could be integrated as a service in unsupervised classifiers based on (t, f ) signal processing. When the presented results are considered, it is evident how commonly used measures capture only some TFD surface complexities per their definitions, leading to inconsistent behaviour and possible deformations in the resulting TFDs, even in the noiseless case. However, the fact is that no single parameter alone can fully describe all of the surface complexities as recognised in the field of surface metrology. In the proposed measure, four indicators derived from the Abbot-Firestone curve, similar to volume parameters in surface metrology, are aggregated into a single metric. When applied in the considered examples, consistent behaviour was demonstrated by correctly capturing and numerically reflecting situations such as a missing component, increased cross-terms and low concentration. As no additional parameters are required, the proposed measure is a good candidate for a black-box TFD optimiser. However, an additional analysis should be performed for noisy signals and other combinations of components with different frequency and amplitude modulations and varying time support. Furthermore, the MTEK could be replaced by the RGK, which would be an even more significant challenge given the number of kernel parameters. However, optimised correctly, a fully adaptable kernel topology could provide even better-performing time-frequency distributions.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/s23135804/s1. presented in this study are available in the supplementary material or can be recreated based on formulae presented in this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: (t, f )