{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Steps. Step by step.\n\nThis EDA is about steps. We will locate and visualize individual steps and construct the \"step rate\" and \"current step duration\" features. We will also find array indices for individual steps that can be used for further feature construction. \n\nUpdate:\nThere is a dataset of mine with discussed step-related features available. They are calculated for both TDCSFoG and DEFoG datasets with the same file structure as the training data.\nhttps://www.kaggle.com/datasets/vrbaryshev/step-data\n\n\n## Table of contents:\n1. [Introduction](#Introduction)\n2. [Discussing problems with AccV and AccAP and explaining why we'll want to use AccMl for step detection.](#Discussing)\n3. [Visualizing steps: explaining the problem and suggesting continious wavelet transformation as a solution.](#Visualizing)\n4. [Locating steps and constructing the \"step rate\" and \"current step duration\" features. We'll do it step by step.](#Locating)\n5. [Summarizing it all in one \"detect_steps\" function.](#Summarizing)\n6. [Testing the function and discussing some issues.](#Testing)\n7. [Conclusion](#Conclusion)","metadata":{}},{"cell_type":"markdown","source":"# 1. Introduction <a name=\"Introduction\"></a>\n\nSteps are the core part of any walking. I believe and that step rate itself should be an important feature for identifying FoG and other features calculated over individual steps can be important too. However, correctly splitting the process of walking into individual steps solely with acceleration data might be a complicated task. In this article I've tried to create, thoroughly explain and discuss a reliable solution for that task.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd \nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nfrom math import ceil\nimport scipy as sp\nfrom scipy import signal\nfrom scipy import fftpack\nimport pywt  # internet says pywt wavelets are better then scipy's. I'll believe that.\nimport tqdm\nimport re\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:02.996367Z","iopub.execute_input":"2023-05-25T15:13:02.998865Z","iopub.status.idle":"2023-05-25T15:13:03.555643Z","shell.execute_reply.started":"2023-05-25T15:13:02.998827Z","shell.execute_reply":"2023-05-25T15:13:03.554330Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#defines\ndatasets_folder = '/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/'\ndata_groups_names = ('tdcsfog', 'defog')\nfrequencies = (128, 100)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:03.557846Z","iopub.execute_input":"2023-05-25T15:13:03.558185Z","iopub.status.idle":"2023-05-25T15:13:03.562466Z","shell.execute_reply.started":"2023-05-25T15:13:03.558155Z","shell.execute_reply":"2023-05-25T15:13:03.561551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We'll read only tdcsfog data, as it'll be sufficient for the purpose of this EDA. I believe that the results are valid for other datasets too. Sorry for converting all the names to snake_case.","metadata":{}},{"cell_type":"code","source":"def read_data_file(test_id, subfolder_name = 'train/tdcsfog/'):\n    data = pd.read_csv(datasets_folder+subfolder_name + str(test_id) + '.csv')\n    data.columns = data.columns.str.replace(r'(?<!^)(?=[A-Z])', '_').str.lower()\n    data = data.rename(columns={\"acc_m_l\": \"acc_ml\", \"acc_a_p\": \"acc_ap\"})\n    return data","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:03.564061Z","iopub.execute_input":"2023-05-25T15:13:03.564366Z","iopub.status.idle":"2023-05-25T15:13:03.579390Z","shell.execute_reply.started":"2023-05-25T15:13:03.564339Z","shell.execute_reply":"2023-05-25T15:13:03.578256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_tdcsfog = pd.read_csv('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/tdcsfog_metadata.csv')\ndf_tdcsfog.columns = df_tdcsfog.columns.str.replace(r'(?<!^)(?=[A-Z])', '_').str.lower()\ndf_tdcsfog = df_tdcsfog.rename(columns={\"acc_m_l\": \"acc_ml\", \"acc_a_p\": \"acc_ap\"})\n\n# adding 'data' column to the whole table, so we can have access to all of the time series's from one table\ndf_tdcsfog['data'] = df_tdcsfog['id'].apply(read_data_file)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:03.582369Z","iopub.execute_input":"2023-05-25T15:13:03.582827Z","iopub.status.idle":"2023-05-25T15:13:24.213448Z","shell.execute_reply.started":"2023-05-25T15:13:03.582793Z","shell.execute_reply":"2023-05-25T15:13:24.212187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Discussing problems with AccV and AccAP and explaining why we'll want to use AccML for step detection. <a name=\"Discussing\"></a>\n* AccV means vertical acceleration sensor data\n* AccAP means forward-backward acceleration sensor data\n* AccML means left-right acceleration sensor data\n\nIntuitively, I began with looking at AccV, as it was supposed to measure vertical movements. I've never worked with biological data before, so I imagined a walking platform that steadily carries the sensors and keeps their orientation. The reality turned out to be different: people seem to lean forward a lot when they walk and they actively use that leaning techique to keep balance. As a result, the back-mounted sensor gets tilted (mostly) forward with constantly changing tilting angle. As a result, most of measured acceleration swings in AccV and AccAP correspond to the leaning (tilting) angle changes. Let us demonstrate that with some data:","metadata":{}},{"cell_type":"code","source":"subject = df_tdcsfog.loc[0,'data']   # we'll use the first patient in the database here for demonstration purposes\ntime = np.linspace(0, len(subject)/frequencies[0], len(subject))   \n\nplt.figure(figsize=(18,5))\nplt.plot(time, subject['acc_ap'])\nplt.title('AccAP time series')\nplt.gca().set_xlim(time[0], time[-1])\nplt.gca().set_ylabel('Acceleration, m/s^2')\nplt.gca().set_xlabel('Time, s')\nplt.show()\n\nprint('Average AccV is',subject['acc_v'].mean())\nprint('Average AccAP is',subject['acc_ap'].mean())","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:24.215095Z","iopub.execute_input":"2023-05-25T15:13:24.215414Z","iopub.status.idle":"2023-05-25T15:13:24.576215Z","shell.execute_reply.started":"2023-05-25T15:13:24.215385Z","shell.execute_reply":"2023-05-25T15:13:24.575189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As AccAP data tell us, sensor measures positive acceleration for almost the whole way. This would've meant constantly increasing speed if it was \"real\" forward acceleration. What we actually see here is projection of Earth's gravity force on the sensor's AP axis that gets there because the sensor gets tilted forward. Correspondigly, average AccV is less the g value because the V axis isn't actually vertical. As a result, AccAP and AccV mostly measure the tilting angle. We can also confirm this by calculating sum of their squares:","metadata":{}},{"cell_type":"code","source":"print('Average (AccV^2+AccAP^2)^0.5 - g =',((subject['acc_v']**2 + subject['acc_ap']**2)**0.5).mean() - 9.81)","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:24.577733Z","iopub.execute_input":"2023-05-25T15:13:24.578745Z","iopub.status.idle":"2023-05-25T15:13:24.587185Z","shell.execute_reply.started":"2023-05-25T15:13:24.578705Z","shell.execute_reply":"2023-05-25T15:13:24.586005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, mean sum of squares of $AccV{^2} + AccAP{^2}$ is very close to $g{^2}$. This also means that the AccML sensor is almost free of gravity force projection. So we'll want use AccML for detecting steps.","metadata":{}},{"cell_type":"markdown","source":"# 3. Visualizing steps: explaining the problem and suggesting continious wavelet transformation as a solution. <a name=\"Visualizing\"></a>\n\nIntuitively, steps seemed to be easy to locate. We have measurements of gait here and with steps being the core mechanic, they surely have to stand out in any measurements. Let us demonstrate some complications with a 5 second time window of AccML signal:","metadata":{}},{"cell_type":"code","source":"t1, t2 = 20*frequencies[0], 25*frequencies[0]\nplt.figure(figsize=(18,2))\nplt.plot(time, subject['acc_ml'])\nplt.title('AccML for 5 seconds')\nplt.gca().set_xlim(time[t1], time[t2])\nplt.gca().set_ylabel('Acceleration, m/s^2')\nplt.gca().set_xlabel('Time, s')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:24.588631Z","iopub.execute_input":"2023-05-25T15:13:24.589058Z","iopub.status.idle":"2023-05-25T15:13:24.831420Z","shell.execute_reply.started":"2023-05-25T15:13:24.589007Z","shell.execute_reply":"2023-05-25T15:13:24.830318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Intuitively, I expected AccML to cross zero back and forth on every step and thought about locating thouse crossings to measure steps. This random 5 second walking interval shows a different picture. Signal has all kind of noises and random oscillations and steps don't visually stand out at all. We are supposed to see about 10 steps in 5 seconds, by the way. But with that kind of signal we'll need some kind of spectral approach to detect them. Let us briefly look at what Fourier transformation results in for the whole session:","metadata":{}},{"cell_type":"code","source":"temp_fft = fftpack.fft(subject['acc_ml'].values)    # this gives the Fourier transformation\nfftfreq = fftpack.fftfreq(len(temp_fft), 1. / frequencies[0])  # this makes frequency values for points in temp_fft\ni = fftfreq > 0\n\nplt.figure(figsize=(18,2))\nplt.plot(fftfreq[i], abs(temp_fft[i]))\nplt.title('AccML spectrum')\nplt.gca().set_xlim(0, 10)\nplt.gca().set_ylabel('Signal density')\nplt.gca().set_xlabel('Frequency, Hz')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:24.833140Z","iopub.execute_input":"2023-05-25T15:13:24.834118Z","iopub.status.idle":"2023-05-25T15:13:25.082998Z","shell.execute_reply.started":"2023-05-25T15:13:24.834085Z","shell.execute_reply":"2023-05-25T15:13:25.081760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As the graph shows, the spectrum isn't smooth and nice at all and doesn't give any outstanding step frequency. We should expect this as the whole time graph has about 35 seconds in it with and includes standing up, acceleration, walking, deceleration and possible FoG events. We'll want to look at smaller time windows to analyze walking as our goal is to locate and describe every step. So we'll use continious wavelet transformation (CWT) to make steps vizible:","metadata":{}},{"cell_type":"code","source":"wavelet='morl'       # we'll use Morlet wavelet\nscales = np.exp(np.arange(3, 6, 0.05))   # wavelet scales define the frequency values we get\n# continious wavelet transformation from pywt library\n# coeff gives how much of each frequency we have at each time moment\n# freq is the corresponding list of frequency values\ncoeff, freq = pywt.cwt(subject['acc_ml'], scales, wavelet)\n\nplt.figure(figsize=(18,2))\nplt.title('CWT of AccML')\nplt.pcolormesh(time, freq*frequencies[0], abs(coeff), cmap='jet')\nplt.gca().set_yscale('log')\nplt.gca().get_yaxis().set_major_formatter(mpl.ticker.ScalarFormatter())\nplt.yticks(ticks=[0.5, 1, 2, 3, 4])\nplt.gca().set_ylabel('Frequency, Hz')\nplt.gca().set_xlabel('Time, s')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:25.084640Z","iopub.execute_input":"2023-05-25T15:13:25.085284Z","iopub.status.idle":"2023-05-25T15:13:25.893540Z","shell.execute_reply.started":"2023-05-25T15:13:25.085246Z","shell.execute_reply":"2023-05-25T15:13:25.892683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And here we see the steps! Well, mostly) Hopefully) They are presented as a chain of green spots around frequency of 1 Hz. It is known that humans walk at about 2 Hz and we see exactly that, because the patient makes and the graph describes two steps per AccML cycle, one with each leg. Now let us locate and count them.","metadata":{}},{"cell_type":"markdown","source":"# 4. Locating steps and constructing the \"step rate\" and \"current step duration\" features. <a name=\"Locating\"></a>\n#### *We'll do it step by step )*\n\nSo we want to extract the green dots from the above graph and present them as features. Let us start with locating the maximum absolute value of CWT at each time slice:","metadata":{}},{"cell_type":"code","source":"coeff_argmax_index = np.argmax(abs(coeff), 0)\n\nplt.figure(figsize=(18,3))\nplt.plot(time, freq[coeff_argmax_index]*128)\nplt.title('Noisy step rate')\nplt.gca().set_xlim(time[0], time[-1])\nplt.gca().set_ylim(0, 4)\nplt.gca().set_ylabel('Frequency, Hz')\nplt.gca().set_xlabel('Time, s')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:25.897361Z","iopub.execute_input":"2023-05-25T15:13:25.897923Z","iopub.status.idle":"2023-05-25T15:13:26.212439Z","shell.execute_reply.started":"2023-05-25T15:13:25.897890Z","shell.execute_reply":"2023-05-25T15:13:26.210959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is our step rate. It is actually quite real, but a bit noisy at this time. It also jumps and drops to wrong frequency values a lot. We'll have to clear it. Now that is a bit tricky point. Slices of CWT image along any horizontal line give us our signal at that frequency. This is  somewhat like a frequency filter. So we'll make a slice of CWT along the line composed by our maximum image value array. We'll believe that it represents the step-related part of our original AccML data:","metadata":{}},{"cell_type":"code","source":"coeff_max = np.array([coeff[coeff_argmax_index[i], i] for i in range(coeff.shape[1])])\n\nfig,ax=plt.subplots(2,1,figsize=(18,7))\nax[0].set_title('Original AccML data')\nax[0].plot(time, subject['acc_ml'])\nax[0].set_xlim(time[0], time[-1])\nax[1].set_title('Processed AccML data')\nax[1].plot(time, coeff_max)\nax[1].set_xlim(time[0], time[-1])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:26.214014Z","iopub.execute_input":"2023-05-25T15:13:26.214372Z","iopub.status.idle":"2023-05-25T15:13:26.704989Z","shell.execute_reply.started":"2023-05-25T15:13:26.214342Z","shell.execute_reply":"2023-05-25T15:13:26.704165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This processed data visibly corresponds to most of the steps and actually crosses zero all the time. We can locate those crossings and count steps, but first we want to correct errors. Partly, the errors are pretty obvious: some steps are way too wide, while the ones next to them are super short.  Also, there are sharp turns of the signal line. We believe that the graph has all the necessary frequency information, but is noisy. We'll clear it again with the same CWT:","metadata":{}},{"cell_type":"code","source":"coeff_2, freq = pywt.cwt(coeff_max, scales, wavelet)\n\nfig,ax=plt.subplots(2,1,figsize=(18,5))\n\nax[0].set_title('CWT of original AccML data')\nax[0].pcolormesh(time, freq*frequencies[0], abs(coeff), cmap='jet')\nax[1].set_title('CWT of processed AccML data')\nax[1].pcolormesh(time, freq*frequencies[0], abs(coeff_2), cmap='jet')\n\nfor axx in ax:\n    plt.sca(axx)\n    axx.set_xlim(time[0], time[-1])\n    axx.set_ylabel('Frequency, Hz')\n    axx.set_yscale('log')\n    axx.get_yaxis().set_major_formatter(mpl.ticker.ScalarFormatter())\n    plt.yticks(ticks=[0.5, 1, 2, 3, 4])\n    axx.set_xlabel('Time, s')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:26.706470Z","iopub.execute_input":"2023-05-25T15:13:26.707100Z","iopub.status.idle":"2023-05-25T15:13:28.356572Z","shell.execute_reply.started":"2023-05-25T15:13:26.707044Z","shell.execute_reply":"2023-05-25T15:13:28.355621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we have a much clearer view. Steps are there and a lot of noise is gone. Let us collect the data along the maximum CWT image line:","metadata":{}},{"cell_type":"code","source":"coeff_argmax_index_2 = np.argmax(abs(coeff_2), 0)\ncoeff_max_2 = np.array([coeff_2[coeff_argmax_index_2[i], i] for i in range(coeff_2.shape[1])])\n\nfig,ax=plt.subplots(2,1,figsize=(18,5))\nax[0].set_title('Noisy step rate')\nax[0].plot(time, freq[coeff_argmax_index_2]*128)\nax[0].set_xlim(time[0], time[-1])\nax[0].set_ylim(0, 4)\nax[0].set_ylabel('Frequency, Hz')\nax[1].set_title('Processed AccML data')\nax[1].plot(time, coeff_max_2)\nfor axx in ax:\n    axx.set_xlim(time[0], time[-1])\n    axx.set_xlabel('Time, s')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:28.357731Z","iopub.execute_input":"2023-05-25T15:13:28.358863Z","iopub.status.idle":"2023-05-25T15:13:28.902294Z","shell.execute_reply.started":"2023-05-25T15:13:28.358827Z","shell.execute_reply":"2023-05-25T15:13:28.901120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The step rate became less noisy and processed ML data contains less errors. We have a problem in those graphs: np.argmax of CWT data pushes us back and forth between frequency bands. We want to stay in one band that represents steps and we believe that it is the most represented band. So we'll let the majority vote for the band by going for median:","metadata":{}},{"cell_type":"code","source":"coeff_argmax_index_21 = np.round(pd.Series(coeff_argmax_index_2).rolling(128, center=True, min_periods=1).median()).astype(int)\ncoeff_max_21 = np.array([coeff_2[coeff_argmax_index_21[i], i] for i in range(coeff_2.shape[1])])\n\nfig,ax=plt.subplots(2,1,figsize=(18,5))\nax[0].set_title('Noisy step rate')\nax[0].plot(time, freq[coeff_argmax_index_21]*128)\nax[0].set_xlim(time[0], time[-1])\nax[0].set_ylim(0, 4)\nax[0].set_ylabel('Frequency, Hz')\nax[1].set_title('Processed AccML data')\nax[1].plot(time, coeff_max_21)\nfor axx in ax:\n    axx.set_xlim(time[0], time[-1])\n    axx.set_xlabel('Time, s')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:28.903696Z","iopub.execute_input":"2023-05-25T15:13:28.904073Z","iopub.status.idle":"2023-05-25T15:13:29.491648Z","shell.execute_reply.started":"2023-05-25T15:13:28.904042Z","shell.execute_reply":"2023-05-25T15:13:29.490405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is much smoother. We can repeat the procedure for more smoothness:","metadata":{}},{"cell_type":"code","source":"coeff_3, freq = pywt.cwt(coeff_max_21, scales, wavelet)\ncoeff_argmax_index_3 = np.argmax(abs(coeff_3), 0)\ncoeff_argmax_index_31 = np.round(pd.Series(coeff_argmax_index_3).rolling(128, center=True, min_periods=1).median()).astype(int)\ncoeff_max_31 = np.array([coeff_2[coeff_argmax_index_31[i], i] for i in range(coeff_2.shape[1])])\n\nfig,ax=plt.subplots(2,1,figsize=(18,5))\nax[0].set_title('Noisy step rate')\nax[0].plot(time, freq[coeff_argmax_index_31]*128)\nax[0].set_xlim(time[0], time[-1])\nax[0].set_ylim(0, 4)\nax[0].set_ylabel('Frequency, Hz')\nax[1].set_title('Processed AccML data')\nax[1].plot(time, coeff_max_31)\nfor axx in ax:\n    axx.set_xlim(time[0], time[-1])\n    axx.set_xlabel('Time, s')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:29.493085Z","iopub.execute_input":"2023-05-25T15:13:29.493448Z","iopub.status.idle":"2023-05-25T15:13:30.184534Z","shell.execute_reply.started":"2023-05-25T15:13:29.493417Z","shell.execute_reply":"2023-05-25T15:13:30.183628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"But, inevitably, we'll need to somehow decide between two (equally) represented frequency bands. Also, as we are picking the brightest dots along each time slice, we are doomed to pick wrong data between steps when our step band is empty. So we'll use find_peaks() from numpy to detect bigger, wider and distant from each other peaks on our processed AccML data line and make another data array that is located along this peak line. We interpolate data in index space here:","metadata":{}},{"cell_type":"code","source":"max_cwt_points = sp.signal.find_peaks(abs(coeff_max_31), distance=20, width=20)[0]\nmax_cwt_points = np.concatenate(([0], max_cwt_points, [len(coeff_max_31)-1]))\n# and interpolating our line in wavelet space between the bigger peaks\nmax_cwt_line_indexes = np.round(np.interp( range(0, len(coeff_max_31)), max_cwt_points, \n                                                    coeff_argmax_index_31[max_cwt_points])).astype(int)\ncoeff_line = np.array([coeff_3[max_cwt_line_indexes[i], i] for i in range(coeff.shape[1])])\n\nplt.figure(figsize=(18,2))\nplt.plot(time, coeff_line)\nplt.title('AccML data along the CWT wider peak line')\nplt.gca().set_xlim(time[0], time[-1])\nplt.gca().set_xlabel('Time, s')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:30.185913Z","iopub.execute_input":"2023-05-25T15:13:30.186532Z","iopub.status.idle":"2023-05-25T15:13:30.451658Z","shell.execute_reply.started":"2023-05-25T15:13:30.186497Z","shell.execute_reply":"2023-05-25T15:13:30.450107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, since we used a smooth interpolated line in index space, our processed AccML data is reasonably smooth. Let us check where it crosses zero and locate the steps:","metadata":{}},{"cell_type":"code","source":"zero_crossings = np.where( np.diff(np.sign(pd.Series(coeff_line).rolling(10, center=True, min_periods=1).mean())))[0]\nzero_crossings = np.concatenate( ([0], zero_crossings, [len(coeff_line)]))    \n#filling each step with its duration while cutting possible outbursts with median\nstep_lengths = []\nfor i in range(1, len(zero_crossings)):\n    step_lengths = np.concatenate( (step_lengths, [zero_crossings[i] - zero_crossings[i-1]]*(zero_crossings[i] - zero_crossings[i-1]) ))\nstep_durations = pd.Series(step_lengths).rolling(32, center=True, min_periods=1).median()\n# making a nice and smooth step_rate array\nstep_rate = pd.Series( 1./step_durations )*frequencies[0]\nstep_rate = step_rate.where(step_rate<5, 0).rolling(frequencies[0], center=True, min_periods=1).mean()\n\nfig,ax=plt.subplots(3,1,figsize=(18,8))\nax[0].set_title('CWT of original AccML data')\nax[0].pcolormesh(time, freq*frequencies[0], abs(coeff), cmap='jet')\nax[0].set_ylabel('Frequency, Hz')\nax[0].set_yscale('log')\nax[0].get_yaxis().set_major_formatter(mpl.ticker.ScalarFormatter())\nplt.sca(ax[0])\nplt.yticks(ticks=[0.5, 1, 2, 3, 4])\nax[1].set_title('Finally (hopefully), clean step rate )')\nax[1].plot(time, step_rate)\nax[1].set_xlim(time[0], time[-1])\nax[1].set_ylim(0, 4)\nax[1].set_ylabel('Frequency, Hz')\nax[2].set_title('Duration of each step')\nax[2].plot(time, step_durations)\nfor axx in ax:\n    axx.set_xlim(time[0], time[-1])\n    axx.set_xlabel('Time, s')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:30.453146Z","iopub.execute_input":"2023-05-25T15:13:30.453499Z","iopub.status.idle":"2023-05-25T15:13:31.643903Z","shell.execute_reply.started":"2023-05-25T15:13:30.453467Z","shell.execute_reply":"2023-05-25T15:13:31.642687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Well, almost. We have steps. Our step rate is smooth and generally corresponds to the image. There are 3 visible problems on that graph:\n* wrong short step at the 14th second\n* wrong long step at the 23rd second\n* extra steps in the beginning\n\nAppearace of extra steps when there is no actual walking is to be solved by some sort of an \"if_walking\" filter. Currently, this is out of the scope of this EDA. But it might be added later. The first two problems are to be discussed after we collect the whole procedure in one function that we can use.","metadata":{}},{"cell_type":"markdown","source":"# 5. Summarizing it all in one \"detect_steps\" function. <a name=\"Summarizing\"></a>\n\n#### Parameters:\n* ***signal*** is a time series to detect steps from. I recommend the AccML data.\n* ***frequency*** is the signal frequency in Hz, which is 128 for our data, but different for the other dataset\n* ***smoothing*** is the number of wavelet smoothing cycles. I've seen artifacts that appear on 1 and disappear on 2.\n* ***faster_please*** is an option to skip the find_peaks stage. It creates artifacts, so I don't recommend it.\n* ***min_scale_log*** and ***max_scale_log*** determine the frequency band for all the CWT's. Wider is for handling far-off steprates, narrower is to avoid band miscalculations like on the graph above.\n\n#### Return value:\n* ***step_rate*** pd.Series with a smooth step rate value for every time point in dataset\n* ***step_durations*** pd.Series with flat plates as high as step duration in seconds as a time series for every time point in dataset\n* ***zero_crossings*** np.array with indices of time point that mark beginnings of new steps\n\nAll parameters except the signal have my suggested default values. Steps are considered to begin and end with two feet on the ground and the body weighting them equally. Some might consider those moments to be middles of steps and recalculate steps in a different way.\n","metadata":{}},{"cell_type":"code","source":"def detect_steps(signal, frequency=128, smoothing=2, faster_please=False, min_scale_log=3.6, max_scale_log=5.4):\n    \n    # scales for wavelets that define frequencies, exp for more even distribution\n    # choosing border values is tricky\n    # too narrow and we lose some slow or fast steprates\n    # too wide and we get the steprate wrong for unusual walking wavelet spectrum cases\n    # check patient #6 in the tdcsfog dataset as an example of an unusual pattern\n    scales = np.exp(np.arange(min_scale_log, max_scale_log, 0.05))\n    wavelet='morl'  # chosing the Morlet wavelet\n    \n    # transforming the signal, preferrably the 'AccMl' one\n    coeff, freq = pywt.cwt(signal, scales, wavelet)\n    # finding the brightest dots on every time slice\n    coeff_argmax_index = np.argmax(abs(coeff), 0)\n    # and collecting their coeff values in a new \"signal\", which already looks like our step rate, but is a bit noisy\n    coeff_max = np.array([coeff[coeff_argmax_index[i], i] for i in range(coeff.shape[1])])\n    \n    # repeating the smoothing procedure\n    for i in range(smoothing):\n        # transforming the noisy step rate again to clear it\n        coeff, freq = pywt.cwt(coeff_max, scales, wavelet)\n        # finding the brightest dots on every time slice\n        coeff_argmax_index = np.argmax(abs(coeff), 0).astype(int)\n        # and smoothing their indices with median to cut some outbursts off\n        coeff_argmax_index = np.round(pd.Series(coeff_argmax_index).rolling(128, center=True, min_periods=1).median()).astype(int)\n        # collecting the values along our line of indices\n        coeff_max = np.array([coeff[coeff_argmax_index[i], i] for i in range(coeff.shape[1])])\n    \n    # this smoothing round is optional\n    if not faster_please:\n        # finding hopefully bigger peaks on our line\n        # hopefully skiping too narrow or too close ones that are likely to be outbursts\n        max_cwt_points = sp.signal.find_peaks(abs(coeff_max), distance=20, width=20)[0]\n        max_cwt_points = np.concatenate(([0], max_cwt_points, [len(signal)-1]))\n        # and interpolating our line in wavelet space between the bigger peaks\n        max_cwt_line_indexes = np.round(np.interp( range(0, len(signal)), max_cwt_points, \n                                                    coeff_argmax_index[max_cwt_points])).astype(int)\n        coeff_max = np.array([coeff[max_cwt_line_indexes[i], i] for i in range(coeff.shape[1])])\n    \n    # finding zeroes on this smooth line to separate individual steps from other steps    \n    zero_crossings = np.where( np.diff(np.sign(pd.Series(coeff_max).rolling(10, center=True, min_periods=1).mean())))[0]\n    zero_crossings = np.concatenate( ([0], zero_crossings, [len(signal)]))    \n    #filling each step with its duration while cutting possible outbursts with median\n    step_lengths = []\n    for i in range(1, len(zero_crossings)):\n        step_lengths = np.concatenate( (step_lengths, [zero_crossings[i] - zero_crossings[i-1]]*(zero_crossings[i] - zero_crossings[i-1]) ))\n    step_durations = pd.Series(step_lengths).rolling(32, center=True, min_periods=1).median()\n    # making a nice and smooth step_rate array\n    step_rate = pd.Series( 1./step_durations )*frequency\n    step_rate = step_rate.where(step_rate<5, 0).rolling(frequency, center=True, min_periods=1).mean()\n    \n    return step_rate, step_durations, zero_crossings","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:31.645840Z","iopub.execute_input":"2023-05-25T15:13:31.646645Z","iopub.status.idle":"2023-05-25T15:13:31.667592Z","shell.execute_reply.started":"2023-05-25T15:13:31.646600Z","shell.execute_reply":"2023-05-25T15:13:31.666219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This function is autonomous and \"copy&paste\" friendly. It is for everyone to use, but please upvote this work if you use it and mention me in credits.","metadata":{}},{"cell_type":"markdown","source":"# 6. Testing the function and discussing some issues. <a name=\"Testing\"></a>\n\nLet us make a tool for displaying test results:","metadata":{}},{"cell_type":"code","source":"def feature_display(subject, faster_please=False, smoothing=2, min_scale_log=3.6, max_scale_log=5.4 ):\n\n    signal = subject['acc_ml']  \n    time = np.linspace(0, len(signal)/frequencies[0], len(signal)) \n    scales = np.exp(np.arange(min_scale_log, max_scale_log, 0.05))\n    wavelet='morl'\n    coeff, freq = pywt.cwt(signal, scales, wavelet)\n    print('Min frequency =', np.min(freq)*frequencies[0])\n    print('Max frequency =', np.max(freq)*frequencies[0])\n    \n    step_rate, step_durations, step_borders = detect_steps(subject['acc_ml'], \n                                                               frequency=frequencies[0], \n                                                               faster_please=faster_please,\n                                                              min_scale_log=min_scale_log,\n                                                              max_scale_log=max_scale_log)\n    \n    fig,ax=plt.subplots(3,1,figsize=(18,7))\n    \n    ax[0].set_title('AccML CWT data')\n    ax[0].pcolormesh(time, freq*frequencies[0], abs(coeff), cmap='jet')\n    ax[0].set_ylabel('Frequency, Hz')\n    ax[0].set_yscale('log')\n    ax[0].get_yaxis().set_major_formatter(mpl.ticker.ScalarFormatter())\n    plt.sca(ax[0])\n    plt.yticks(ticks=[0.5, 1, 2, 3, 4])\n    \n    ax[1].set_title('Step rate')\n    ax[1].plot(time, step_rate)\n    ax[1].set_xlim(time[0], time[-1])\n    ax[1].set_ylim(0, 5)\n    ax[1].set_ylabel('Frequency, Hz')\n    \n    ax[2].set_title('Step durations')\n    ax[2].plot(time, step_durations/frequencies[0])\n    ax[2].set_xlim(time[0], time[-1])\n    ax[2].set_ylim(0, 1.5)\n    ax[2].set_ylabel('Time, s')\n    \n    plt.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:31.669090Z","iopub.execute_input":"2023-05-25T15:13:31.669522Z","iopub.status.idle":"2023-05-25T15:13:31.691920Z","shell.execute_reply.started":"2023-05-25T15:13:31.669491Z","shell.execute_reply":"2023-05-25T15:13:31.690955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can test the function. Let us begin with our initial case, the first patient in our data:","metadata":{}},{"cell_type":"code","source":"feature_display(df_tdcsfog.loc[0,'data'])","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:31.693466Z","iopub.execute_input":"2023-05-25T15:13:31.694843Z","iopub.status.idle":"2023-05-25T15:13:33.190499Z","shell.execute_reply.started":"2023-05-25T15:13:31.694807Z","shell.execute_reply":"2023-05-25T15:13:33.189262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, this step rate is correct and both miscalculated steps that were present in the \"step by step\" chapter are now gone. This is because the frequency band we used here is narrower. Choice of that band is delicate. Wide band leads to mistakes and narrow one screws at high or low step rate series. Let us test some more series:","metadata":{}},{"cell_type":"code","source":"feature_display(df_tdcsfog.loc[100,'data'])","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:33.192020Z","iopub.execute_input":"2023-05-25T15:13:33.192403Z","iopub.status.idle":"2023-05-25T15:13:34.778874Z","shell.execute_reply.started":"2023-05-25T15:13:33.192370Z","shell.execute_reply":"2023-05-25T15:13:34.777624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_display(df_tdcsfog.loc[11,'data'])","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:34.780270Z","iopub.execute_input":"2023-05-25T15:13:34.780657Z","iopub.status.idle":"2023-05-25T15:13:37.069686Z","shell.execute_reply.started":"2023-05-25T15:13:34.780620Z","shell.execute_reply":"2023-05-25T15:13:37.068653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_display(df_tdcsfog.loc[14,'data'])","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:37.071014Z","iopub.execute_input":"2023-05-25T15:13:37.071332Z","iopub.status.idle":"2023-05-25T15:13:38.918445Z","shell.execute_reply.started":"2023-05-25T15:13:37.071305Z","shell.execute_reply":"2023-05-25T15:13:38.917375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_display(df_tdcsfog.loc[200,'data'])","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:38.920204Z","iopub.execute_input":"2023-05-25T15:13:38.920539Z","iopub.status.idle":"2023-05-25T15:13:40.537044Z","shell.execute_reply.started":"2023-05-25T15:13:38.920510Z","shell.execute_reply":"2023-05-25T15:13:40.535618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_display(df_tdcsfog.loc[307,'data'])","metadata":{"execution":{"iopub.status.busy":"2023-05-25T15:13:40.538694Z","iopub.execute_input":"2023-05-25T15:13:40.539033Z","iopub.status.idle":"2023-05-25T15:13:41.945886Z","shell.execute_reply.started":"2023-05-25T15:13:40.539003Z","shell.execute_reply":"2023-05-25T15:13:41.944666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All of those seem to work correctly if we don't look at the parts where there is no walking. Some points with step_rate>3 are suspicios. I honestly tried to walk at step_rate=4. It's possible, but hard. But we have time series with that kind of frequencies being clearly visible, we have to identify them as step rate.","metadata":{}},{"cell_type":"markdown","source":"# 7. Conclusion <a name=\"Conclusion\"></a>\n\n* A resonably reliable procedure for locating individual steps is developed\n* \"Step rate\" and \"step durations\" features are constructed\n* A tool for visualizing steps is presented\n\n#### Additionally\n* created CWT images might help constructing other features\n* some parts of those images actually correlate with FoG marks","metadata":{}}]}