Solution to dy/dan's Fan Slowing Problem
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
% 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