{"cells":[{"metadata":{"_uuid":"a2fb38d1d89a3a958c51b979cba3b6e7fd3f03cb"},"cell_type":"markdown","source":"# Exploratory Data Analysis of the acustic signal on the frequency domain\n\nIn this notebook we will explore the shape and behavior of the frequency spectrum of the acoustic signals (training data only).\nWe will explore the following:\n* The average power spectrum of the training data\n* The average power spectrum as a function of the time to failure.\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n%matplotlib inline\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"sampling_frequency = 4e6\ntime_step = 1./sampling_frequency\n\nchunksize = 150_000\n\ntrain = pd.read_csv('../input/train.csv', iterator=True, \n                    chunksize=chunksize, \n                    dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e77727c284a0ba06cb31260360ce912771944779"},"cell_type":"markdown","source":"First, we compute the frequency power spectrum for every chunk. Since the spectrum for each chunk is too big, we will downscale the data to reduce the memory requirements by:\n* Smoothing the data using a 10-points moving window.\n* Keeping one every 10 samples\n* Keeping only frequencies < 450khz (larger frequencies have a negligible contribution).\n\nKeeping low memory requirements facilitates the display of the data using interactive plots."},{"metadata":{"trusted":true,"_uuid":"8ca04f23ca93343291046df38bc1b7cae0f68155"},"cell_type":"code","source":"from scipy import signal\n\ndef get_spectrum(input_signal):\n    \"\"\"\n    Get a pandas Series with the fourier power spectum for a given signal segment.\n    \"\"\"\n    input_signal = np.asarray(input_signal.values, dtype='float64')\n    \n    # Remove the mean  \n    input_signal -= input_signal.mean()  \n    \n    # Estimate power spectral density using a periodogram.\n    frequencies , power_spectrum = signal.periodogram(input_signal, sampling_frequency, scaling='spectrum')    \n    \n    # Run a running windows average of 10-points to smooth the signal.\n    power_spectrum = pd.Series(power_spectrum, index=frequencies).rolling(window=10).mean()        \n    \n    return pd.Series(power_spectrum)\n\ndef get_segment_spectrum(segment_df):\n    \"\"\"\n    Get the fourier power spectrum of a given segment.\n    \n    Returns the quake_time, frequencies, and power_spectrum\n    \"\"\"\n    \n    quake_time =segment_df['time_to_failure'].values[-1]\n    \n    _power_spectrum = get_spectrum(segment_df['acoustic_data']).dropna() \n    \n    # Keep only frequencies < 450khz (larger frequencies have a negligible contribution).\n    _power_spectrum = _power_spectrum[_power_spectrum.index<450_000]\n    \n    # Keep one every 10 samples\n    power_spectrum=_power_spectrum.values[::10]\n    frequencies=_power_spectrum.index.values[::10]    \n    \n    return quake_time, frequencies, power_spectrum","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1527d267293715293bb383429ed3c16db5660398"},"cell_type":"markdown","source":"Let's compute the frequency power spectrum for every segment."},{"metadata":{"trusted":true,"_uuid":"40c586ea224f7a279d2f337904c8e01df495f0bc"},"cell_type":"code","source":"quake_times = list()       \npower_spectrums = list()\n\nfor df in train:    \n    quake_time, _frequencies, power_spectrum = get_segment_spectrum(df)    \n    if df.shape[0]<chunksize:\n        continue\n    \n    frequencies=_frequencies\n    quake_times.append(quake_time)    \n    power_spectrums.append(power_spectrum)\n\npower_spectrums = np.stack(power_spectrums, axis=0)\nquake_times = pd.Series(quake_times)    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f6cd603a1393c3d6998d2bdbf403f5e2f3b182e8"},"cell_type":"code","source":"print(\"power_spectrums.shape:\",power_spectrums.shape)\nprint(\"power_spectrums size:\",power_spectrums.nbytes/(1024**2.),\"[mb]\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"543647778ca8075d6d15cd18bd05e642c4cf2f58"},"cell_type":"markdown","source":"## Explore the average power spectrum\n\nFirst, let's see how the mean power spectrum looks like. We will display the results using bokeh so they can be explored interactively."},{"metadata":{"trusted":true,"_uuid":"2f65283426c7a33ae3d7429ec20c826ea2fd0722"},"cell_type":"code","source":"from bokeh.layouts import column, row\nfrom bokeh.plotting import Figure, show\nfrom bokeh.io import output_notebook\nfrom bokeh.models import PrintfTickFormatter\nfrom bokeh.models import LinearAxis, Range1d\n\noutput_notebook() # Display Bokeh plots inline in a classic Jupyter notebooks\n\n# Compute the average spectrum\naverage_power_spectrum = pd.Series(power_spectrums.mean(axis=0))\naverage_power_spectrum = average_power_spectrum.rolling(window=10).mean()\n\n# \"@foo{(0.00 a)}\" # formats 1230974 as: 1.23 m\nTOOLTIPS = [ (\"x\", \"@x{(0.00 a)}Hz\"), (\"y\", \"$y\")   ]\n\npl = Figure(plot_width=800, plot_height=400,title=\"Average power spectrum\",tooltips=TOOLTIPS)\npl.line(frequencies, average_power_spectrum, line_color=\"navy\")\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency [hz]\"\npl.yaxis.axis_label = \"Power [V**2]\"\nshow(pl)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b49844dea58f3c3f51593710f7c782b1c7250240"},"cell_type":"markdown","source":"## Normalized average spectrum\nThe latter plot provides only a general picture of the signal spectrum. Signals with a very low total power (power spectrum integral) are average together with spectrums with a very large value. In consequence, this may result in certain features or peaks present in the spectrum being masked out.\n\nLet's now plot the mean value of the spectrums, but normalized by their total energy.  This will show us the mean shape of the spectrum, independent the total energy of the segments signals."},{"metadata":{"trusted":true,"_uuid":"1372388e77bcfec33aa3938ce3b1b42d1d6edfab"},"cell_type":"code","source":"normalized_spectrums = power_spectrums/(power_spectrums.sum(axis=1)[:,np.newaxis])\naverage_power_spectrum = pd.Series(normalized_spectrums.mean(axis=0))\naverage_power_spectrum = average_power_spectrum.rolling(window=10).mean()\n\nTOOLTIPS = [ (\"x\", \"@x{(0.00 a)}Hz\"), (\"y\", \"$y\")   ]\npl = Figure(plot_width=800, plot_height=400, title=\"Average of normalized power-spectrums\",\n            tooltips=TOOLTIPS)\npl.line(frequencies, average_power_spectrum, line_color=\"navy\")\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency [hz]\"\npl.yaxis.axis_label = \"Normalized Power []\"\nshow(pl)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a42fd2128ea4b16d099989ee3ef7f7f98bfac60f"},"cell_type":"markdown","source":"We see that on average, the power spectrum have 5 regions where the power spectrum reaches a local maximum:\n* Small maxima : 35-60 khz\n* Medium peak : 70-110 khz\n* Two contiguous maximums:  120-156khz and 156-200khz\n* Large maximum: 200-300khz"},{"metadata":{"_uuid":"3560b8bda93dedc87049b4198603a12cf3b437e8"},"cell_type":"markdown","source":"# Time to failure composites\n\nIn the last section, we show the average value and the average shape of the spectrum in the entire train dataset. That analysis masked out the dependence between the shape of spectrum shape and the remaining time to failure (quake time).\n\nIn this section, we explore the relationship between the spectrum shape with the remaining time to failure. For that, we will compute the normalize spectrum average over different \n"},{"metadata":{"trusted":true,"_uuid":"66711cb2a37a2fa3902b28cec193ac2582cbd6de","_kg_hide-input":false},"cell_type":"code","source":"# Create normalized spectrum composites for mean, std, min and max.\n\ndt = 0.2 # Interval of time to aggregate a composite\nmax_quake_time = quake_times.max()\nquake_times_intervals = np.arange(dt/2, max_quake_time+dt/2+0.01, dt)    \n\ncomposite_spectrums_mean = np.zeros((quake_times_intervals.size, frequencies.size))\ncomposite_spectrums_std = np.zeros((quake_times_intervals.size, frequencies.size))\ncomposite_spectrums_min = np.zeros((quake_times_intervals.size, frequencies.size))\ncomposite_spectrums_max = np.zeros((quake_times_intervals.size, frequencies.size))\n\nnormalized_spectrums = power_spectrums/(power_spectrums.sum(axis=1)[:,np.newaxis])\n\nfor n , quake_time in enumerate(quake_times_intervals):\n    t0 = quake_time-dt/2\n    t1 = quake_time+dt/2\n    segments = quake_times[(quake_times>=t0) & (quake_times<t1)].index.values\n    if len(segments) > 0:\n        composite_spectrums_mean[n] = normalized_spectrums[segments,:].mean(axis=0)\n        composite_spectrums_std[n] = normalized_spectrums[segments,:].std(axis=0)\n        composite_spectrums_min[n] = normalized_spectrums[segments,:].min(axis=0)\n        composite_spectrums_max[n] = normalized_spectrums[segments,:].max(axis=0)\n\nprint(\"composite_spectrums_mean.shape:\",composite_spectrums_mean.shape)\nprint(\"composite_spectrums_mean size:\",composite_spectrums_mean.nbytes/(1024**2.),\"[mb]\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8e3002dff82d46c18e44c915870b2d20d973ac88"},"cell_type":"markdown","source":"Next, we will make an interactive visualization of the mean normalized spectrum for each quake time interval."},{"metadata":{"trusted":true,"_uuid":"4723b4647712f642701ef220e026c749274e5074"},"cell_type":"code","source":"from bokeh.models import CustomJS, Slider\nfrom bokeh.plotting import figure, output_file, show, ColumnDataSource\n\n_quake_time = quake_times_intervals[0]\n\nTOOLTIPS = [ (\"x\", \"@frequencies{(0.00 a)}Hz\"), (\"y\", \"$y\")   ]\n\npl = Figure(plot_width=800, plot_height=400, title= f\"{_quake_time-dt/2} <= Quake time < {_quake_time+dt/2}\",\n           y_range=(0, 4e-3), tooltips=TOOLTIPS)\n\ndata_dict = dict()\nfor i in range(quake_times_intervals.size):    \n    data_dict[str(i)] = composite_spectrums_mean[i]\nall_data = ColumnDataSource(data=data_dict)\nsource = ColumnDataSource(data=dict(frequencies=frequencies, composite_mean=composite_spectrums_mean[0]))\n\npl.line('frequencies', 'composite_mean', line_color=\"navy\", source=source)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\n\n                            \ncallback = CustomJS(args=dict(source=source, \n                              all_data=all_data,plot=pl,\n                              quake_times_intervals=quake_times_intervals,\n                              dt=dt), \n                    code=\"\"\"\n    var data = source.data;\n    var interval = slider.value;\n    var composite_means = all_data.data[String(interval)];    \n    var y = data['composite_mean']\n    for (var i = 0; i < y.length; i++) {\n        y[i] = composite_means[i];\n    }\n    var _quake_time = quake_times_intervals[interval]\n    console.log(String(_quake_time-dt/2));\n    console.log(String(_quake_time+dt/2));\n    plot.title.text = String((_quake_time-dt/2).toFixed(1)) + \" <= Quake time < \" + String((_quake_time+dt/2).toFixed(1));\n    source.change.emit();\n\"\"\")\n\ntime_slider = Slider(start=0, end=composite_spectrums_mean.shape[0], value=0, step=1,\n                     title=\"Quaketime interval number\", callback=callback)\ncallback.args[\"slider\"] = time_slider\n                            \nlayout = column(time_slider,  pl)\nshow(layout)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f5794920c4233c1ed31422854d6a06570d52ad8d"},"cell_type":"markdown","source":"# Final remarks\n\nIn this notebook, we explored the Fourier spectrum characteristics. We found out that the main contributions come from 5 regions of the frequency spectrum. I have tried computing the common features over those frequency intervals but only marginal improvement over the usual computations was observed.\n\nIn addition,  the power spectra (normalized) have a similar shape for different quake times, suggesting that it is independent of the remaining time to failure. The motivation of this analysis was finding unusual trends or relationships that can be used as new features. Nevertheless, I couldn't find any interesting information to be used for new features. But probably I missed something.\n\nI hope that you found this kernel useful and or at least interesting.  \n\nAndrés"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}