%% Fan Switchoff RPM Profile
% Estimate the RPM profile of a fan blade as it slows down after switchoff.
%
% An industrial roof fan spinning at 20,000 rpm is turned off. Air resistance (with a negligible contribution from bearing friction) causes the fan rotor to stop in approximately 6 seconds. A high-speed camera measures the x-coordinate of one of the fan blades at a rate of 1 kHz.

fs = 1000;
t = 0:1/fs:6-1/fs;

rpm0 = 20000;

%
% Idealize the fan blade as a point mass circling the rotor center at a radius of 50 cm. The blade experiences a drag force proportional to speed, resulting in the following expression for the phase angle:
% \phi=2\pi f_0T(1-e^{-t/T}),
% where f_0 is the initial frequency and T=0\ldotp 75 second is the decay time.

a = 0.5;
f0 = rpm0/60;
T = 0.75;

phi = 2*pi*f0*T*(1-exp(-t/T));

%
% Compute and plot the x- and y-coordinates of the blade. Add white Gaussian noise of variance 0\ldotp 1^2.

x = a*cos(phi) + randn(size(phi))/10;
y = a*sin(phi) + randn(size(phi))/10;

plot(t,x,t,y)

%
% Use the rpmtrack function to determine the RPM profile. Type
% rpmtrack(x,fs)
% at the command line to open the interactive figure.
%
% Use the slider to adjust the frequency resolution of the time-frequency map to about 11 Hz. Assume that the signal component corresponds to order 1 and set the end time for ridge extraction to 3.0 seconds. Use the crosshair cursor in the time-frequency map and the Add button to add three points lying on the ridge. (Alternatively, double-click the cursor to add the points at the locations you choose.) Click Estimate to track and extract the RPM profile.
%
% Verify that the RPM profile decays exponentially. On the Export tab, click Export and select Generate MATLAB Script. The script appears in the Editor.

% MATLAB Code from rpmtrack GUI

% Generated by MATLAB 9.4 and Signal Processing Toolbox 8.0

% Generated on 18-Dec-2017 19:00:35

% Set sample rate
fs = 1000.0000;
% Set order of ridge of interest
order = 1.0000;
% Set ridge points on ridge of interest
ridgePoints = [...
    0.4501 179.8246;...
    0.9944 88.5965;...
    2.4161 11.4035];
% Estimate RPM
[rpmOut,tOut] = rpmtrack(x,fs,order,ridgePoints,...
    'Method','stft',...
    'FrequencyResolution',11.1612,...
    'PowerPenalty',Inf,...
    'FrequencyPenalty',0.0000,...
    'StartTime',0.0000,...
    'EndTime',3.0000);

% Run the script. Display the RPM profile in a semilogarithmic plot.

semilogy(tOut,rpmOut)
ylim([500 20000])

%
%
%
%
% Copyright 2012 The MathWorks, Inc.