# Estimate “speed” of repeating signal pattern

I have a digitized signal containing a repeating pattern. One could call it periodic, but the time for one period is not constant. The "speed" with which the period is advanced can change over time. How can i estimate/calculate this "speed" function if i know how one period with constant speed looks like?

An example: Take one period of a sine wave. I know how this period of a sine wave looks like when plotted with a constant frequency. Now i record a signal where this sine wave changes frequency constantly. How can i calculate the frequency for each point?

My concrete example: The signal represents the IR light reflected by the metallic plate in an energy meter (Ferraris type). There is a red colored area on this plate which i can see as a dip in signal intensity. One rotation of the plate is 1/75 kWh consumed. The interesting thing: My IR detector is so sensitive that i can see more patterns on this metal plate than just the dip from the red area. The "speed" at which this pattern evolves is basically the current power draw, which i want to calculate. If i just count the passing of the red area, i get only a very course resolution (one minute average) of the power consumption, but if a manage to look at the whole period in total i think i can get down to second-resolution. I hope this makes some sense.

• Are you restricting yourself to sines or to generally periodic signals? The answer to the first case is significantly easier, but might not help you if you're not dealing with a sine. – Marcus Müller Nov 16 '19 at 0:45
• The signal i want to analyze is not sinusoidal. This was just an example. :-) – Lazarus535 Nov 16 '19 at 0:51
• It would be useful to detail how the pattern repeats: with varying amplitude, allowed percentage of possible truncation, varying in scale, etc. Writing a "model" can help a lot – Laurent Duval Nov 16 '19 at 0:55
• The pattern varies for 30s to 2min length and always has the same amplitude and scale. It is recorded with a fairly high signal-to-noise ratio. What you mean by truncation i don't know. – Lazarus535 Nov 16 '19 at 1:01
• Consider "scale" as a stretch. If the frequency changes, then periodicity per se may loose meaning. By truncation, I suspected that the pattern might repeat "not fully". Here, you seem to have a kind of basic rotation? – Laurent Duval Nov 16 '19 at 1:21

You may be able to use a modified transform to compute change in a time stretch parameter (I will use $$\zeta$$) versus time (t). Your success in being able to do this depends on the cross correlation properties of your pattern with time stretched versions of itself. Let me explain:

First notice that the Fourier Transform is a correlation over all possible frequencies (here using the simpler form with complex tones $$e^{j\omega t}$$, once you know that $$e^{j\omega t}$$ represents a single tone in frequency this form is so much cleaner than expressing it with sines and cosines; since a cosine is two complex tones: $$2cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$$).

$$F(j\omega) = \int_{t=-\infty}^\infty x(t) e^{-j\omega t}dt$$

Observe that in this process, you set $$\omega$$ to every possible value, and for each $$\omega$$ you compute the correlation of your waveform with that frequency, where what I refer to as "correlation" generally is a complex conjugate multiplication and integration process (and digitally with the DFT we do a multiply and accumulate over a finite length):

$$Corr = \int x(t) y^*(t)$$

Similarly, the short term FT (STFT) we introduce a time variable $$\tau$$ where we sweep a window function (which selects a portion of our waveform) to access the change in frequency versus $$\tau$$:

$$F(\tau, j\omega) = \int_{t=-\infty}^\infty x(t)w(t-\tau) e^{-j\omega t}dt$$

Where $$w(t)$$ is a window function, and we have the trade space of time resolution and frequency resolution to contend with. (The simplest window is a "rectangular window" which is all one over a specified duration of time). If our signal is changing very slowly with time, we can use longer windows in the STFT to achieve fine frequency resolution. but if it is changing quickly, and we want to observe that then we must use shorter windows resulting in courser frequency resolution. Within the choice of window itself we have the trade space of resolution versus sensitivity (dynamic range): a rectangular window has the finest frequency resolution but very poor sensitivity: the sidebands roll-off at 1/f. Higher performance windows (Kaiser, Hamming, Blackman etc) can have much better sideband rejection (better sensitivity) but always at the expense of the main lobe width (resolution).

