how to obtain the TSA standard deviation?

I am performing a TSA of a measurement of a rotating square bar at 500 rpm for 45 revolutions. the problem, as you can see in the picture, is that there is a noticeable deviation (calculated by iterations and visual inspection i found that the best average speed is 499.508 rpm). i would like to automatize this process by minimizing the std of the TSA but i am unable to find any option to obtain this information.

[Meg Noah]

I think what you might need is called 'Welford’s method'. Conceptualized by B. P. Welford and published widely by Knuth in The Art of Computer Programming (Volume 2, Seminumerical Algorithms, 3rd edn., p. 232.).

Take a look at: https://www.geeksforgeeks.org/web-tech/expression-for-mean-and-variance-in-a-running-stream/

where:

Answer 1

[William Rose]

@Antonio,

You can esitmate the mean frequency by counting zero crossings. Depending on the level of noise in the data, you may want to do a little bit of lowpass filtering before detecting zero crossings. If your goal is to find the right frequency for timke-synchronous averaging, I think this is easier than minimizing the stadard deviation of the TSA. Demonstrate with example:

Make a signal:

fs=1000;            % sampling rate (Hz)
f0=499.5/60;        % rotation frequency (Hz) (499.5 rpm)
Tmax=45/f0;         % time for 45 revolutions
t=0:1/fs:Tmax;      % time vector
y=cos(4*2*pi*f0*t)+0.05*randn(size(t));  % sinusoid with noise

The y=cos() includes a factor of 4 because the object has 4 corners, so it produces 4 peaks per revolution, as is seen in your original plot.

Use the measured data to estimate the frequency of revolution:

izc=find(diff(sign(y-mean(y)))==2);                 % find indices of upward zero crossings
f0est=(length(izc)-1)/(t(izc(end))-t(izc(1)))/4;    % esitmate the rotation rate (Hz)

The division by 4 above is because the bar produces 4 peaks per cycle as it rotates.

fprintf('Est. rotation rate=%.3f rpm.\n',f0est*60) 

Est. rotation rate=499.535 rpm.

This estimate of the rotation frequency is quite good. You will get a more accurate estimate by

a) sampling at a faster rate

b) interpolating to find the first and last zero crossing times more exactly

c) smoothing the noisy data (slightly) before finding zero crossings

Use this estimated frequency to calculate the time-synchronous average signal.

tCyc=mod(t,1/f0est)*f0est;  % time of each data point computed as fraction of a cycle
plot(tCyc,y,'rx'); 
grid on; xlabel('Time (fraction of Cycle)'); ylabel('Amplitude')

●

The timing of the samples within each cycle is not exactly the same, from one cycle to the next. Thererfore, if you wish to estimate standard deviation at various times throughout the cycle, you make want to interpolate to a common time-within-the-cycle.

[William Rose]

@Antonio, As I mentioned above, if you wish to estimate average or standard deviation or both at various times throughout the cycle, you will need to interpolate to a common time-within-the-cycle. Therefore I will make a signal and esitmate the rotation frequency, as above.

Make a signal. This signal is slightly more complicated than the one constructed previously: the four peaks per cycle have slightly unequal average heights. This varation in peak height is smaller than the st.dev. of the noise, so it will not be easy to detect. @Antonio's original figure hinted at this possibility.

fs=1000;            % sampling rate (Hz)
f0=499.5/60;        % rotation frequency (Hz) (499.5 rpm)
Tmax=45/f0;         % time for 45 revolutions
t=0:1/fs:Tmax;      % time vector
y=cos(4*2*pi*f0*t)+0.02*sin(2*pi*f0*t)+0.05*randn(size(t));  % sinusoid with noise

Estimate the frequency of revolution from the data. To improve the accuracy of the frequency estimate, smooth the data with a 5-point moving average first, and interpolate to estimate first and last zero crossing times. These are two of the methods for improving the estimate which I recommended in my answer above.

ys=smooth(y,5);                         % smooth with 5-point moving average
izc=find(diff(sign(ys-mean(ys)))==2);   % indices of upward zero crossings
tzc1=(t(izc(1))*ys(izc(1)+1)-t(izc(1)+1)*ys(izc(1)))/(ys(izc(1)+1)-ys(izc(1))); % tzc1 interp
tzce=(t(izc(end))*ys(izc(end)+1)-t(izc(end)+1)*ys(izc(end)))/(ys(izc(end)+1)-ys(izc(end))); % tzce interp
f0est=(length(izc)-1)/(tzce-tzc1)/4;    % estimate the rotation rate (Hz)
% Divide by 4 above because the bar produces 4 peaks per cycle as it rotates.
fprintf('Est. rotation rate = %.4f Hz = %.2f rpm.\n',f0est,f0est*60) 

Est. rotation rate = 8.3248 Hz = 499.49 rpm.

Use the estimated frequency to calculate the time within each cycle.

tCyc=mod(t,1/f0est)*f0est;  % time of each data point computed as fraction of a cycle

Plot the data, and plot zoomed-in views of the 2nd and 4th peaks.

figure
subplot(211), plot(tCyc,y,'rx'); 
grid on; xlabel('Time (fraction of cycle)'); ylabel('Amplitude')
subplot(223), plot(tCyc,y,'rx'); 
grid on; xlabel('Time (fraction of cycle)'); ylabel('Amplitude');
xlim([.2,.3]); ylim([.8,1.2])
subplot(224), plot(tCyc,y,'rx');
grid on; xlabel('Time (fraction of cycle)'); ylabel('Amplitude')
xlim([.7,.8]); ylim([.8,1.2])

The variation in peak heights is not obvious on the figure above. Since the data point timing is different on each cycle, resample the data using a sampling frequency that is an integer multiple of the estimated frequency. That will produce samples that are at the same time in each cycle. Then we will be able to compute the time-syncronous average.

nsc=240;                        % number of samples per cycle
fsre=f0est*nsc;                 % resampling frequency
tre=t(1):1/fsre:t(end);         % resamping time vector
yre=interp1(t,y,tre,'pchip');   % interpolate to the new sampling rate

Use the estimated frequency to calculate the time within each cycle.

tCycre=mod(tre,1/f0est)*f0est;  % time of each data point computed as fraction of a cycle

Reshape the signal into a matrix where each row is a complete cycle. Then compute the time-synchronous-average signal.

lastElem=nsc*floor(length(yre)/nsc);  % last element of last complete cycle
yreshape=reshape(yre(1:lastElem),nsc,[])';

Compute time-synchronous average

yTSA=mean(yreshape);
tTSA=(0:nsc-1)/nsc;  % time base for yTSA, for plotting

Plot the resampled and mean data

figure
subplot(211), plot(tCycre,yre,'rx',tTSA,yTSA,'-k.'); 
grid on; xlabel('Time (fraction of cycle)'); ylabel('Amplitude'); legend('Data','TSA')
% Plot zoomed-in views of the 2nd and 4th peaks.
subplot(223), plot(tCycre,yre,'rx',tTSA,yTSA,'-k.'); 
grid on; xlabel('Time (fraction of cycle)'); ylabel('Amplitude');
xlim([.2,.3]); ylim([.8,1.2]); legend('Data','TSA')
subplot(224), plot(tCycre,yre,'rx',tTSA,yTSA,'-k.');
grid on; xlabel('Time (fraction of cycle)'); ylabel('Amplitude')
xlim([.7,.8]); ylim([.8,1.2]); legend('Data','TSA')

The difference in the peaks heights is obvious in the time-synchronous average signal.