Safety check scheme of piezo control waveform

[kdw503]

Last time, we decided to change our safety check method from checking control speed to checking resonance frequency components.
I added some codes analyzing frequency components of control pulses using FFTW library.
Also, I added a function which integrates frequency components around(bandwidth) the resonance frequency(resonancefreq), compares the value with a certain threshold and issues error message.

Here, we need to decide those values : resonancefreq, bandwidth, threshold.

And, Imagine gets control input from external .json control file but can generate those controls internally using GUI setup now.
Then, does this means Imagine needs to equip notch filter to make piezo pulse be smoother?
If so, should it be IIR or FIR filter? Any good library related this filter?
(Until I found is https://github.com/berndporr/iir1)
Any idea would be appreciated.

ISSUE

1. Decide resonancefreq, bandwidth and threshold value for piezo safety check.
2. If fail, should we provide the user an option to filter out those coompoents in his/her piezo control waveform?
3. Does Imagine needs notch filter?

[timholy]

@HaoyangRong, what's the result of your analysis? We could use your help on those 3 parameters.

I think all Imagine should do is to check and then reject; it shouldn't try to fix the problem, the user should do that her/himself.

[kdw503]

I changed the method for checking resonance frequency component as below.

ratio = (power around(bandwidth) resonance frequency(resonanceFreq))/(total power)
if (ratio > threshold)  then error

And, I just hard-coded some parameters roughly with these.

bandwidth = 100hz
threshold = 1% (triangular waveform generated by Imagine typically has around 0.1% when resonanceFreq = 5kHz)

But, resonanceFreq should be recalculated with below equation whenever the load attached to the actuator is changed. And, this value could be effective by just updating f_res field in imagine.js file and re-executing Imagine without compile.

resonanceFreq = f0res * sqrt(meff/(meff+M))

f0res : Resonant frequency of the actuator used in ocpi-2
meff : Effective mass of the actuator used in ocpi-2
M : Mass attached to the actuator (ex. objective lens...)

[Cody-G]

bandwidth = 100hz

Why such a large bandwidth? If the resonance frequency is 100Hz then this means you're calculating the power between 50 and 150Hz? I think the range should be smaller. What about a bandwidth of 30Hz? 30Hz is still arbitrary, but I think 100Hz is too conservative.

when resonanceFreq = 5kHz

Is 5kHz the number you're getting? That seems much too high. That may explain why you chose the 100Hz bandwidth I just mentioned.

Also I think we should do the power calculations in a "windowed" manner. If, for example, the piezo command is a slow sinusoid for 99% of the recording period but there are a few moments with abrupt (dangerous) steps in voltage then we should still consider the command to be unsafe. It's quite likely that someone will try to run such a waveform by mistake due to a bug in their script. The precise duration of the moving window probably doesn't matter that much but I think it should be fairly short, like a couple of seconds.

[kdw503]

I set the window size as 8192 samples(this can catch 8 pulses of 100Hz signal in 100000hz sample rate). And, I didn't move the window one sample by one sample because of the long computation time. I moved the window with a step size which is the half of the window size. I thought that this step size will works ok if the total sample size is quite long enough.
About bandwidth, you are right. I just roughly decided when resonanceFreq = 5kHz.
Then let's set this as 30Hz temporarily.
About dangerous steps, we don't have speed limit any more. This step could be catched as high frequency components in frequency analysis but this is not quite related with resonance frequency.
That means we don't have any screening method for the step pulses right now. Any idea about this?
@Cody-G Could you let me know these parameters below?

f0res : Resonant frequency of the actuator used in ocpi-2
meff : Effective mass of the actuator used in ocpi-2
M : Mass attached to the actuator in ocpi-2 currently

[Cody-G]

And, I didn't move the window one sample by one sample because of the long computation time. I moved the window with a step size which is the half of the window size.

Sounds good. I don't see any advantage to stepping sample-by-sample.

One more thought: In the case that the safety check fails could we print out some information that helps the user find the problematic part of the command? Maybe we could print the start and stop time of the window that fails the test?

About dangerous steps, we don't have speed limit any more. This step could be catched as high frequency components in frequency analysis but this is not quite related with resonance frequency.
That means we don't have any screening method for the step pulses right now. Any idea about this?

I think we shouldn't allow any commands that exceed the resonant frequency. Voltage steps definitely exceed resonance frequency, so you're right that they should show up in the high frequency parts of the spectrum. What if instead of thinking of this as a bandpass filter we think of it as highpass? Like this:

ratio = (power @ frequencies >= (resonanceFreq-0.5*bandwidth))/(total power)

You may still run into edge cases. For example, if the command is flat throughout the window except that it changes by a single integer value (say from an Int16 value of 33 to 34) then 100% most of the power will show up at high frequencies even though this is a safe waveform. You could only run the check when power is above some minimum value. But I wonder if it will be better to ignore the highest-frequency parts of the spectrum and implement a velocity check separately. I'm not sure, you may just have to test and see.

The piezo resonance calculation can be found here in Equations 10 and 11:
http://www.piezo.ws/piezoelectric_actuator_tutorial/Piezo_Design_part3.php
(These equations together are a bit different than the one you're using. Where did yours come from?)

Here's the data sheet for our piezo, the NanoSX800 CAP.

https://www.piezosystem.com/fileadmin/redakteure/datenblatt/piezo_nanopositionierung/piezo_aktoren_nano_positioniersysteme/nanox_hochleistungs_positionierer/serie_nanosx_800/en/nanoSX_800_ds_Rev05_2015_01_15_MKr.pdf

It provides the stiffness and the "unloaded" resonance frequency (separately for each axis). Take the lowest resonant frequency provided (210g) and use that. I think the number given in the datasheet is not the perfect "unloaded" value because it includes the weight of the hardware attached to the piezo stacks, but hopefully it's close enough for our purposes.

The currently attached mass M (with the small objective lens) for OCPI2 is 264g

[kdw503]

One more thought: In the case that the safety check fails could we print out some information that helps the user find the problematic part of the command? Maybe we could print the start and stop time of the window that fails the test?

Now, we are just showing pop-up error message saying,

anxial piezo control has excessive frequency components (20% > Threshold =1%) around resonance 
frequency (5kHz). Please check 'Frequency component spectrum'

(I added 'Frequency component spectrum' tab in the Imagine.)

What if instead of thinking of this as a bandpass filter we think of it as highpass?

Sound good to me.

if the command is flat throughout the window except that it changes by a single integer value (say from an Int16 value of 33 to 34) then 100% most of the power will show up at high frequencies even though this is a safe waveform.

In this case, it has some DC component about 33.5. So, ratio will not be that big in our current calculation. By the way, I realized that the current calculation method is problematic in some sense for example:

control1 = 10(DC component) + 1(resonance frequency) --> ratio_control1 = 10% --> error!
control2 =  200(DC component) + 1(resonance frequency) --> ratio_control2 = 0.5% --> no error

I think both controls are equally dangerous though, right? We need to remove some components around 0Hz before calculating ratio. Now go back to your example, the ratio of the command might be higher under this new calculation. In this case, we need to run the check only when power(except DC component) is above some minimum value as you said.

But I wonder if it will be better to ignore the highest-frequency parts of the spectrum and implement a velocity check separately.

In this case, we need to set piezo speed limit again. Do you have any value in mind?

Where did yours come from?

From Haoyang. Chapter 5. in https://www.piezosystem.com/piezopedia/piezotheory/
I think yours is same as this except that yours has one more equation for calculating f0res.

It provides the stiffness and the "unloaded" resonance frequency (separately for each axis). Take the lowest resonant frequency provided (210g)

210Hz right? Then, resonanceFreq = 163.8hz because f0res = 210hz, meff=410g and M=264g.
What happen if we need to change objective lens. Those lenses have similar weight?
Thank for your ideas and all the information I needed.

[Cody-G]

Now, we are just showing pop-up error message saying,

I like the error message as-is, but I think it could be even better to include some information about which time window failed the safety check. Especially since, as we mentioned, the problem may not be obvious when looking at the spectrum for the entire command (I assume that's what the new plot window is showing?)

I think both controls are equally dangerous though, right? We need to remove some components around 0Hz before calculating ratio.

Agreed

In this case, we need to set piezo speed limit again. Do you have any value in mind?

I'm still not sure the velocity check will be needed. First I would see if we're getting good results with the spectral method (subtracting DC and computing the ratio only when the total power exceeds a minimum power threshold). Actually now that I think about it the "minimum power threshold" is related to velocity. I think the total power (without DC) within the window should be the square of the total distance traveled by the piezo. If we take the square root of that power and divide that by the duration of the window then we get the average velocity within the window, so at least the velocity check can be performed within the same framework. As for how to set the velocity/power threshold, I think we can use something similar to what we used before? Or if that was too conservative we can raise the limit a bit.

From Haoyang. Chapter 5. in https://www.piezosystem.com/piezopedia/piezotheory/
I think yours is same as this except that yours has one more equation for calculating f0res

Ah yes I think they're the same. I just didn't see the 1/2pi term in your equation.

210Hz right?

Oops, yes!

Then, resonanceFreq = 163.8hz because f0res = 210hz, meff=410g and M=264g.

410g is the total mass of the piezo scanner, but it's not the same as effective mass. The data sheet doesn't provide meff directly but we can use Equation 10 from the link that I sent and calculate it from the resonance frequency (210Hz) and the stiffness (0.2N/um). I'm getting 115g for meff and 115.7Hz for the resonance frequency of the loaded piezo.

(Effective mass)
julia> (0.2*1e6)/((210*2pi)^2)
0.11487662544482741 (kg) --> 115g

(Resonance frequency with a load of 264g)
julia> 210 * sqrt(115/(115+264))
115.67742429171378 (Hz)

What happen if we need to change objective lens. Those lenses have similar weight?

All of our small Olympus objectives (10x, 20x, 40x UMPL-- series) have similar weight so we don't need to worry about that. If we switch to one of our heavier objectives then it will matter much more, so it would still be good to have some flexibility in the software.

Just let me know the final parameters you choose (window size, ratio, power threshold) and I can add a similar safety check to ImagineInterface.jl

[kdw503]

I think it could be even better to include some information about which time window failed the safety check. Especially since, as we mentioned, the problem may not be obvious when looking at the spectrum for the entire command (I assume that's what the new plot window is showing?)

Agreed.

If we take the square root of that power and divide that by the duration of the window then we get the average velocity within the window, so at least the velocity check can be performed within the same framework. As for how to set the velocity/power threshold, I think we can use something similar to what we used before? Or if that was too conservative we can raise the limit a bit.
Just let me know the final parameters you choose (window size, ratio, power threshold) and I can add a similar safety check to ImagineInterface.jl

Here is a fake code summarizing all the posts until now.
There are some parameters temporarily decided and undecided.
How about tuning the parameters and this code after my experiments in Imagine and your experiments in ImagineInterface?

# parameters
sample_rate = 100000
resonanceFreq = 115.7 # Hz
bandwidth = 50 # Hz
threshold = 0.01 # resonance frequency power ratio threshold 
mininum_power_th = 100 # not decided yet, need some experiment. I'll let you know later.
maxspeed = 2000 # piezo maxspeed in Imagine
velocity_th = maxspeed*2 # 200% increase from maxspeed
window_size = 8192 # samples
step_size = window_size/2 # samples

# safety check routine
frequency_resolution = sample_rate/window_size
window_size_in_time = window_size/sample_rate
for i in 1:step_size:total_samples
	fftCoeff = real_fft(piezo_control[i:i+window_size]) # I ignored BoundError for this fake code
	total_power = real(sum(fftCoeff.*conj(fftCoeff))) 
	DC_power = real(fftCoeff[1]*conj(fftCoeff[1]))
	AC_power = totol_power - DC_power
	if AC_power > minimum_power_th
		# check power around resonance frequency
		lower_bound = max(1, Int((resonanceFreq-bandwidth/2)/frequency_resolution))
		upper_bound = min(length(fftCoeff),Int((resonanceFreq+bandwidth/2)/frequency_resolution))
		resonance_power = real(sum(fftCoeff[lower_bound:upper_bound].*conj(fftCoeff[lower_bound:upper_bound])))
		ratio = resonance_power/AC_power
		if ratio > threshold
			error("anxial piezo control has excessive frequency components ($ratio > Threshold = $threshold %)
				around resonance frequency ($resonanceFreq Hz) in sample number from $i to $(i+window_size).
				Please check 'Frequency component spectrum'.")
		end
		# check velocity
		velocity = sqrt(AC_power)/window_size_in_time 
		if velocity > velocity_th
			error("anxial piezo control is too fast ($velocity > Threshold = $velocity_th)
				in between sample number from $i to $(i+window_size).
				Please check 'Waveform in time'.")
		end
	end
end

[Cody-G]

How about tuning the parameters and this code after my experiments in Imagine and your experiments in ImagineInterface?

Sounds good. But now I'm having second thoughts about this statement I made:

I think the total power (without DC) within the window should be the square of the total distance traveled by the piezo.

It's true that the total (DC-subtracted) power has units of squared distance, but I think it's not equal to (distance traveled)^2. Now I don't see a worthwhile shortcut for extracting velocity from the spectrum. But I think the velocity computation is cheap enough that maybe it doesn't matter:

distance_in_window = sum(abs.(diff(signal_in_window)))
velocity_in_window = distance_in_window / window_duration

[kdw503]

But I think the velocity computation is cheap enough that maybe it doesn't matter:

I agreed. And, we already have velocity check function also.

[kdw503]

I updated the algorithm and parameters.

# parameters
sample_rate = 100000
resonanceFreq = 115.7 # Hz
bandwidth = 20. # updated but not decided yet.
threshold = 0.02 # updated but not decided yet.
mininum_power_th = 100 # not decided yet.
maxspeed = 4000 # updated piezo maxspeed in Imagine
window_size = sample_rate # to make frequency_resolution be 1hz
step_size = window_size/2

### safety check routine
# 1. check velocity
error = velocity_check(piezo_control, maxspeed, strt, stop) # The function we already had
if error == true
    error("'axial piezo' control is too fast in between sample index $strt and $stop.")
end

# 2. check resonance frequency components
frequency_resolution = sample_rate/window_size
window_size_in_time = window_size/sample_rate
for i in 1:step_size:total_samples
	fftCoeff = real_fft(piezo_control[i:i+window_size]) # I ignored BoundError for this fake code
	total_power = real(sum(fftCoeff.*conj(fftCoeff))) 
	DC_power = real(fftCoeff[1]*conj(fftCoeff[1]))
	AC_power = totol_power - DC_power
	if AC_power > minimum_power_th
		# check power around resonance frequency
		lower_bound = max(1, Int((resonanceFreq-bandwidth/2)/frequency_resolution))
		upper_bound = min(length(fftCoeff),Int((resonanceFreq+bandwidth/2)/frequency_resolution))
		resonance_power = real(sum(fftCoeff[lower_bound:upper_bound].*conj(fftCoeff[lower_bound:upper_bound])))
		ratio = resonance_power/AC_power
		if ratio > threshold
			error("'axial piezo' control has excessive frequency components ($ratio > Threshold = $threshold %)
				around resonance frequency ($resonanceFreq Hz) in between sample index $i and $(i+window_size).
				Please check 'Frequency component spectrum'.")
		end
	end
end