Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

Topological data analysis for true step detection in periodic piecewise constant signals

    Abstract

    This paper introduces a simple yet powerful approach based on topological data analysis for detecting true steps in a periodic, piecewise constant (PWC) signal. The signal is a two-state square wave with randomly varying in-between-pulse spacing, subject to spurious steps at the rising or falling edges which we call digital ringing. We use persistent homology to derive mathematical guarantees for the resulting change detection which enables accurate identification and counting of the true pulses. The approach is tested using both synthetic and experimental data obtained using an engine lathe instrumented with a laser tachometer. The described algorithm enables accurate and automatic calculations of the spindle speed without any choice of parameters. The results are compared with the frequency and sequency methods of the Fourier and Walsh–Hadamard transforms, respectively. Both our approach and the Fourier analysis yield comparable results for pulses with regular spacing and digital ringing while the latter causes large errors using the Walsh–Hadamard method. Further, the described approach significantly outperforms the frequency/sequency analyses when the spacing between the peaks is varied. We discuss generalizing the approach to higher dimensional PWC signals, although using this extension remains an interesting question for future research.

    1. Introduction

    Piecewise constant (PWC) signals, an important subclass of piecewise continuous data, occur in a variety of applications such as bioinformatics, astrophysics, geophysics, molecular biosciences and digital imagery [15]. Figure 1ac shows several forms of PWC signals, both deterministic and stochastic. The PWC waveform also occurs in the output of several sensors, which include tachometers. Sensors with PWC waveform are widely used and implemented in applications that involve rotating machinery such as investigating pressure fluctuations in bearings [6], flat plate motions [7], analysing the powertrain vibrations of vehicles [8] and for the experimental investigations of the dynamic performance of wind turbines [9]. One important example of PWC waveform sensors are laser tachometers, which are often used in machining process applications, including monitoring rotation velocity in the presence of spindle speed modulation [10] or vibration-assisted cutting in turning [11], defect detection for friction stir welding [12], measuring the feed error in tapping processes [13] and for stroboscopic sampling of vibrations in milling [14].

    Figure 1.

    Figure 1. An overview of the categorization (left) and analysis (right) of the type of signals considered in this study. The paper investigates periodic PWC signals with two-state trajectories and random digital ringing (b). (Online version in colour.)

    The output of laser tachometers assumes two logic levels, e.g. low/high or on/off, based on the reflectivity or contrast of the target material when subjected to a laser source. These tachometers often output transistor–transistor–logic (TTL) pulses triggered by the change in the reflected laser beam. The resulting two-state trajectory can be collected using digital or analogue channels. If analogue channels are used to collect the laser tachometer signal, then noise will be superimposed on the pulse train. Generally, using linear filtering for noise removal from PWC signals is inefficient because both the noise and the PWC signal have a broadband frequency spectrum. However, since for laser tachometers the signal to noise ratio is typically large the underlying digital signal can be recovered by hard-thresholding an adequately sampled time series. However, hard-thresholding is not suitable, in general, for stochastic PWC signals such as the case shown in figure 1c.

    Assuming that a representative, noise-free sample of the laser tachometer signal has been recovered, another challenge is the phenomenon where spurious pulses appear in the signal near the transitions of reflective/non-reflective target, as shown in figure 1b. We call these peaks digital ringing. We observed these peaks while collecting laser tachometer signals during a turning experiment; see figure 2 for a sample of real data showing this phenomenon. These spurious peaks can occur due to (i) chips or cooling fluid interrupting the laser beam in manufacturing applications, (ii) vibrations of the laser tachometer holder near regions on the target with different contrast, (iii) vibration of the target away from or towards the laser tachometer, or (iv) unintended reflections of the laser beam on non-target surfaces due to the motion of the machine. These spurious peaks can occur at the rising or falling edges but they do not occur as a repetitive pattern, and they may not occur at some peaks or at all combinations of cutting parameters. Hence, they are considered a random variable and are grouped under the umbrella of stochastic systems in figure 1. The challenge in obtaining the once-per-revolution pulse or in counting the number of pulses per unit time is exacerbated by these spurious peaks. Moreover, significant complications are introduced if the spacing between the pulses is further randomly modulated either due to external disturbances or due to some noise superimposed on intentional, regular variations of the rotational speed. In both these scenarios, the interest is in finding either (i) the rising edge, (ii) the peak, or (iii) the falling edge where a true transition in the reflectivity of the target occurs. In tachometer signals, it is common to detect falling or rising edges, so, in this study, we choose to investigate the rising edges of true transitions.

    Figure 2.

    Figure 2. Data from the laser tachometer experiment (§5), with zoom showing the digital ringing inherent in the data. (Online version in colour.)

    While finding true peaks from a signal contaminated with noise may seem like an easy problem, the existence of spurious peaks precludes using traditional, step detection methods for peak detection and counting. For example, peak-finding algorithms are not useful here because true peaks can be incorrectly counted multiple times, thus giving false counts and consequently yielding exaggerated spindle speeds. In addition, setting width thresholds in the peaks algorithm requires the user to have a priori knowledge of the existing noise, which may not be feasible or practical. The peak-finding algorithm also catastrophically fails if the spacing between peaks randomly varies due, for example, to randomness in the spindle drive or variations in the cutting load of machine tools. Further, since the signal is corrupted by digital ringing, typical filtering techniques for PWC signals [3,4], such as total variation denoising, hidden Markov chains and wavelets, are ineffective. Traditional clustering methods which require knowing the number of clusters desired in advance (i.e. k-means clustering) are not useful for our need to automatically count the number of pulses. A more unsupervised clustering approach such as DBScan [15] could be used, but requires expert knowledge input to determine input parameters.

    It is worth mentioning that identifying true peaks in the presence of digital ringing as described in this paper is a task where humans outperform traditional computer algorithms. However, while humans can often identify true peaks, counting them is prone to error, particularly when the number of pulses is large per unit time. Moreover, it is impractical to rely on human interpretation of laser tachometer signals, especially in high-speed or real-time operations.

    Therefore, the objective of this paper is to introduce a simple yet powerful approach for automatically detecting the true steps in a periodic PWC two-state square wave subject to digital ringing. The resulting algorithm is completely automatic as it requires no input parameters and comes, under reasonable noise assumptions, with a correctness guarantee. The specific application is obtaining the spindle speed of a machine tool with a laser tachometer whose signal is corrupted by spurious on/off pulses (figure 2). We also extend the approach to the case where the spacing between the pulses varies randomly, discuss an extension of the theory to higher dimensional analogues, and derive mathematical guarantees for the resulting change detection which enables accurate identification and counting of the true pulses. The approach we use based on topological data analysis (TDA) [1622] is outlined in §3.

    The prevalence of TDA has exploded in recent years due to its use in many disparate domains. Arguably the most prevalent application is that of TDA to time-series analysis and signal processing, having spawned some of the earliest results in the field [2325]. This subfield of TDA is often referred to as topological signal processing (TSP) [26]. From the insights gained from the beginnings of TSP came one of the most prominent tools in TDA, namely persistent homology [27,28], which quantifies shape and structure in data.

    Now, TSP has become a mature field in its own right. Much of the work stems from using persistent homology in conjunction with delay coordinate embeddings to give topological quantification of attractors of dynamical systems [2933]. Recent work has come in the form of classification and quantification of periodicity and quasi-periodicity [3438]. Other applications using TSP include wheeze detection [39], computer performance [40], market prices [41,42], sonar [43], image processing [44], mode hunting [45], ice core analysis [46], machining dynamics [4750] and gene expression [51,52].

    In this paper, we will use one of the simplest examples of persistent homology, namely zero-dimensional persistence defined on points in R. Viewing our data in this way gives access to the powerful theory built up for understanding persistent homology with respect to noise. In particular, the existence of a metric on the space of persistence diagrams leads to the powerful notion of stability [53,54]. This knowledge in conjunction with our assumptions on noise is the theoretical basis for the algorithm we develop in §3. While there is a growing collection of ever faster code for computation of persistent homology [55], our restricted setting also gives rise to simplified algorithms, thus making the analysis quite fast. An additional perk of this viewpoint is that it also allows for generalization to higher dimensional problems as outlined in §6a.

    The approach is described and tested using both synthetic and experimental data. The experimental apparatus includes an engine lathe instrumented with a laser tachometer which detects the change in reflectivity of a tape adhered to the circumference of the spindle. The resulting mean spindle speed is then reported with appropriate error bounds. We also compare the results with the output of Fourier analysis (frequency domain) and the Walsh–Hadamard transform (sequency analysis). It is found that the described approach provides comparable results to Fourier analysis for regular pulse waves with digital ringing; however, the numerical calculations show that the described approach outperforms both Fourier and Walsh–Hadamard analysis when the spacing between the pulses is varied. Furthermore, we found that digital ringing has a strong negative influence on the Walsh–Hadamard analysis and that both Fourier and persistence approaches show better results, with the latter showing the best performance as the digital ringing becomes more and more severe.

    2. Background

    This section provides the necessary background for motivating and presenting the new approach for pulse counting using persistence diagrams. Section 2a discusses existing methods for pulse counting and spindle speed calculations. Sections 2b–d present the background theory on persistent homology and describe how it relates to the current work.

    (a) Current methods for counting and rotational speed calculation

    Sensors for detecting rotary motion include proximity sensors, photoelectric sensors and encoders. These sensors output pulses that can be counted and used to find the speed of the shaft, typically in units of revolutions per minute (RPM). The pulse signals are PWC functions with two logic levels: high or on, and low or off. The quality of the output data depends on the number of pulses per revolution, which affects the data resolution, as well as the symmetry of the pulses, which influences the accuracy and consistency of the data.

    Once the pulses are obtained, there are generally two techniques for determining the corresponding RPM: (i) frequency measurement approach (calculate RPM from pulse count and (a) pulse frequency or (b) pulse sequency), and (ii) period measurement approach (calculate RPM from pulse count and pulse period). The frequency approach involves transforming the signal into the frequency domain using a Fourier transform or into the sequency domain using the Walsh–Hadamard transform; however, the Fourier transform of a PWC signal can have slow convergence [56], and both Fourier and Walsh–Hadamard transforms often do not perform well for sparse data [57], such as the data output from a laser tachometer. Other counting algorithms include local maxima or peak detection with wavelet transforms or other methods [3]. However, the interest in this study is in detecting true peaks, not all peaks. Therefore, using conventional methods for peak detection must be combined with a threshold for rejecting false peaks or retaining true ones. It may be tempting to use a statistical measure such as the variance or standard deviation of the pulse duration. However, since the data can vary from bimodal to uniform depending on the amount of noise introduced, statistical dispersion measures are generally not effective.

    When there is noise superimposed on the signal (for example, due to collecting the signal on an analogue channel) this noise can be removed by hard thresholding. If the noise component is so large that the pulse structure can no longer be distinguished as shown in figure 1c, the denoising becomes more difficult because the traditional approach of low-pass filtering typically introduces large spurious oscillations in PWC signals [56]. In this case, either the viewpoint of piecewise constant smoothing or that of level-set recovery can be used for determining the location of the jumps [3]. In this paper, we focus on signals where the small amount of noise superimposed on the pulse amplitude can be removed by thresholding, but where spurious random pulses occur at the rising or falling edges and where the spacing between peaks can randomly vary. This type of noise cannot be removed by thresholding, or by traditional filtering techniques. Varying the period length between two consecutive pulses significantly complicates the analysis, and without a proper method for change detection can lead to poor RPM calculations. Therefore, there is still a need for new, robust tools for pulse detection in PWC signals and this paper presents a method for reliable pulse counting using zero-dimensional persistent homology.

    (b) Zero-dimensional persistent homology

    Persistent homology1 [23,24,27,28], a tool arising from TDA [21,22], seeks to quantify shape and structure in datasets. Algebraic topology [58,59] is a field of mathematics which quantifies qualitative similarities in the structure of spaces. One such method for quantification is homology, which, given a topological space, provides a vector space2 for each dimensional structure being studied. This paper will only focus on dimension 0, which quantifies connected components.

    The intuition behind persistent homology for a point cloud dataset is to increase a connectivity parameter, and quantify how the topological structure changes. This powerful method can find interesting, higher dimensional structure, using each dimension of homology; however, we will look at the simplest version for the purposes of our problem. Zero-dimensional persistent homology quantifies how the clusters change when viewed at different scales. In fact, zero-dimensional persistent homology is closely related to classical clustering methods such as single-linkage hierarchical clustering, dendrograms and minimal spanning trees [60]. Here, we adopt the view of the procedure as a restricted case of persistent homology in order to gain understanding and predictive power in our analysis, thus resulting in theorem 3.1.

    In the general setting, assume we are given a point cloud χRD with |χ| = n. We can define a function fχ:RDR by fχ(x) = 2∥x − χ∥, where ∥x − A∥ = infyAx − a∥ for any set ARD. The set of points for which fχ(x) ≤ r is the union of D-dimensional balls of radius 12r centred at the points of χ; we write this as fχ−1( − ∞, r]. If we allow r to increase from 0, initially, f−1χ( − ∞, r] has n distinct connected components. However, as r is increased, these components will merge together until we are finally left with only one connected component. In particular, these mergings happen at the instant two discs touch, and thus when r is equal to the distance between the associated points. We say that a connected component dies when it merges with another connected component; that is, a death occurs any time two clusters merge. We can keep track of the function values at which these deaths occur in the following manner.

    A function value r is a homological critical value if the number of connected components decreases at function value r. The multiplicity of a critical value is the net decrease in the number of connected components. We define dgm(χ)R to be the collection of homological critical values, with the number of copies equal to the multiplicity. The resulting set of values is called a zero-dimensional persistence diagram, notated dgm = {d1 ≤ d2 ≤  · s ≤ dk}.3 One observation that will be useful for interpretation is that the number of connected components of f−1( − ∞, r] is one more than the number of di∈dgm with di > r.

    We will also use the slightly more general formulation of persistence for any subset ARD. If fA:RDR is given by fA(x) = 2∥x − A∥, dgm(A) is the collection of function values for which the number of connected components of f−1( − ∞, r] changes, again with multiplicity. In order to ensure that the persistence diagrams are finite, we assume that A has finitely many connected components.

    (c) Stability of persistence diagrams

    The main reason for thinking of this information as a persistence diagram is that these structures come with a metric known as the bottleneck distance. Here, we describe the metric in the restricted zero-dimensional persistence diagram setting; see, for example, [17] for the full definition. Given two zero-dimensional persistence diagrams, dgm1 = {d1 ≤ d2 ≤  · s ≤ dk} and dgm2 = {e1 ≤ e2 ≤  · s ≤ e}, a partial matching η is a bijection between subsets of the two diagrams η:A → B, A⊂dgm1, B⊂dgm2. The cost of a matching is defined to be

    c(η)=max({|aη(a)|}aA{a2}adgm1A{b2}bdgm2B),
    and the bottleneck distance dB(dgm1,dgm2)=minηc(η) is the minimum cost of the possible matchings. Note that, because we are restricting our diagrams to finitely many points, the set of matchings is finite, thus this minimum is always achieved. In addition, it is possible that multiple matchings achieve the minimum, so we will reference such a matching as a min-cost matching.

    Consider, for example, a diagram consisting of a single point dgm1 = {d1} versus a diagram with two points dgm2 = {e1, e2}. Without loss of generality, assume e1e2. The only possible matchings are where η1(d1) = e1, η2(d1) = e2, or η3, which matches nothing. The scores of these are c(η0)=max{d1/2,e1/2,e2/2}, c(η1)=max{|d1e1|,e2/2} and c(η2)=max{|d1e2|,e1/2}. In the example of the two persistence diagrams in figure 3b, the score is lowest using η1, so the distance between the diagrams is max{|d1e1|,e2/2}.

    Figure 3.

    Figure 3. (a) An example point cloud with its zero-dimensional persistence diagram drawn as a histogram at the centre. Two points in the diagram above the split translates to having three clusters in the point cloud. (b) An example of the pairing for two persistence diagrams used to compute the bottleneck distance. (Online version in colour.)

    The bottleneck distance is particularly useful due to the stability theorem [53]. Recall that the Hausdorff distance between sets A,BRD is

    dH(A,B)=max{supaAaB,supbBbB}.
    Then the stability theorem is as follows.

    Theorem 2.1 ([53])

    Under mild assumptions4on the setsA,BRD,

    dB(dgm(A),dgm(B))dH(A,B).

    This theorem is particularly useful when we take the view, as we will need to in §3, that B is a noisy point cloud approximation of A. In this case, dH(A, B) is small, so the resulting persistence diagrams will be close as well.

    Let us now explore what it means for a diagram to be close to another diagram. Assume we again have a diagram with a single point dgm1 = {d1}. For ϵ small relative to d1, we can think of the set of diagrams within distance ϵ of dgm1 = {d1}. Say the second diagram is dgm2 = {e1, …, en}, so the options for matchings are ηi matching d1 to ei for i = 1, …, n, and η0 matching nothing. Assume dB(dgm1, dgm2) ≤ ϵ and let ϵ < d1/3. If a minimum cost matching is ηi, this implies |d1 − ei| ≤ ϵ, and ej ≤ 2ϵ for all ji. This means that there must be exactly one point, ej, within distance ϵ of d1, and all remaining points are in [0, 2ϵ]. As shown in the rightmost portion of figure 3, this means there is exactly one point in the top green region, and any remaining points are in the bottom green region.

    What will be useful in §3b is the set of diagrams close to the diagram which has n copies of the point d1, e.g. dgmn = {d1, d1, …, d1}, . In this case, any diagram within bottleneck distance ϵ (again sufficiently small ϵ < d1/3) will have exactly n points in [d1 − ϵ, d1 + ϵ], with all remaining points less than 2ϵ.

    (d) Computation of the persistence diagram using minimal spanning trees

    In order to do computations, we convert this information into combinatorial structures. Following [23,24], we compute persistence using a minimal spanning tree. A graphG = (V (G), E(G)) consists of a finite list of vertices, V (G), and a set of edges, E(G), between them. A graph is complete if every pair of vertices has an edge between them. A subgraphAG is a subset of G, which is itself a graph. That is, A = (V (A), E(A)), where V (A)⊆V (G), E(A)⊆E(G).

    A path in G is a sequence of vertices v0, …, vn such that there is an edge between every adjacent pair: (vi, vi+1)∈E(G) ∀i. A path is a cycle if v0 = vn. A graph is connected if there is a path between every pair of vertices. A graph is a tree if it is connected and there are no cycles. A subgraph TG is called a spanning tree if it is a tree and V (T) = V (G).

    A weighted graph G = (V (G), E(G), ω) is a graph with a real value, called a weight, associated with each edge: ω:E(G)R. The total weight of a spanning tree of a weighted graph is the sum of the weights on the edges. A spanning tree is called a minimal spanning tree (MST) if it has minimal total weight among the set of spanning trees. If G is not connected, no MST exists as there is no connected subgraph using all the vertices. A minimal spanning forest (MSF) is an acyclic subgraph such that the restriction to each connected component of G is an MST.

    Given a finite set of points χRD, the Čech graph for parameter r > 0 is the complete graph Č(χ, r) with the vertex set in 1-1 correspondence to the points in χ, with all edges (x, y) for which ∥x − y∥ ≤ r, and with edge weight equal to the distance between the associated points. First, note that, if r ≤ s, Č(χ, r)⊆Č(χ, s). Second, if r is larger than the diameter of χ (which is finite since χ is finite), Č(χ, r) is a complete graph and Č(χ, r) = Č(χ, s) for all sr. Thus, we denote this graph as Č(χ, ∞).

    Our first intuition for understanding the connected components of the union of discs as described above is to watch the connected components of Č(χ, r) change as r increases. It is a consequence of the celebrated nerve lemma [58] that the connected components of Č(χ, r) match up with the connected components of the union of discs of radius r/2; equivalently with the connected components of fχ( − ∞, r] in the notation of the previous section. So, we just need to determine the list of function values for which the connected components of Č(χ, r) change.

    Let T be an MST for Č(χ, ∞)5 . It is an immediate consequence of [24] that the list of weights on the edges of T is exactly dgm(χ). Restated in our notation, this is the following.

    Lemma 2.2 ([24, Lem. 3])

    Let T be an MST for Č(χ, ∞), and let Tr be the restriction to the set of edges of weight at most r so that Tr ⊂ Č(χ, r]. Then Tr is an MSF for Č(χ, r].

    Since an MSF of a graph has the same connected components as the graph itself, this means that the changes in the connected components happen exactly at the values of the weights in T. Putting this all together, we have the following theorem.

    Theorem 2.3

    For a finite set of pointsχ, dgm(χ) is the set of weights on the edges of an MST of Č(χ, ∞).

    The final thing to note is that we will largely be interested in points in R (so D = 1). Say χ = {a0, …, ak}. In this case, it is obvious that the MST is exactly the graph with edges {(ai, ai+1)}k−1i=0, and so the diagram dgm(χ) = {(ai+1 − ai)}k−1i=0.

    3. Method

    (a) Basic model assumptions

    The initial data are a time series X defined on {a = t0 < t1 < · s < tN = b} given by {X(t0), X(t1), …, X(tN)}. We assume that these data approximate the pulse wave PτT given by

    PTτ(t)={1if 0(tmodT)τ,0else,
    and shown in figure 5. We will assume that the duty, δ = τ/T, is less than 0.5 since otherwise we could swap our analysis to check for the pulses at 0 instead of the pulses at 1.

    We first assume that our noisy data are given by

    X(t)=PTτ(t+δx)+δy,3.1
    where δx∼unif( − α · τ, α · τ) and δy∼unif( − β, β) with α[0,12] and β∈[0, 1]. For simplicity, we assume the noise in y, δy, separates the high and low values, so β12. The noise in x, δx, serves to give incorrect pulse values when we near the beginning or end of a pulse by looking forward and backward in time to reach an output. Since α is given as a percentage of τ, and we will not be able to see anything if α is large enough to entirely erode the off-cycle with digital ringing, we assume α ≤ (T − τ)/3τ. For practical applications with a small τ, this is satisfied for any α[0,12].

    It will be important later to note that, with these assumptions on β and α, the Hausdorff distance between M and χ={tX(t)12} for any realization of {X(t)} dense enough in t will be bounded by ατ.

    (b) Counting pulses and revolutions per minute calculation using persistence

    Counting the number of pulses in a range [A, B] can be thought of as counting the number of connected components of M = (PτT)−1({1})∩[A, B]. With this in mind, we use a special case of zero-dimensional persistence to determine this from the point cloud approximation of M, namely χ={tiX(ti)>12}. See figure 4 for an example.

    Figure
4.

    Figure 4. (a) A pulse train, with T = 2, τ = 0.8, α = 0.2, β = 0.1. The corresponding sets χ and M are drawn above. The choice of domain for the RPM computation is shown by the vertical (red) lines near 4 and 6 in the pulse train. (b) The histogram for the zero-dimensional diagram. There is a large split between the non-empty bins; K = 2 points in the diagram above the split implies that there are K + 1 = 3 pulses seen in the full window, and K − 1 = 1 pulse seen within the vertical red lines in the pulse train. (Online version in colour.)

    Figure 5.

    Figure 5. (a) An ideal pulse train PτT. (b,c) Noisy instances of the pulse train, both with T = 2, τ = 0.5, α = 0.05 and β = 0.05. Panels (b,c) have ϵ = 0 and ϵ = 0.5, respectively. Note the uneven pulse spacing for ϵ≠0. (Online version in colour.)

    In order to determine the persistence-based RPM ΩP, we need to find the range on the t-axis containing the full periods, and count the number of full periods there. In some sense, this procedure is measuring the voids between the components of M, rather than measuring the components of M themselves. The following theorem essentially outlines an algorithm, which, given controlled enough noise, dense enough sampling and enough visible pulses, calculates the number of pulses in that region using the persistence diagram dgm(χ).

    Theorem 3.1

    Assume a time seriesX(0), X(t1), …, X(tN) is data drawn from the model of equation (3.1) withα<15(T/τ1),andβ < 0.5,and with times {A = t0 < t1 < · s < tN = B} evenly spaced such thatti − ti−1 < τ(1 − 2α)/2 and (B − A)/T > 3. Setχ={tiX(ti)>12}={a1am}and denote the resulting persistence diagram dgm(χ) = {b1 ≤ · s ≤ b}. Chooseμ∈(dj, dj+1) forj = argmaxk{dk+1 − dk}. Then |{d∈dgm(χ)|d > μ}| − 1 is the number of pulses in the range

    A<Alow:=(AτT+1)T<BTT=:AhighB.3.2

    Proof.

    First, note that, if M has m-connected components, there will be m − 1 mergings of components that occur exactly at distance T − τ. Thus, the persistence diagram of M, dgm(M), has m − 1 copies of d1 = T − τ. By definition of the noise parameters α and β, we know that the Hausdorff distance between M and χ is at most ατ. So, by the stability theorem, theorem 2.1, we have that dB(dgm(M), dgm(χ)) < ατ.

    Thus, we turn our attention to computation of the zero-dimensional persistence of χ. As discussed previously, the MST for points in R has edges {(ai, ai+1)}k−1i=0. Hence, the set of pairwise distances {a1 − a0, a2 − a1, …, an − an−1} is exactly the zero-dimensional persistence diagram dgm(χ).

    The set of possible zero-dimensional diagrams within bottleneck distance ατ of dgm(M) have exactly m − 1 points in [(T − τ) − ατ, (T − τ) + ατ], and all remaining points are in [0, 2ατ]. Requiring α<15(T/τ1) means that there will be a distinct split in dgm(χ); in particular that the largest difference between points in the diagram will occur between a point in [0, 2ατ] and a point in [(T − τ) − ατ, (T − τ) + ατ]. These points are then notated as dj and dj+1 where j = argmaxk{dk+1 − dk}.

    We choose a threshold μ between the collections, and count #{d∈dgm(χ)|dμ} + 1; that is to say, the number of points in the diagram in the interval [(T − τ) − ατ, (T − τ) + ατ], so this is m − 1.

    While it is an annoying case study in number theory to determine exactly how many pulses appear in the interval [A, B], what we can say is that if we cut off the first (possibly partial) pulse seen, and the last (possibly partial) pulse seen by looking in the interval [Alow, Ahigh] defined in equation (3.2) then the number of pulses is m − 2. Thus, the theorem follows. ▪

    In particular, since Alow and Ahigh can be approximated within ατ by alow=min{aiχaiai1>μ}andahigh=max{aiχaiai1>μ}, this theorem allows us to approximate the RPM in the following sections as

    ΩP=#{ddgm(χ)dμ}1ahighalow.3.3
    Putting this together with theorem 3.1 gives the algorithm.

    Inline Graphic

    (c) Counting pulses and revolutions per minute calculation using Fourier transform

    The pulse counting and RPM calculation can be obtained by transforming the signal to the frequency domain. Specifically, the one-sided Fourier transform of the signal is obtained after subtracting the mean value from it to remove the DC component from the spectrum. Then, the frequency f in Hertz corresponding to the first amplitude A that satisfies A > Amax/w is found, where w is a scalar, and the Fourier-based RPM ΩF is calculated according to

    ΩF=f1×60.3.4
    In this study, we chose w = 3 since it provided accurate results for the studied cases with constant peak spacing. However, other values of w did not improve the accuracy of the method for the cases with non-constant peak spacing, and therefore the deterioration in the performance of the Fourier-based approach for these cases was not related to our choice of w = 3.

    (d) Counting pulses and revolutions per minute calculation using the Walsh–Hadamard transform

    If the data take only two values similar to the PWC signals in this study, then the Fourier transform may not be the optimal choice since it requires an infinite number of terms to faithfully capture the signal. Therefore, in the case of these binary signals the Walsh–Hadamard transform can be more advantageous because it represents the signal using rectangular bases functions [57,61,62]. Formally, the discrete Walsh–Hadamard transform XW of a discrete signal x of N points is given by

    XW[m]=1Nn=0N1x[n]Wm[n],m=0,1,2,,N1,3.5
    where Wm is the Walsh sequence of order m. In contrast to Fourier bases, the rectangular Walsh bases Wm generally have humps of unequal length; therefore, they provide sequency information instead of frequency. Sequency is a generalization of frequency and it either denotes the (average) number of zero crossings of a Walsh basis or the number of humps in a unit of time [61]. In this study, we will use sequency to denote the number of humps in the basis function. While frequency provides a measure of the periodicity of goniometric functions, sequency provides a measure of the periodicity of Walsh functions.

    After transforming the data using the discrete Walsh transform, we can obtain a sequency PW spectrum which is analogous to the power spectrum of the Fourier transform. This spectrum can be obtained from the Walsh–Hadamard transform XW, ordered by sequency, according to [62]

    PW[0]=XW2[0],3.6
    PW[m]=XW2[2m1]+XW2[2m],m=1,2,,N213.7
    andPW[N2]=XW2[N1].3.8

    Then, by examining the sequency with the largest amplitude we can glean information about the number of pulses in the signal. Dividing the number of pulses by the corresponding time interval gives the Walsh–Hadamard RPM ΩW. Note that, in this study, we performed the Walsh transform on the signal after subtracting its mean to remove the DC component (constant function value), which corresponds to sequency 0.

    (e) Accordion model assumptions

    We now create a generalized model for testing purposes where we have data which do not maintain regular pulses. This model will generate data on a window [0, W]. Assume Q1, …, QK are drawn iid from unif((1 − ϵ)T, (1 + ϵ)T). We assume ϵ∈[0, 1] and K > W/(1 − ϵ)T is large enough to ensure that i=1KQiW. We will generate pulses so that the length of the ith period is Qi.

    For the sake of notation, let Q0 = 0. For any s∈[0, W], define σ(s)=max{j[0,,K]i=0jQis}. The reparameterization function is given by

    ϕ(s)=(Tσ(s))+TQσ(s)+1(si=0σ(s)Qi).
    Then the noisy time series is given as
    X(s)=PTτ(ϕ(s)+δx)+δy,3.9
    where δx and δy are noise as defined in §3a. See figure 5 for an example. Note that this model simplifies to equation (3.1) when ϵ = 0.

    4. Numerical simulations and robustness analysis

    This section describes the numerical simulations and shows the results of the robustness analysis. The basic idea is to take a standard pulse train with a small duty cycle and introduce spurious peaks at both the rising and the falling edge of the data using a uniformly distributed random variable as defined in equation (3.1). This noise is found to mimic the digital ringing noted in the laser tachometer experimental data discussed in §5. Furthermore, we simulate data from a signal where the spacing between the peaks is non-constant as defined in equation (3.9). The signal with modulated peak separation simulates, for example, a heavy cutting process where the spindle speed varies during the cut. Alternatively, this spindle speed variation might be intentionally introduced to mitigate chatter and improve the cutting process [6366]. The basic parameters used in the simulation are shown in table 1.

    Table 1.The parameters used in the numerical simulations.

    parameterdescriptionvalue(s) studied
    δnominal duty cycle for the pulse5%
    Ω0nominal RPMs[30, 24 000]
    Tnominal period of the pulse in seconds60/Ω0
    τtime in on cycleδ · T
    αuniform noise in x as % of τ[0, 0.5]
    βnoise in y as % of amplitude[0, 0.5]
    ϵpeak spacing uniform noise amplitude/T[0.02, 0.65]
    nnumber of periods to simulate32
    moversampling factor32
    Nnumber of simulation pointsN = m(Ω0/60)nT

    In order to test the robustness of the persistence-based algorithm, the described numerical simulations were used in a number of analyses, which include:

    • (i) Investigating the described approach as a function of digital ringing. This is controlled by the parameter α, which specifies, as a percentage of τ, the interval around the rising and falling edges that will be affected by digital ringing. The digital ringing is introduced using a uniformly distributed variable defined on the interval described by α. Therefore, as α gets larger, the digital ringing gets more severe.

    • (ii) Studying the performance of the method when varying the spacing between the peaks (the accordion effect). This is controlled by the parameter ϵ, which represents, as a percentage of T, the width of the interval that will be modulated by uniform noise. The larger the ϵ, the more likely it is for the peaks to get closer.

    • (iii) Studying the performance of the method when varying the amount of noise in y. This is controlled by the parameter β. A large value of β increases the probability that the 0-portion of the pulse will be split into two.

    • (iv) Analysing the results of our approach versus the well-established Fourier transform method. With Fourier the frequency is found and the RPM or the count of pulses per unit time are then deduced from the prominent frequency components.

    • (v) Benchmarking the computational time of the persistence-based approach versus the heavily optimized Fourier transform approach.

    All simulations and related analyses were performed using Python's scipy stack. The numerical simulations consisted of a 100 × 100 uniform grid in each of the (Ω0, α), (Ω0, ϵ) and (Ω, β) planes. The grid boundaries are given in table 1. For all simulations, 100 replicates were generated at each grid point. The starting seed for the random number generator was 48 824, and was incremented by 1 for each subsequent iteration. For each grid point, the RPM was calculated using both persistence (equation (3.3)) and Fourier (equation (3.4)). The computation time for each method was also recorded using the performance counter from Python's time package.

    The sample mean was used as the point estimator on each grid point. The relative error between the output of each algorithm (Fourier, and persistence) and the nominal RPM Ω0 was computed. A comparison of results between Fourier and persistence is summarized in figure 6a,c and b,d, respectively. Figure 6a,b shows the results for varying the accordion parameter ϵ, while figure 6c,d corresponds to varying the digital ringing parameter α with colour maps fixed for each row.

    Figure 6.

    Figure 6. Heat maps showing a comparison of the relative errors between the calculated RPM using persistence (a,c) and Fourier (b,d). (a,b) The comparison in the (Ω0, ϵ) plane, while (c,d) compares the two algorithms in the (Ω0, α) plane. The colour scale is fixed for each row to facilitate the comparisons.

    (a) α versus Ω0

    Figure 6c shows that persistence performed significantly better than Fourier. Specifically, the relative error of the persistence algorithm is below 7% for α < 30%. By contrast, the error in Fourier is as high as 24% for α < 30% and goes to at least 40% for α > 35%. A closer look at the heat map for persistence in figure 7 shows how the error in the persistence algorithm is between 15% and 25% but only for α > 30%. These are large values that correspond to a signal with very strong variations between the peaks. Even though we show the results for α∈[0, 0.5], we could deal with the case of α > 0.5 by reversing the persistence algorithm from looking at the gaps between peaks (logic 0) to looking at logic 1 (high values in the two-level digital signal). Therefore, the user has the flexibility to adjust the algorithm to recover good performance, as needed.

    Figure 7.

    Figure 7. (a) The heat map of relative error in the persistence-based RPM as a function of Ω0 and α showing more detail (equivalent to figure 6c). (b,c,d) The calculated RPM using persistence and Fourier algorithms for ϵ = 25% with nominal RPMs 3178, 12136 and 20126, respectively, which are marked using horizontal black lines. The 68% confidence bands are shown.

    Figure 7 further shows the plot of the average RPM versus α calculated using both persistence and Fourier for ϵ = 25% at three nominal RPMs: 3178, 12 136 and 20 126. The curves for the calculated persistence and Fourier RPM values also include the 68% confidence band computed using an empirical bootstrap distribution of the sample mean. It can be seen that, for α < 20%, Fourier consistently underestimates all three nominal RPM values. Moreover, as α is increased beyond 20%, the Fourier estimate starts to decrease further below the nominal RPM. By contrast, persistence gives close estimates of the nominal RPM for α < 25%, and only tends to consistently undershoot the nominal RPM for α > 25%. Thus, we can conclude that the persistence algorithm is robust to the width of the digital ringing.

    (b) ϵ versus Ω0

    Figure 6a,b shows that for ϵ < 7% the two algorithms produce comparable results with less than 10% relative error. However, as ϵ is increased, the performance of Fourier progressively deteriorates from 10% to greater than 40% as evidenced by the increase in the relative error across approximately four horizontal bands with a width of 10% each. These bands are roughly given by ϵ ≤ 15%, 18% < ϵ ≤ 30%, 30% < ϵ ≤ 40% and ϵ > 40%. By contrast, the heat map in figure 6a shows that persistence consistently maintains less than 10% error for almost all values of ϵ; however, errors start to increase for ϵ > 42%. The deterioration of the performance of the described approach is better observed from the heat map in figure 8, where higher relative errors start occurring towards the upper limit of the ϵ values.

    Figure 8.

    Figure 8. (a) The heat map of relative error in the persistence-based RPM as a function of Ω0 and ϵ showing more detail (equivalent to figure 6a). (b,c,d) The calculated RPM using persistence and Fourier algorithms for α = 10% with nominal RPMs 3178, 12136 and 20126, respectively, which are marked using horizontal black lines. The 68% confidence bands are shown.

    In addition to the relative error as a function of ϵ and Ω0, figure 8 shows the plot of the average RPM versus ϵ calculated using both persistence and Fourier for α = 10% at three nominal RPM values: 3178, 12 136 and 20 126. The curves for the calculated persistence and Fourier RPMs also include the 68% confidence band computed using an empirical bootstrap distribution of the sample mean. We can see that for both algorithms the confidence bands are narrow. Further, for ϵ < 5% both algorithms are close to the nominal value for all three RPM choices. However, for ϵ > 5%, the Fourier algorithm quickly diverges away from the nominal RPM while persistence remains close to the nominal value. As ϵ is increased to 45%, we start to see persistence undershooting the nominal value, which indicates that, at high values of ϵ, persistence undercounts the true peaks. Therefore, we can conclude that the persistence algorithm is robust to the variation in the spacing between the peaks.

    (c) β versus Ω0

    For all previous numerical simulations, we have set β = 0 in order to discount the effect of noise in y. Theorem 3.1 assumed the noise in y was drawn from a uniform distribution. The proof relies on the control of β in order to promise correct RPM computation. In particular, the only thing needed for that theorem is that δy is bounded. Still, it is reasonable in practice to expect that δy comes from a Gaussian distribution, or at least a distribution with heavier tails than uniform.

    To this end, we have tested the persistence-based method as we vary β. In this test, we assume δy is drawn from a normal distribution N(0, β). The code for persistence-based RPM returns not a number (NaN) if no split is found in the persistence diagram. As can be seen in figure 9a, there is a sharp transition where the persistence code goes from returning an answer with low relative error to returning NaN. In figure 9a, the results are averaged over all 100 replicates, so the result is NaN if any of the replicates returned NaN. This behaviour is to be expected for the algorithm since a single t, particulary one near the middle of the 0-portion of a pulse, for which the noise in y is large enough to flip the outcome to a 1 after thresholding, and thus including t in χ, will replace one point in the persistence diagram above the wide split into two points right in the middle of the split. The more of these occurences that happen, the less likely that the persistence method will be able to automatically differentiate between the subsets of points in the diagram.

    Figure 9.

    Figure 9. The effect of changing β on relative error, where β is the standard deviation of a normal distribution, δyN(0, β2). (a) White space means any one of the tests returned RPM as NaN, (b) the average relative error over those tests which return an answer and (c) the computed RPM using Fourier.

    In figure 9b, we look only at the average relative error of those replicates which returned an answer. We can see that the results are similar to Fourier until approximately β = 0.2, but increase more steeply after that. Thus, we can conclude that, so long as the variance is small, no issue is to be expected with the persistence-based algorithm in the face of noise in y. However, algorithm ?? does not fail gracefully when confronted with large variance.

    (d) Runtime comparison

    In terms of theoretical runtime, the algorithm for computing ΩP is simple, requiring only a sequence of pairwise subtractions and sortings. Since these take O(n) and O(nlogn), respectively, the theoretical worst case runtime for the persistence-based algorithm is O(nlogn). Meanwhile, the fast Fourier transform (FFT) algorithm computes the Fourier transform in O(nlogn) as well, so the worst-case analysis of the two methods is the same.

    Runtime benchmarks for Fourier and persistence algorithms were measured using the performance counter from Python's time package. The machine used for the simulations was a standard desktop running Windows 10 with 16 GB ram and a 4.2 GHz Intel i7 processor. A total of 200 runs for each algorithm was performed at each value of the nominal RPM. The nominal RPM is kept as a variable of runtime because it controls the length of the simulated pulse (table 1), and, consequently, the runtime for each algorithm. The different runs were then averaged for each algorithm, and the results are plotted in figure 10. The 68% confidence interval obtained by bootstrap sampling of the replicates is included in figure 10, but is too small to be visible. Note the small y-axis scale in figure 10, indicating that both algorithms run relatively fast. However, the figure shows that persistence runs slightly faster than Fourier, and that the runtime is approximately constant across the RPM range. By contrast, the runtime for Fourier varies more than its persistence counterpart for the different values of the nominal RPM.

    Figure 10.

    Figure 10. Average runtime for all trials plotted on a logarithmic scale as a function of the nominal RPM, i.e. the signal length. (Online version in colour.)

    5. Experimental verification

    The experimental apparatus is shown in figure 11. It consists of a Terahertz Technologies LT-880 Laser Tachometer mounted against the spindle of a Clausing-Gamet 33 cm (13 inch) engine lathe. The laser beam is emitted onto a point along the circumference of the spindle, which is half covered with white tape and half covered with black tape. When the emitted light beam is reflected back to the tachometer's receiving lens, a 5 V TTL pulse is registered. These pulses are captured by an NI USB-6366 data acquisition box using Matlab. The data are collected using an analogue channel to emulate a common practice in industry where sensory inputs in machining processes (including those of laser tachometers) are all collected using analogue channels of the same data acquisition device. The signal to noise ratio for laser tachometers is large and it is easy to recover a digital signal by hard-thresholding. Using the portions of the signals before the first pulse close to 5, the noise in y was determined to have standard deviation 2 × 10−5 as a percentage of pulse height. In this study, additive normal noise was removed by setting the threshold at 2.5 V. Thus, any signal above the threshold was set to 1, while signals below the threshold were set to zero.

    Figure 11.

    Figure 11. The experimental apparatus consisting of a laser tachometer mounted against a rotating spindle. (a) Back view, (b) the side view. Half of the spindle is fitted with white tape whereas the other half is covered with black tape. (Online version in colour.)

    The experiment was performed using several nominal spindle speeds, which are shown in table 2. Three trials were performed for each speed and data collection started after the initial spindle ramp up ended. The data were sampled at 80 kHz and the duration of the data collection was chosen such that for each trial the recorded time series captured 30–32 true peaks in the tachometer's pulse train (table 2).

    Table 2.The nominal spindle speeds (top row) and the corresponding duration of the data collection for each trial (bottom row). Three records were collected at each nominal RPM.

    nominal RPM
    30405472981301752353204255707701030
    62503530201512975433
    duration of data record (s)

    The reason for choosing a sampling rate of 80 kHz is twofold. First, it is now becoming standard to use digital signal processing (DSP) for low/high/band-pass filtering instead of in-line analogue filters. A conservative rule of thumb in DSP is to oversample by a factor of 16, digitally filter the signal, and then downsample to the frequency range of interest. As many signals in cutting processes are collected to study or prevent chatter, which can occur at frequencies of a few thousand hertz, we wanted to make sure that our results apply to cases where by the time all the signals being collected are downsampled by a factor of 16 we can recover a noise-free signal within the range of 5 kHz; so 16 × 5 = 80 kHz. The second reason is that the frequency roll-off of the laser tachometer we used is 40 kHz. Therefore, using 80 kHz ensured that we were using the full dynamic range of the sensor. This eliminates the possibility of attributing the digital ringing we observed in the signal to insufficient sampling. Although we report the results when sampling at 80 kHz, other, lower, oversampling ratios showed the same behaviour.

    Figure 12a shows the results of calculating ΩP, ΩF and ΩW using the collected signals. Note that the experimental data were collected under a no-load condition on the lathe and therefore it is free of the accordion parameter, i.e. these data are most closely related to equation (3.9) with ϵ = 0. It can be seen that Fourier and persistence results are almost indistinguishable. Further, close examination of the two inset panels in figure 12 shows that the two curves can cross at several points and trade places. Nevertheless, they remain close and parallel for all practical purposes. On the other hand, ΩW follows the nominal RPM curve closely until Ω0 = 98. Beyond this point, ΩW shows a larger deviation from Ω0, ΩF and ΩP. However, ΩW sporadically approaches the nominal values again at Ω0∈{235, 425, 570, 1030} RPM. A grid is superimposed on the figure to show the expected points on a curve that represents perfect RPM calculation versus Ω0, the nominal RPM set on the machine. It can be seen from the figure that both ΩP and ΩF match Ω0 for Ω0 < 800 RPM. However, for Ω0≥800, both ΩP and ΩF undershoot Ω0. This was expected since the machine for these spindle speeds was operating at its high range and, therefore, spindle slipping and inaccuracies in the Ω0 were expected. The registered ΩP and ΩF, which are considered more accurate than the nominal value set on the machine, confirm this expectation.

    Figure 12.

    Figure 12. (a) Comparison of ΩP, ΩF, ΩW and the nominal spindle speeds using experimental data. The inset panels show that both ΩP and ΩF are close to each other and to the nominal speed even at the sections where they deviate the most. By contrast, ΩW shows larger deviations for Ω0 > 98 RPM. The figure shows that both ΩP and ΩF give a speed that is lower than its nominal value for Ω0 > 800 RPM. This is attributed to the chuck slip effects at the higher range of speeds. (b) Average runtime for the experimental data on a semi-log scale where three trials were used for each RPM (table 2). (Online version in colour.)

    Figure 12b shows the runtime associated with calculating ΩP, ΩF and ΩW from the experimental data on a semi-log scale. Note that ΩF was computed using the optimized FFT function from Python's scipy.fftpack package while no optimizations were attempted for the ΩP and ΩW algorithms. Therefore, the runtime for the Walsh–Hadamard approach, which is the worst in this case, is shown for completeness, but until optimized algorithms are tested no concrete conclusions about its runtime can be made. On the other hand, the runtime for both Fourier and persistence is small and comparable (note the y-scaling), although persistence is consistently faster than Fourier.

    6. Conclusion and discussion

    This paper described a new approach for periodic pulse detection using tools from TDA, specifically, zero-dimensional persistent homology. This viewpoint allows us to provide guarantees via theorem 3.1 for counting pulses even in the presence of digital ringing and stochasticity in period length. The approach involves computing the persistence diagram for the set of times when the time series is above a threshold, and using the widest split in this diagram to determine how many true pulses were seen. The method is further extended to estimate RPM.

    The described approach was verified using numerical and experimental studies. Simulated pulse trains with small duty cycles were generated using two independent uniform noise components: one that varied the spacing between the pulses ϵ, and one that introduced spurious peaks near the rising and falling edges of the signal α. The α noise component simulates the digital ringing observed in the experimental data. Figure 6a,b shows that, while both Fourier and persistence methods are resilient to ϵ < 7%, persistence is better suited for higher ϵ values. Specifically, for the Fourier-based approach the relative error progressively grows to about 40% across four horizontal bands of ϵ intervals while the relative error of the persistence-based approach remains below 10%.

    Figure 6c,d shows that the persistence-based approach gives results with a relative error below 7% for α < 30%. The Fourier-based approach gives errors as high as 24% for the same α range. Further, the figure shows that, while the persistence-based approach gives a relative error between 15% and 25% for α > 30%, the error in the Fourier-based approach is at least 40%. Beyond 42%, the described approach can regain its accuracy by detecting high logic between valleys as opposed to detecting low logic between peaks.

    A close examination of sample trajectories of ϵ versus ΩP and ΩF in figure 8 shows that both methods have tight error bounds. However, the persistence-based method remains close to Ω0 while ΩF deviates rapidly from Ω0 for ϵ > 5%. Similarly, figure 7 shows that persistence remains close to Ω0 with noticeable deviations from Ω0 starting at α > 30%. The same figure shows that the Fourier-based approach deviates from Ω0 over the whole 0.02 < α < 0.5 range.

    Computationally, persistence and Fourier have the same theoretical runtime. However, the benchmark in figure 10 shows that, while both runtimes are close, persistence is consistently faster. Further, the figure shows larger variations in the runtime of Fourier in comparison with persistence, and that these variations become larger as N, the number of grid points in the frequency spectrum, is increased from 1024 to 8192. The 68% confidence band around the two curves is extremely small, which indicates high confidence in the mean as an estimator of the true runtime. A computation of the runtime for the RPM calculation using the experimental data in figure 12b confirms the numerical results and shows that, while both Fourier and persistence are fast, persistence is generally faster. This figure also shows the runtime for the Walsh–Hadamard approach, although our non-optimized code, in this case, is most probably the reason for the discrepancy between Fourier and Walsh. With improved code, we expect the runtime for the frequency and the sequency methods to be comparable with even a better advantage for the latter because it only involves arithmetic operations that do not include multiplication with complex exponentials.

    The experimental investigations depicted in figure 11a further established the reliability of the described approach. Specifically, figure 12 shows that both persistence-based and Fourier-based calculated RPMs are very close. For the specific lathe used in the experiment, both methods gave spindle speeds close to the nominal value for Ω0 < 800 RPM. However, for Ω0 > 800 RPM, both ΩP and ΩF yielded RPM values smaller than the nominal. This discrepancy at higher speeds is attributed to operating the machine at its high range where slipping in the spindle drive can cause the actual RPM to be smaller than the set RPM.

    On the other hand, Walsh-based calculations showed mixed results. Specifically, while ΩW initially matched ΩP and ΩF for Ω0 < 100 RPM, it generally deviated from Ω0 for the rest of the studied range. The deviation was large for certain speeds, as shown by the triangular regions near Ω0∈{175, 320, 770} in figure 12, while for other speeds ΩW would sporadically approach Ω0. We also numerically studied the RPMs computed using the Walsh-based approach relative to the Fourier and persistence methods. These figures were not included in the paper because, although slight improvements were found by using Walsh versus Fourier, the RPMs calculated using the Walsh-based approach closely followed their Fourier counterparts. Consequently, the conclusions regarding the numerical comparison between Fourier and persistence extend to the Walsh-based approach as well, and persistence was found to perform better than the other two methods.

    We conclude that, in general, a sequency-based approach is not reliable for true step detection especially in the presence of digital ringing, and thus we have chosen to keep the article brief by only including Fourier results. Although for deterministic pulses Walsh successfully detects the true RPM, only a small amount of digital ringing was sufficient to deteriorate its performance. The explanation for this is that Walsh–Hadamard does not assume periodicity but rather gives a measure of the zero crossings (or alternatively the number of humps) of the Walsh basis functions. When the signal is free of digital ringing, the number of humps in the prominent Walsh basis gives an accurate representation of the underlying RPM. However, when digital ringing is introduced, the sequency spectrum fills up with ‘spurious’ sequencies that lead to erroneous results. By contrast, the Fourier approach assumes periodicity and, therefore, shows more resilience in the face of digital ringing. This is because Fourier breaks the signal into the summation of periodic functions with the main signal component having a large contribution in the power spectrum, while the digital ringing adds power mostly at the higher end of the spectrum. Nevertheless, both the sequency and the frequency methods require setting a threshold parameter to pick the frequency/sequency that represents the most prominent component above the noise floor in the spectrum. This makes the presented persistence approach more attractive because in addition to its demonstrated robustness and speed it does not require choosing any threshold parameter, as long as the mild assumptions of theorem 3.1 are satisfied.

    One caveat we mention here is that it is important to pick a span of time that contains a reasonable number of true peaks. If ϵ is too large and the time span is too long, the resulting distribution of the gaps may fill up the histogram of the points in the persistence diagram (figure 3), which can make distinguishing the true peaks difficult or impossible when using the accordion model; see equation (3.9).

    (a) Higher dimensional extensions

    In essence, we have taken a difficult problem—clustering points in Rd—to a rather simplified setting—d = 1. This restriction allows for very fast algorithms, particularly owing to the fact that one knows the minimal spanning tree in advance. Still, algorithms for computing zero-dimensional persistent homology for points in higher dimensions are still quite fast: an implementation using a modified union-find algorithm runs in O((n)), where n is the number of points and α is the notoriously slow-growing inverse Ackermann function.

    One could, however, imagine higher dimensional signals (i.e. images) where there are sufficient guarantees on the pulse rate and noise inputs so as to create higher dimensional analogues of theorem 3.1. Consider, for example, the image of figure 13 generated by Y (s, t) = X1(s) · X2(t) for two realizations of equation (3.1). The zero-dimensional persistence diagram is shown in figure 13b, and, as with the previous data shown in this paper, there is a clear distinction between the points. Since there are 35 points above a threshold of 8 in the histogram, we can conclude that there are 36 clusters in the image, which can be confirmed by inspection. Thus, we hypothesize that this method is quite useful in higher dimensional applications, although such an application is still an open question for future research.

    Figure 13.

    Figure 13. An example of a higher dimensional analogue of the given method for pulse counting. The image (a) is a product of two realizations of equation (3.1) and the histogram (b) shows the resulting zero-dimensional persistence diagram. The distinct split in the diagram can be used to count the number of clusters seen in (a). (Online version in colour.)

    Data accessibility

    The data files, code and experimental parameters can be found at the gitlab repository https://gitlab.msu.edu/TSAwithTDA/TDA-for-true-step-detection-in-PWC-signals.

    Author's contributions

    F.A.K. and E.M. wrote the paper. F.A.K. collected and analysed the synthetic and the experimental data. E.M. wrote the code for generating the synthetic data and derived the mathematical theory.

    Competing interests

    We declare we have no competing interests.

    Funding

    This material is based upon work supported by the National Science Foundation under grant nos. CMMI-1759823 and DMS-1759824 with PI F.A.K., and CMMI-1800466 and DMS-1800446 with PI E.M.

    Disclaimer

    The views and opinions expressed in this article are those of the authors.

    Acknowledgments

    The authors acknowledge the help of David Petrushenko in setting up the lathe for the experiments. They also thank the four anonymous reviewers whose feedback improved the paper, and Basil Nabulsi for helpful discussions.

    Footnotes

    1 This is often colloquially shortened to ‘persistence’.

    2 We work with field coefficients, most typically Z2, so the homology group is, in fact, a vector space.

    3 Experts will note that, in general, persistence diagrams are given as collections of points in R2. Since for zero-dimensional persistence on point clouds, all components are born at 0, the true persistence diagram would be the collection {(0, d1), …, (0, dk)}. We elect to drop the repeated first coordinate in order to simplify the description.

    4 The distance functions fAandfBmust be what is known as tame in the TDA literature. As all sets considered in this paper have finitely many connected components inR, we will always satisfy this assumption.

    5 This is sometimes called a ‘Euclidean spanning tree’.

    Published by the Royal Society. All rights reserved.

    References

    • 1
      Sowa Y, Rowe AD, Leake MC, Yakushi T, Homma M, Ishijima A, Berry RM. 2005Direct observation of steps in rotation of the bacterial flagellar motor. Nature 437, 916–919. (doi:10.1038/nature04003) Crossref, PubMed, ISIGoogle Scholar
    • 2
      Little MA, Jones NS. 2010Sparse Bayesian step-filtering for high-throughput analysis of molecular machine dynamics. In Proc. 2010 IEEE Int. Conf. on Acoustics, Speech and Signal Processing, Dallas, TX, 14–19 March 2010, pp. 4162–4165. New York, NY: IEEE. Google Scholar
    • 3
      Little MA, Jones NS. 2011Generalized methods and solvers for noise removal from piecewise constant signals. I. Background theory. Proc. R. Soc. A 467, 3088–3114. (doi:10.1098/rspa.2010.0671) LinkGoogle Scholar
    • 4
      Little MA, Jones NS. 2011Generalized methods and solvers for noise removal from piecewise constant signals. II. New methods. Proc. R. Soc. A 467, 3115–3140. (doi:10.1098/rspa.2010.0674) LinkGoogle Scholar
    • 5
      Nirody JA, Sun Y-R, Lo C-J. 2017The biophysicist's guide to the bacterial flagellar motor. Adv. Phys. X 2, 324–343. (doi:10.1080/23746149.2017.1289120) ISIGoogle Scholar
    • 6
      Youssef A, Matthews D, Guzzomi A, Pan J. 2017Measurement of pressure fluctuations inside a model thrust bearing using PVDF sensors. Sensors 17, 878. (doi:10.3390/s17040878) Crossref, ISIGoogle Scholar
    • 7
      Jin Y, Ji S, Liu B, Chamorro L. 2016On the role of thickness ratio and location of axis of rotation in the flat plate motions. J. Fluids Struct. 64, 127–137. (doi:10.1016/j.jfluidstructs.2016.05.004) Crossref, ISIGoogle Scholar
    • 8
      Albers A, Schille F, Behrendt M. 2016Method for measuring and analyzing transient powertrain vibrations of hybrid electric vehicles on an acoustic roller test bench. In Proc. 9th Int. Styrian Noise, Vibration & Harshness Congress: The European Automotive Noise Conference, Graz, Austria, 22–24 June 2016. Warrendale, PA: SAE International. Google Scholar
    • 9
      Rahman M, Salyers T, Ahmed M, ElShahat A, Soloiu V, Maroha E. 2016Investigation of aerodynamic performance of helical shape vertical-axis wind turbine models with various number of blades using wind tunnel testing and computational fluid dynamics. In Proc. ASME International Mechanical Engineering Congress and Exposition, Tampa, FL, 3–9 November 2017, vol. 7: Fluids engineering. New York, NY: ASME. Google Scholar
    • 10
      Urbikain G, Olvera D, de Lacalle LL, Elías-Zúñiga A. 2016Spindle speed variation technique in turning operations: modeling and real implementation. J. Sound Vib. 383, 384–396. (doi:10.1016/j.jsv.2016.07.033) Crossref, ISIGoogle Scholar
    • 11
      Amini S, Aghaei M, Lotfi M, Hakimi E. 2017Analysis of linear vibration in rotary turning of AISI 4140 steel. Int. J. Adv. Manuf. Technol. 91, 4107–4116. (doi:10.1007/s00170-017-0108-5) Crossref, ISIGoogle Scholar
    • 12
      Das B, Bag S, Pal S. 2016Defect detection in friction stir welding process through characterization of signals by fractal dimension. Manuf. Lett. 7, 6–10. (doi:10.1016/j.mfglet.2015.11.006) CrossrefGoogle Scholar
    • 13
      Wan M, Ma Y-C, Feng J, Zhang W-H. 2017Mechanics of tapping process with emphasis on measurement of feed error and estimation of its induced indentation forces. Int. J. Mach. Tools Manuf. 114, 8–20. (doi:10.1016/j.ijmachtools.2016.12.003) Crossref, ISIGoogle Scholar
    • 14
      Honeycutt A, Schmitz TL. 2017Surface location error and surface roughness for period-n milling bifurcations. J. Manuf. Sci. Eng. 139, 061010–061018. (doi:10.1115/1.4035371) Crossref, ISIGoogle Scholar
    • 15
      Ester M, Kriegel H-P, Sander J, Xu X. 1996A density-based algorithm for discovering clusters in large spatial databases with noise, pp. 226–231. Portland, OR: AAAI Press. Google Scholar
    • 16
      Ghrist R. 2014Elementary applied topology. Scotts Valley, CA: CreateSpace Independent Publishing Platform. Google Scholar
    • 17
      Edelsbrunner H, Harer J. 2010Computational topology: an introduction. Providence, RI: American Mathematical Society. Google Scholar
    • 18
      Much E. 2017A user's guide to topological data analysis. J. Learn. Anal. 4, 47–61. (doi:10.18608/jla.2017.42.6) CrossrefGoogle Scholar
    • 19
      Kaczynski T, Mischaikow K, Mrozek M. 2004Computational homology. Berlin, Germany: Springer. CrossrefGoogle Scholar
    • 20
      Zomorodian A. 2005Topology for computing. Cambridge Monographs on Applied and Computational Mathematics, vol. 16. Cambridge, UK: Cambridge University Press. CrossrefGoogle Scholar
    • 21
      Ghrist R. 2008Barcodes: the persistent topology of data. Bull. Am. Math. Soc. 45, 61–75. (doi:10.1090/s0273-0979-07-01191-3) Crossref, ISIGoogle Scholar
    • 22
      Carlsson G. 2009Topology and data. Bull. Am. Math. Soc. 46, 255–308. (doi:10.1090/S0273-0979-09-01249-X) Crossref, ISIGoogle Scholar
    • 23
      Robins V, Meiss JD, Bradley E. 1998Computing connectedness: an exercise in computational topology. Nonlinearity 11, 913–922. (doi:10.1088/0951-7715/11/4/009) Crossref, ISIGoogle Scholar
    • 24
      Robins V, Meiss J, Bradley E. 2000Computing connectedness: disconnectedness and discreteness. Physica D 139, 276–300. (doi:10.1016/S0167-2789(99)00228-6) Crossref, ISIGoogle Scholar
    • 25
      Robins V, Rooney N, Bradley E. 2004Topology-based signal separation. Chaos 14, 305–316. (doi:10.1063/1.1705852) Crossref, PubMed, ISIGoogle Scholar
    • 26
      Robinson M. 2014Topological signal processing. Berlin, Germany: Springer. CrossrefGoogle Scholar
    • 27
      Edelsbrunner H, Letscher D, Zomorodian A. 2002Topological persistence and simplification. Discrete Comput. Geom. 28, 511–533. (doi:10.1007/s00454-002-2885-2) Crossref, ISIGoogle Scholar
    • 28
      Zomorodian A, Carlsson G. 2004Computing persistent homology. Discrete Comput. Geom. 33, 249–274. (doi:10.1007/s00454-004-1146-y) Crossref, ISIGoogle Scholar
    • 29
      Garland J, Bradley E, Meiss JD. 2016Exploring the topology of dynamical reconstructions. Physica D 334, 49–59. (doi:10.1016/j.physd.2016.03.006) Crossref, ISIGoogle Scholar
    • 30
      MacPherson R, Schweinhart B. 2012Measuring shape with topology. Russ. J. Math. Phys. 53, 073516. (doi:10.1063/1.4737391) Crossref, ISIGoogle Scholar
    • 31
      Alexander Z, Bradley E, Meiss JD, Sanderson NF. 2015Simplicial multivalued maps and the witness complex for dynamical analysis of time series. SIAM J. Appl. Dyn. Syst. 14, 1278–1307. (doi:10.1137/140971415) Crossref, ISIGoogle Scholar
    • 32
      Maletić S, Zhao Y, Rajković M. 2016Persistent topological features of dynamical systems. Chaos 26, 053105. (doi:10.1063/1.4949472) Crossref, PubMed, ISIGoogle Scholar
    • 33
      Pereira CM, de Mello RF. 2015Persistent homology for time series and spatial data clustering. Expert Syst. Appl. 42, 6026–6038. (doi:10.1016/j.eswa.2015.04.010) Crossref, ISIGoogle Scholar
    • 34
      Perea JA. 2016Persistent homology of toroidal sliding window embeddings. In Proc. 2016 IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP), Shanghai, China, 20–25 March 2016. New York, NY: IEEE. Google Scholar
    • 35
      Perea JA, Harer J. 2015Sliding windows and persistence: an application of topological methods to signal analysis. Found. Comput. Math. 15, 799–838. (doi:10.1007/s10208-014-9206-z) Crossref, ISIGoogle Scholar
    • 36
      Tralie CJ, Perea JA. 2018(Quasi)Periodicity quantification in video data, using topology. SIAM J. Imaging Sci. 11, 1049–1077. (doi:10.1137/17M1150736) Crossref, ISIGoogle Scholar
    • 37
      Robinson M. 2015Universal factorizations of quasiperiodic functions. In Proc. 2015 Int. Conf. on Sampling Theory and Applications (SampTA), Washington, DC, 25–29 May 2015, pp. 588–592. New York, NY: IEEE. Google Scholar
    • 38
      de Silva V, Skraba P, Vejdemo-Johansson M. 2012Topological analysis of recurrent systems. In Proc. NIPS 2012 Workshop on Algebraic Topology and Machine Learning, Lake Tahoe, NV, 3–8 December 2018. Red Hook, NY: Curran Associates Inc. Google Scholar
    • 39
      Emrani S, Gentimis T, Krim H. 2014Persistent homology of delay embeddings and its application to wheeze detection. IEEE Signal Process. Lett. 21, 459–463. (doi:10.1109/LSP.2014.2305700) Crossref, ISIGoogle Scholar
    • 40
      Alexander Z, Meiss JD, Bradley E, Garland J. 2012Iterated function system models in data analysis: detection and separation. Chaos 22, 023103. (doi:10.1063/1.3701728) Crossref, PubMed, ISIGoogle Scholar
    • 41
      Gidea M, Katz Y. 2018Topological data analysis of financial time series: landscapes of crashes. Physica A 491, 820–834. (doi:10.1063/1.3701728) Crossref, ISIGoogle Scholar
    • 42
      Gidea M. 2017Topological data analysis of critical transitions in financial networks, pp. 47–59.Cham, Switzerland: Springer International Publishing, Google Scholar
    • 43
      Robinson M. 2013Multipath-dominant, pulsed Doppler analysis of rotating blades. IET Radar Sonar Navig. 7, 217–224. (doi:10.1049/iet-rsn.2012.0117) Crossref, ISIGoogle Scholar
    • 44
      Vejdemo-Johansson M, Pokorny F, Skraba P, Kragic D. 2015Cohomological learning of periodic motion. In Applicable algebra in engineering, communication and computing, pp. 1–22. Heidelberg, Germany: Springer. Google Scholar
    • 45
      Bauer U, Munk A, Sieling H, Wardetzky M. 2015Persistence barcodes versus Kolmogorov signatures: detecting modes of one-dimensional signals. Found. Comput. Math. 17, 1–33. (doi:10.1007/s10208-015-9281-9) Crossref, ISIGoogle Scholar
    • 46
      Berwald JJ, Gidea M, Vejdemo-Johansson M. 2014Automatic recognition and tagging of topologically different regimes in dynamical systems. Discontinuity Nonlinearity Complexity 3, 413–426. (doi:10.5890/DNC.2014.12.004) CrossrefGoogle Scholar
    • 47
      Khasawneh FA, Munch E. 2016Chatter detection in turning using persistent homology. Mech. Syst. Signal Process. 70–71, 527–541. (doi:10.1016/j.ymssp.2015.09.046) Crossref, ISIGoogle Scholar
    • 48
      Khasawneh FA, Munch E. 2014Exploring equilibria in stochastic delay differential equations using persistent homology. In Proc. of the ASME 2014 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Buffalo, NY, 17–20 August 2014. Paper no. DETC2014/VIB-35655. New York, NY: American Society of Mechanical Engineers. Google Scholar
    • 49
      Khasawneh FA, Munch E. 2014Stability determination in turning using persistent homology and time series analysis. In Proc. of the ASME 2014 Int. Mechanical Engineering Congress & Exposition, Montreal, Canada, 14–20 November 2014. Paper no. IMECE2014-40221. New York, NY: American Society of Mechanical Engineers. Google Scholar
    • 50
      Khasawneh FA, Munch E. 2017Utilizing topological data analysis for studying signals of time-delay systems, pp. 93–106. Cham, Switzerland: Springer International Publishing. Google Scholar
    • 51
      Deckard A, Anafi RC, Hogenesch JB, Haase SB, Harer J. 2013Design and analysis of large-scale biological rhythm studies: a comparison of algorithms for detecting periodic signals in biological data. Bioinformatics 29, 3174–3180. (doi:10.1093/bioinformatics/btt541) Crossref, PubMed, ISIGoogle Scholar
    • 52
      Perea JA, Deckard A, Haase SB, Harer J. 2015SW1PerS: sliding windows and 1-persistence scoring; discovering periodicity in gene expression time series data. BMC Bioinf. 16, 257. (doi:10.1186/s12859-015-0645-6) Crossref, PubMed, ISIGoogle Scholar
    • 53
      Cohen-Steiner D, Edelsbrunner H, Harer J. 2007Stability of persistence diagrams. Discrete Comput. Geom. 37, 103–120. (doi:10.1007/s00454-006-1276-5) Crossref, ISIGoogle Scholar
    • 54
      Cohen-Steiner D, Edelsbrunner H, Harer J, Mileyko Y. 2010Lipschitz functions have lp-stable persistence. Found. Comput. Math. 10, 127–139. (doi:10.1007/s10208-010-9060-6) Crossref, ISIGoogle Scholar
    • 55
      Otter N, Porter MA, Tillmann U, Grindrod P, Harrington HA. 2017A roadmap for the computation of persistent homology. EPJ Data Sci. 6, 17. (doi:10.1140/epjds/s13688-017-0109-5) Crossref, PubMed, ISIGoogle Scholar
    • 56
      Stéphane M. 2009Chapter 2—the Fourier kingdom. In A wavelet tour of signal processing (ed. M Stéphane), 3rd edn, pp. 33–57. Boston, MA: Academic Press. Google Scholar
    • 57
      Broadbent HA, Maksik YA. 1992Analysis of periodic data using Walsh functions. Behav. Res. Methods Instrum. Comp. 24, 238–247. (doi:10.3758/BF03203501) CrossrefGoogle Scholar
    • 58
      Hatcher A. 2002Algebraic iopology. Cambridge, UK: Cambridge University Press. Google Scholar
    • 59
      Munkres J. 2000Topology, 2nd edn. Upper Saddle River, NJ: Pearson. Google Scholar
    • 60
      Murtagh F, Contreras P. 2011Algorithms for hierarchical clustering: an overview. In Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, vol. 2, pp. 86–97. Hoboken, NJ: Wiley. Google Scholar
    • 61
      Beauchamp KG. 1984Applications of Walsh and related functions: with an introduction to sequence theory. Cambridge, MA: Academic Press. Google Scholar
    • 62
      Ashrafi A. 2017Walsh–Hadamard transforms: a review. In Advances in imaging and electron physics, pp. 1–55. San Diego, CA: Elsevier. Google Scholar
    • 63
      Sexton B, Stone JS. 1978The stability of machining with continuously varying spindle speed. Ann. CIRP 27, 321–326. Google Scholar
    • 64
      Yilmaz A, AL-Regib E, Ni J. 2002Machine tool chatter suppression by multi-level random spindle speed variation. J. Manuf. Sci. Eng. 124, 208–216. (doi:10.1115/1.1378794) Crossref, ISIGoogle Scholar
    • 65
      Lin S, DeVor R, Kapoor S. 1990The effects of variable speed cutting on vibration control in face milling. J. Eng. Ind. 112, 1–11. (doi:10.1115/1.2899290) CrossrefGoogle Scholar
    • 66
      Radulescu R, Kapoor SG, DeVor RE. 1997An investigation of variable spindle speed face milling for tool-work structures with complex dynamics. Part 1: Simulation results. J. Manuf. Sci. Eng. 119, 266–272. (doi:10.1115/1.2831103) Crossref, ISIGoogle Scholar