To the degree that your base pattern, what I will call m(t), is sufficiently uncorrelated to time stretched versions of itself ($$m(t/\zeta)$$), you can modify this transform to be what you are looking for. "Sufficiently" here defines the sensitivity range of your approach. Assuming it is sensitive enough for the actual time stretching you are experiencing, you can then use this appraoch to determine how your actual signal $$x(t)$$ correlates to to different time stretching of your base pattern. You will have the same trade space as in the STFT: in this case $$\zeta$$ resolution versus time resolution. If the $$\zeta$$ is changing very slowing then you can integrate over a long duration and get a fine resolution for $$\zeta$$ up to the maximum cross correlation properties your base function allows, but if changing quickly and you want to measure that fast change then you will need to use shorter integration lengths.

That said if you want to go down this path I suggest:

Step 1) Analyze the cross correlation property of your signal $$m(t)$$ with time stretched versions of itself using:

$$rr_{t \zeta} = \int_0^T m^*(t)m(t/\zeta)d(t)$$

Where $$\zeta$$ is the parameter stretching your time, when $$\zeta=1$$ you get the autocorrelation at time=0, for values of $$\zeta$$ less than 1, your waveform $$m(t/\zeta)$$ will be stretched in time relative to $$m(t)$$, and for values of $$\zeta$$ greater than 1, the waveform will be compressed in time. If your waveform is real then you can ignore the complex conjugate, however this is critical to include for complex waveforms, so I included it in the formula.

The amount of cross correlation will indicate the usability of this approach and your measurement sensitivity and precision. Check this with increasing values for T (integration time) to access usability range for T (which would then provide guidance on window lengths). Assuming the cross-correlation is low enough for your use, proceed to step 2.

Step 2) Use the modified transform to determine $$\zeta$$ versus time for your actual capture waveform:

$$F(\tau, \zeta) = \int_{t=-\infty}^{\infty} x(t)w(t-\tau)m^*(t/\zeta)dt$$

Where: $$x(t)$$ is your received signal

$$w(t)$$ is a window function (rectangular, Kaiser, Hamming etc)

$$\tau$$ is your time parameter for the rate of change in the stretching of your base pattern

$$\zeta$$ is the time stretch parameter for your base pattern.

The result is a 3 dimensional plot with $$\zeta$$ on the horizontal axis, $$\tau$$ on the vertical axis, and the magnitude result for your correlation along the z axis. Where the magnitude is maximum indicates the best estimate for $$\zeta$$ at that time $$\tau$$.

• Thanks for the answer! In the last paragraph, do you mean 2D (contour) plot, instead of 3D plot? I will try and implement your suggestion. :-) – Lazarus535 Nov 16 '19 at 17:45
• Sure yes a 2D contour plot would be fine-- I did mean a 3d plot as there is an x, y and z axis (magnitude at each t, \tau location but a contour plot would be fine. Do the cross correlation verification first and please post your results of that as it may illuminate other solutions or limitations. – Dan Boschen Nov 16 '19 at 17:47
• Thanks for clarifying! I will post my findings once i get to it. Still might be some time tough... – Lazarus535 Nov 16 '19 at 18:02

The principal idea behind my recommended approach is to identify characteristic points in your cycle. You can then make a mapping from time within the cycle (as it is occuring) and the location of that point within your representative cycle. Once you have a set of mapping points, create a best fit function. This will then give you a mapping function of real time to representative time. Since you already know the relationship between representative time and current draw, you can then calculate the relationship between real time and current draw quite accurately.

So, how to find the characteristic spots most easily?

I would do a smoothing/differenting pass on your data, for instance:

$$y[n] = x[n] - 2 *x[n-d] + x[n-2d]$$

The bigger the $$d$$, the smoother it is. If you are more comfortable without a phase shift:

$$y[n] = x[n+d] - x[n+1] - x[n-1] + x[n-d]$$

This is possible with post processing.

The zero crossings of this will then make good characteristic values.

For finer resolution, if your signal is bumpy enough, you can repeat the process for a second differentiation and use its zero crossings also.

On the other hand, if your signal (or the differentiated signal) is fairly smooth, you can partition it into rising and falling pieces and define an inverse function on each piece. More characteristic value mappings can then be generated by finding the time location of various heights along each interval.