Extract rpm from tachometer signal

Hi guys… I need to extract a tachometer signal (pulses) and convert it to rpm. Would anyone know how to help write this script?

In this case, the script is used within the Sigview software, which allows the use of the Julia language.

Thanks

What is the format of the signal?

Paulo… the signal is pulse type (step) generated by an optical sensor and reflective tape (0 - 5 V) .

Follow an example of the signal type.

The teensyduino library has some nice algorithms for this.

First of all you need to find the tach crossing times or indices. It appears that you have a digital signal. If amp is your signal then one way is:

trigamp = 2.5
tachcrossindice=findall((amp[2:end] .>= trigamp) .& (amp[1:end-1] .< trigamp))

It may well be that a for loop is faster. If your signal has been filtered then use linear interpolation to obtain a more accurate tach crossing time. You can adapt this for negative going pulses or if your duty cycle is 0.5, for both positive and negative going pulses.

Then you need to get the speed from the times. You have two choices, pick a time and count the number of pulses, or count the time between pulses. I use the later. I then normally smooth the rotational speed with the package “SmoothingSplines”, which seems to work well.

I just watched the Julia scripting video on https://sigview.com/. How do you like the Sigview software?

Sigview for me is very practical for some applications. I use National Instruments boards that it has driver support. However, some functions need to be done via Julia Script.

Take a look at DSP.jl:

julia> using DSP, Unitful

julia> function squarewave(rate, duration, amplitude, freq)
           tt = range(0.0, duration, length = rate * duration)
           s =  [amplitude * (sin(2π*freq*t) > 0) for t in tt]
           return s, tt
       end
squarewave (generic function with 1 method)

julia> s, tt = squarewave(4096, 10, 4.98, 75.5)
([0.0, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98, 4.98  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 0.0:0.00024414658561000026:10.0)

julia> pg = periodogram(s, fs=4096);

julia> pg.power[1] = 0
0

julia> pg.freq[argmax(pg.power)]
75.5

julia> ans*u"Hz" |> u"minute^-1"
4530.0 minute^-1
1 Like

@stillyslalom, very nice code. Just a thought about the square wave, which could also be written more compactly as:

squarewave(t, A₀, freq) = A₀ * sign(sin(2π*freq*t))

and then broadcasted over the input time range t

While this method works well for random signals, for very narrow band signals it might have some issues if the starting and end point of data acquisition is not very well synchronized with the signal itself. Spectral leakage or something like that. I learned this the hard way while measuring Strouhal number from laminar 2D flow numerical simulations.

For very well behaved signals, just counting when the signal passes a threshold is much more accurate.

@Paulo_Jabardo The tachometer signal is very uniform. But I am not able to make the script to perform the count when the value crosses the limit. Below is an example of how

Sigview works:

The signal is the vector input_Windows[1].samples[n]

#Example (multiplying all input samples with 2)

function process( input_windows::Array{SigviewWindow}, output_signal::SigviewSignalWindow )
   output_signal.samples = input_windows[1].samples .* 2
   output_signal.samplingRate = input_windows[1].samplingRate
   return true
end

My head is already melting :slight_smile: :exploding_head:

Thanks for all

If you use the code provided by Jake and assuming that t is a range or array with the times in seconds, then the rotations per minute are simply:

rpm = 60/mean(diff(t[tachcrossindice]))   # 4530.08 in the example

Play around with this function:

"""
`t1, t2, np = pulse_count(fs, x, xmean, up)`

Counts the number of pulses in a well behaved signal. A pulse begins when the signal crosses the value `xmean` either going up (`up=true`) or down (`up=false`).

The function returns where the first pulse starts (`t1`), 
when the last pulse ends (`t2`) and the number of pulses between (`np`).

Arguments:
 * `fs` Sampling rate
 * `x` vector containing the signal
 * `xmean` The trigger level to start counting pulses
 * `up` If true, pulses begin when the signal goes up. 

To calculate the mean pulse frequency, just compute

`(t2 - t1) / npulses`

"""
function pulse_count(fs, x, xmean, up=true)
	n = length(x)
	x0 = x .- xmean # Center the signal

	# Check where x0 changes signal 
	sigchg = x0[1:end-1] .* x0[2:end] .< 0
	# Check whether the signal is going up or down
	if up
		slope = x[2:end] .- x[1:end-1] .> 0
	else
		slope = x[2:end] .- x[1:end-1] .< 0
	end
	pulses = (sigchg .& slope) 
	npulses = sum(pulses)-1
	idx = (1:n-1)[pulses]
	# Interpolate the time where the first pulse starts
	i = idx[1]
	y1 = x0[i]
	y2 = x0[i+1]
	t1 = 1/fs * (i-1 - y1/(y2-y1))
	
	# Interpolate the time where the last pulse ends
	k = idx[end]
	z1 = x0[k]
	z2 = x0[k+1]
	t2 = 1/fs * (k-1 - z1 / (z2-z1))
	return t1, t2, npulses
end

Some nice algorithms here. If you want to see how the speed is changing with time for transient signals, you can find the speed at each pulse. I like working with rotational speed in Hz for easy comparison with spectra in Hz. If your tach pulse signal has been sampled with a DSA there will be a slope at the tach pulse crossing points due to the filtering so you can linearly interpolate the crossing points to obtain more accuracy. Then to get rid of the jitter for sampling at too low a rate you can use SmoothingSplines.jl. Finally if the tach pulse is not nice and clean like yours, there are ways to improve it. If you want to go back to phase, use cumsum to integrate the signal.

using Statistics: mean
using Plots

function squarewave(rate, duration, amplitude, freq)
    tt = range(0.0, duration, length = rate * duration)
    s =  [amplitude * (sin(2π*freq*t) > 0) for t in tt]
    return s, tt
end

amp, t = squarewave(4096, 10, 4.98, 75.5)

fs = 1/t[2]

# put the below into a function
trigamp = mean(extrema(amp))/2.0

tachcrossindice=findall((amp[2:end] .>= trigamp) .& (amp[1:end-1] .< trigamp))
tachcrosstime = tachcrossindice/fs

midtim = (tachcrosstime[2:end] .+ tachcrosstime[1:end-1])/2.0
rotspd = 1.0 ./ diff(tachcrosstime)
rotrpm = 60.0 .* rotspd

plot(tim,spd, label=false)

Very nice. Same approach as mine but much simpler. I will rewrite my code.

The thing to be careful about is that, most often, when generating an “artificial” signal, we generate full cycles. And any of the approaches above work well. But when using the output of an instrument, the signal can begin or end at any point of the cycle. In this case using the fft might result in a bias. Of course, if the number of pulses is very large, this is minimized.

On the other hand, if the number of pulses (or cycles in general) is small counting them directly as your approach does is the best approach and interpolation can also help a lot.

Paulo

Guys, I would like to thank you immensely for all the help that is being provided.

All codes sent are very useful and I have learned a lot.

The biggest problem is that it’s very difficult to make some tools in Sigview.

Sigview is a very practical software for some applications, but it still lacks some tools that must be developed via Julia Script.

I will keep trying

Thanks again to all!

Claudio