Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Solution to dy/dan's Fan Slowing Problem
% Kenneth Finnegan, 2012
% kennethfinnegan.blogspot.com
%
% Predict end of spinning object
% Solution to http://blog.mrmeyer.com/?p=13245
clear all;
close all;
% Single step through video and pull blade hits on bottom of frame
rawdata = [ 0.36, 0.53, 0.70, 0.88, 1.06 1.23 1.42 ...
1.61 1.80 1.98 2.16 2.37 2.56 2.76 2.96 3.16 ...
3.37 3.58 3.79 4.00 4.21 4.43 4.65 4.88 5.10 ...
5.33 5.57 5.83 6.07 6.31 6.56 6.80 7.06 7.30 ...
7.57 7.83 8.09 8.35 8.62 8.90 9.17 9.46 9.74 ...
10.03 10.32 10.61 10.91 11.21 11.51 11.83 12.14 ...
12.47 12.81 13.13 13.46 13.81 14.14 14.51 14.83 ...
15.19 15.54 15.93 16.31 16.71 17.08 17.49 17.90 ...
18.31 18.74 19.16 19.60 20.05 20.48 20.95 21.42 ...
21.88 22.36 22.85 23.35 23.86 24.39 24.92 25.45];
% Convert blade hit times to rotational velocity / 4
% The /4 is because I counted all four blades hitting bottom
% of the frame, not just one, to get more data
for index = 1:(length(rawdata) - 1)
rotvel(index) = 1 / (rawdata(index + 1) - rawdata(index));
end
% Clean up the data a little bit due to frame rate aliasing and
% me being human
for index = 1:(length(rotvel) - 2)
rotvel(index) = mean(rotvel(index:index+2));
end
figure(1);
plot(rawdata(1:length(rawdata)-1), rotvel);
xlabel('Time [seconds]');
ylabel('Rotational Velocity [blades per second]');
% Take derivative to find acceleration
for index = 1:(length(rotvel) - 1)
rotaccel(index) = (rotvel(index + 1) - rotvel(index)) / ...
(rawdata(index + 1) - rawdata(index));
end
figure(2);
plot(rawdata(1:length(rawdata)-2), rotaccel);
xlabel('Time [seconds]');
ylabel('Acceleration [blades / sec / sec]');
figure(3);
hold on;
scatter(rotvel(1:(length(rotvel) - 1)), rotaccel);
xlabel('Rotational Velocity [blades / sec]');
ylabel('Rotational Acceleration [blades /sec / sec');
% Considering that the acceleration (a = F / m) is somehow
% related to the velocity of the fan, we'll say that it's
% probably a second order polynomial, because air drag is N^2,
% and we might be able to get away with calling friction N + C.
% (Friction is really probably a mess, and we're ignoring the
% motor's induction implications)
accelfromvel = polyfit(rotvel(1:(length(rotvel) - 1)), rotaccel, 2);
plot(1.5:0.1:6, polyval(accelfromvel, 1.5:0.1:6), 'r');
figure(4); hold on;
% Now that we have this acceleration(velocity) = N^2 + N + C model,
% we will take each velocity point we got from the original video
% and extrapolate them to the end
stepsize = 0.01;
for index = 1:(length(rotvel) - 1)
clear model_v model_t model_index
model_index = 1;
model_t(1) = rawdata(index);
model_v(1) = rotvel(index);
% Given those intial conditions, walk out to the future
while model_v(model_index) > 0
model_t(model_index + 1) = model_t(model_index) + stepsize;
model_v(model_index + 1) = model_v(model_index) + ...
polyval(accelfromvel, model_v(model_index)) * stepsize;
model_index = model_index + 1;
end
endtimes(index) = model_t(length(model_t));
plot(model_t, model_v, 'r');
end
plot(rawdata(1:length(rawdata)-1), rotvel, 'bd-');
axis([0 (max(endtimes) + 5) 0 6]);
xlabel('Time [seconds]');
ylabel('Rotational Velocity [blades per second]');
figure(5);
hist(endtimes, linspace(min(endtimes), max(endtimes), 10));
xlabel('Modeled End Time [seconds]');
ylabel('Number of simulations [#]');
PredictedFinished = mean(endtimes)
% Result for 0.01 stepsize with 3 point smoothing = 46.7 seconds
%
% Due to aliasing of the data, The model oscillates between ~46 sec
% for odd smoothing and ~62 seconds for even smoothing
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.