{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Summary\nTrying to explore some features in the frequency domain.\n\n### Reference:\n* [EDA on the frequency domain (FFT)](https://www.kaggle.com/aperezhortal/eda-on-the-frequency-domain-fft)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","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%matplotlib inline\nimport matplotlib.pyplot as plt\n\nfrom scipy import signal # for fft\nfrom numpy import fft\nfrom tqdm import tqdm_notebook# progress bar","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# bokeh visualization for spectrum analysis\nfrom bokeh.layouts import column, row\nfrom bokeh.plotting import figure, show, output_file, ColumnDataSource\nfrom bokeh.io import output_notebook\nfrom bokeh.models import LinearAxis, Range1d, CustomJS, Slider, PrintfTickFormatter\n\noutput_notebook() # Display Bokeh plots inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain = pd.read_csv('../input/train.csv', \n                    dtype={'acoustic_data': np.int16, \n                           'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rows = 150000\nnum_segments = int(np.floor(train.shape[0] / rows))\ndt = 0.0375/rows\n# sampling_frequency is the upper bound of the frequency\nsampling_frequency = np.ceil(1/dt)\nprint(\"Total {0:d} segments, time step {1:.3e} seconds, sampling frequency {2:.0f}hz.\"\\\n      .format(num_segments, dt, sampling_frequency) )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Power spectrum estimation\n`scipy.signal`'s [periodogram](https://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.signal.periodogram.html) uses [Welch's method in the code](https://github.com/scipy/scipy/blob/v0.13.0/scipy/signal/spectral.py#L141) to estimate the power spectrum of a time signal, which is applying FFT on a real signal with a rolling averaging window. The function used in Welch's method is [fftpack.rfft](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.rfft.html).\n\nReference: P. Welch, \"The use of the fast Fourier transform for theestimation of power spectra: A method based on time averaging over short, modified periodograms\", IEEE Trans. Audio            Electroacoust. vol. 15, pp. 70-73, 1967."},{"metadata":{"trusted":true},"cell_type":"code","source":"# sampling_frequency is the upper bound of the frequency\nsampling_frequency = np.ceil(1/dt)\n\ndef get_spectrum(input_signal, window=10):\n    \"\"\"\n    returns the Fourier power spectrum for a given signal segment\n    output is a pandas Series \n    output.index is the frequencies\n    output.values is the amplitudes for each frequencies\n    default moving average window is 10\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(\n        input_signal, sampling_frequency, scaling='spectrum')    \n    \n    # Run a running windows average of 10-points to smooth the signal (default). \n    power_spectrum = pd.Series(power_spectrum, index=frequencies).rolling(window=window).mean()        \n    \n    return pd.Series(power_spectrum)\n\ndef get_segment_spectrum(segment_df, \n                         max_frequency = 300000, min_frequency = 100,\n                         step = 10, window = 10):\n    \"\"\"\n    get the Fourier power spectrum of a given segment.\n    returns the quake_time, frequencies, and power_spectrum\n    \"\"\"\n    \n    quake_time = segment_df['time_to_failure'].values[-1]\n#     quake_time = segment_df['time_to_failure'].values[-5000:].mean()\n    \n    _power_spectrum = get_spectrum(segment_df['acoustic_data'], window=window).dropna() \n    # drop the null values\n\n    # Keep only frequencies < max_frequency (larger frequencies have a negligible contribution).\n    # and > min_frequency (some low frequency is outlier)\n    _power_spectrum = _power_spectrum[_power_spectrum.index < max_frequency]\n    _power_spectrum = _power_spectrum[_power_spectrum.index > min_frequency]\n    \n    # Keep one every 10 samples by default\n    power_spectrum=_power_spectrum.values[::step]\n    frequencies=_power_spectrum.index.values[::step]    \n    \n    return quake_time, frequencies, power_spectrum","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Some examples"},{"metadata":{},"cell_type":"markdown","source":"## Example around 1st earthquake \nExplanations of the functions above, first the frequency power spectrum for this segment is obtained through\n* Smoothened data using a 10-points rolling window;\n* Keeping one every 10 samples;\n* Keeping only frequencies < 350khz and > 100hz.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# this is the segment that has the first earth quake, length = 150000\nfirst_quake_time = 5656574\nwindow = 10\nsegment = train.iloc[first_quake_time-rows//2:first_quake_time+rows//2]\nx_tr = np.asarray(segment['acoustic_data'].values, dtype='float64')\nx_tr -= x_tr.mean()\n\nfrequencies , power_spectrum_1 = signal.periodogram(x_tr, sampling_frequency, scaling='spectrum')    \n\npower_spectrum_1 = pd.Series(power_spectrum_1, index=frequencies)\npower_spectrum_1 = power_spectrum_1.rolling(window=window).mean().dropna()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# for this segment the time_to_failure\nmax_frequency = 350000\nquake_time = segment['time_to_failure'].values[-1]\npower_spectrum_1 = power_spectrum_1[power_spectrum_1.index < max_frequency]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore the power spectrum of the first segment\n\nReference: [Bokeh gallery](https://bokeh.pydata.org/en/latest/docs/gallery.html)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Bokeh tooltips\n# \"@foo{(0.00 a)}\" # formats 1230974 as: 1.23 m\nTOOLTIPS = [ (\"x\", \"@x{(0 a)}Hz\"), (\"y\", \"$y\")   ]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"freq1 = power_spectrum_1.index\nspec1 = power_spectrum_1.values\n\npl1 = figure(plot_width=750, plot_height=400,title=\"Power spectrum of 1 segment\",tooltips=TOOLTIPS)\npl1.line(freq1, spec1)\npl1.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl1.xaxis.axis_label = \"Frequency\"\npl1.yaxis.axis_label = \"Power (L^2 norm squared)\"\nshow(pl1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"window_tmp = 100\nfreq2 = power_spectrum_1.index\nspec2 = power_spectrum_1.rolling(window=window_tmp).mean()\n\npl2 = figure(plot_width=750, plot_height=400,\n             title=f\"Power spectrum of 1 segment (sampling {100/step}% data)\",tooltips=TOOLTIPS)\npl2.line(freq2, spec2)\npl2.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl2.xaxis.axis_label = \"Frequency\"\npl2.yaxis.axis_label = \"Power (L^2 norm squared)\"\nshow(pl2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Retrieving the power spectrum for all segments"},{"metadata":{"trusted":true},"cell_type":"code","source":"## perform above for all segments\n# since the frequencies retrieved is changing, initialization is none\nwindow = 20\nstep = 5\nquake_times = np.array([])       \npower_spectrums = np.array([])\nmax_frequency = 350000\nmin_frequency = 100\n\nfor k in tqdm_notebook(range(num_segments)):\n    segment = train.iloc[k*rows : k*rows+rows]\n    quake_time, frequencies, power_spectrum = \\\n    get_segment_spectrum(segment, \n                         max_frequency=max_frequency,\n                         min_frequency=min_frequency,\n                         step=step,window=window)    \n    quake_times = np.append(quake_times, quake_time)\n    power_spectrums = np.append(power_spectrums, power_spectrum)\n    \npower_spectrums = power_spectrums.reshape(num_segments,-1)\nquake_times = quake_times.reshape(num_segments,-1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_spectrums.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"frequencies.size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"window = 5\n\nnormalized_spectrums = power_spectrums/(np.sum(power_spectrums, axis=1)[:,np.newaxis])\n# since power_spectrums are already L^2 norm squared\n# normalization is done for each row by dividing its total energy for total 4194 rows\n\naverage_power_spectrum = pd.Series(np.mean(normalized_spectrums,axis=0)) # average all 4194 spectrum\naverage_power_spectrum = average_power_spectrum.rolling(window=window).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pl = figure(plot_width=750, plot_height=400,title=f\"Normalized average power spectrum with rolling window {window}\",\n            tooltips=TOOLTIPS)\npl.line(frequencies, average_power_spectrum)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency\"\npl.yaxis.axis_label = \"Normalized Power\"\nshow(pl)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"On average, the normalized energy has 5 local maximum around 5 peaks among the following frequency ranges.\n\nHovering the mouse on the diagram can reveal the values of the coordinates.\n\n* 45-60 khz\n* 80-100 khz\n* 130-150khz \n* 165-185khz\n* 230-270khz"},{"metadata":{"trusted":true},"cell_type":"code","source":"freq_intervals = np.array([0, 7e3, 5e4, 1e5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"frequencies.max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X_fft = pd.DataFrame(index=range(num_segments), dtype=np.float64)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for seg_id in tqdm_notebook(range(num_segments)):\n    train_X_fft.loc[seg_id, 'energy_mean_in_7khz'] = \\\n    normalized_spectrums[seg_id, frequencies<7e3].mean()\n    train_X_fft.loc[seg_id, 'energy_mean_in_50khz_100khz'] = \\\n    normalized_spectrums[seg_id, (frequencies>5e4) & (frequencies<1e5)].mean()\n    train_X_fft.loc[seg_id, 'energy_mean_in_130khz_150khz'] = \\\n    normalized_spectrums[seg_id, (frequencies>1.3e5) & (frequencies<1.5e5)].mean()\n    train_X_fft.loc[seg_id, 'energy_mean_in_165khz_185khz'] = \\\n    normalized_spectrums[seg_id, (frequencies>1.65e5) & (frequencies<1.85e5)].mean()\n    train_X_fft.loc[seg_id, 'energy_mean_in_230khz'] = \\\n    normalized_spectrums[seg_id, (frequencies>2.3e5)].mean()\n    \n    train_X_fft.loc[seg_id, 'energy_q90_in_7khz'] = \\\n    np.quantile(normalized_spectrums[seg_id, frequencies<7e3], 0.9)\n    train_X_fft.loc[seg_id, 'energy_q90_in_50khz_100khz'] = \\\n    np.quantile(normalized_spectrums[seg_id, (frequencies>5e4) & (frequencies<1e5)], 0.9)\n    train_X_fft.loc[seg_id, 'energy_q90_in_130khz_150khz'] = \\\n    np.quantile(normalized_spectrums[seg_id, (frequencies>1.3e5) & (frequencies<1.5e5)], 0.9)\n    train_X_fft.loc[seg_id, 'energy_q90_in_165khz_185khz'] = \\\n    np.quantile(normalized_spectrums[seg_id, (frequencies>1.65e5) & (frequencies<1.85e5)], 0.9)\n    train_X_fft.loc[seg_id, 'energy_q90_in_230khz'] = \\\n    np.quantile(normalized_spectrums[seg_id, (frequencies>2.3e5)], 0.9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_freq = train_X_fft.shape[1]\n_, ax1 = plt.subplots(num_freq,1, figsize=(15,5*num_freq))\nax1 = ax1.reshape(-1)\n\nfor i, col in enumerate(train_X_fft.columns):\n    ax1[i].plot(train_X_fft[col], color='blue')\n    ax1[i].set_title(\"Feature: {:s} hz\".format(col))\n    ax1[i].set_ylabel(\"Energy\", color='b')\n\n    ax2 = ax1[i].twinx()\n    ax2.plot(quake_times, color='g', linewidth = 2)\n    ax2.set_ylabel('time_to_failure', color='g')\n    ax2.grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualizing energy vs time_to_failure for all segments\n\nWe try to explore `time_to_failure`'s relation to the spectrum. First `time_to_failure` is divided into intervals. We look at the average spectrum for the segments of which the `time_to_failure` falls into the range. "},{"metadata":{"trusted":true},"cell_type":"code","source":"Dt = dt*(50085878-5656574)/20  # 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/4)\nnum_freq = frequencies.size\nnum_intervals = quake_times_intervals.size\nprint(\"The compositite spectrum is analyzed for %d time_to_failure intervals, has %d frequencies.\" \n      %(num_intervals,num_freq))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"composite_spectrums_mean = np.zeros((num_intervals, num_freq))\ncomposite_spectrums_std = np.zeros((num_intervals, num_freq))\ncomposite_spectrums_min = np.zeros((num_intervals, num_freq))\ncomposite_spectrums_max = np.zeros((num_intervals, num_freq))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\nThe normalized spectrum can be used as features of all the segments."},{"metadata":{"trusted":true},"cell_type":"code","source":"for n , quake_time in enumerate(quake_times_intervals):\n#     t0 = quake_time-Dt/2\n#     t1 = quake_time+Dt/2\n    t0 = quake_time-Dt*0.75 # overlapping\n    t1 = quake_time+Dt*0.75\n    segments = ((quake_times>=t0) & (quake_times<t1)).reshape(-1)\n    if len(segments) > 0:\n        # averaging the spectra for those who fall into the same quake time interval\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        \n_quake_time = quake_times_intervals[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pl = figure(plot_width=800, plot_height=400, \n            title= f\"Avg spectra for {_quake_time-Dt*3/4} <= time_to_failure < {_quake_time+Dt*3/4}\",\n            y_range=(0, 3e-3), tooltips=TOOLTIPS)\n\ndata_dict = dict()\nfor i in range(quake_times_intervals.size):    \n    data_dict[str(i)] = composite_spectrums_mean[i]\n\nall_data = ColumnDataSource(data=data_dict)\nsource = ColumnDataSource(data=dict(frequencies=frequencies, \n                                    composite_mean=composite_spectrums_mean[0]))\n\npl.line('frequencies', 'composite_mean', source=source)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency\"\npl.yaxis.axis_label = \"Normalized Power (L^2-norm squared)\"\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*3/4));\n    console.log(String(_quake_time+Dt*3/4));\n    plot.title.text = String(\"Avg spectra for \"+(_quake_time-Dt*3/4).toFixed(3)) + \" <= time_to_failure < \" + String((_quake_time+Dt*3/4).toFixed(3));\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":{},"cell_type":"markdown","source":"# The power spectrum for subdividing a training segment\n\nDividing the several training segments (150000 time steps) into sub-segments."},{"metadata":{"trusted":true},"cell_type":"code","source":"# for random K samples\nK = 5\nseg_spectrums = np.array([])\n# seg_quake_times = np.array([])\nsegs = np.random.choice(num_segments, size=K, replace=True)\n\nwindow = 10\nsampling_step = 10\nnum_subsegs = 10\nlen_subsegs = int(rows/num_subsegs)\nfor n in tqdm_notebook(segs):\n    for k in range(num_subsegs):\n        sub_seg = train.iloc[n*num_subsegs+k*len_subsegs : n*num_subsegs+k*len_subsegs+len_subsegs]\n        _, frequencies, seg_spectrum = get_segment_spectrum(sub_seg, \n                                                            max_frequency=max_frequency,\n                                                            step=sampling_step, window=window)\n        seg_spectrums = np.append(seg_spectrums, seg_spectrum)\n#         seg_quake_times = np.append(seg_quake_times, seg_quake_time)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seg_spectrums = seg_spectrums.reshape(K*num_subsegs,-1)\nnormalized_seg_spectrums = seg_spectrums/(np.sum(seg_spectrums, axis=1)[:,np.newaxis])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seg_quake_times = train['time_to_failure'].values[rows::rows]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(seg_quake_times.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"average_seg_spectrum = pd.Series(np.mean(normalized_seg_spectrums,axis=0)) \n# average all 10 sub-segments\naverage_seg_spectrum = average_seg_spectrum.rolling(window=window).mean()\n\npl = figure(plot_width=750, plot_height=400,\n            title=f\"Average power spectrum with rolling window {window} for Segment {segs[0]}\",\n            tooltips=TOOLTIPS)\npl.line(frequencies, average_seg_spectrum)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency\"\npl.yaxis.axis_label = \"Normalized Power\"\nshow(pl)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(normalized_seg_spectrums.shape)\nnormalized_seg_spectrums = normalized_seg_spectrums.reshape(K,num_subsegs,-1)\nprint(normalized_seg_spectrums.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## plot the subsegment power spectrum\nsubseg_intervals = np.arange(0, 150000+0.01, len_subsegs)\n\npl = figure(plot_width=750, plot_height=400, \n            title= f\"Training Seg No {segs[0]}: normalized spectrum for {0} <= time step <= {len_subsegs}\",\n            y_range=(0, 1e-1), tooltips=TOOLTIPS)\n\ndata_dict = dict()\nfor i in range(num_subsegs):    \n    data_dict[str(i)] = normalized_seg_spectrums[0,i,:]\n\nall_data = ColumnDataSource(data=data_dict)\nsource = ColumnDataSource(data=dict(frequencies=frequencies, \n                                    seg_spectrum=normalized_seg_spectrums[0]))\n\npl.line('frequencies', 'seg_spectrum', source=source)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency\"\npl.yaxis.axis_label = \"Normalized Power\"\n\n                            \ncallback = CustomJS(args=dict(source=source, \n                              all_data=all_data,plot=pl,\n                              subseg_intervals=subseg_intervals,\n                              len_subsegs=len_subsegs), \n                    code=\"\"\"\n    var data = source.data;\n    var interval = slider.value;\n    var seg_spectrum = all_data.data[String(interval)];    \n    var y = data['seg_spectrum']\n    for (var i = 0; i < y.length; i++) {\n        y[i] = seg_spectrum[i];\n    }\n    var _seg_time = subseg_intervals[interval]\n    console.log(String(_seg_time));\n    console.log(String(_seg_time+len_subsegs));\n    plot.title.text = String(\"Training Seg No.1: Normalized spectrum for \"+(_seg_time).toFixed(0)) + \" <= time step <= \" + String((_seg_time+len_subsegs).toFixed(0));\n    source.change.emit();\n\"\"\")\n\ntime_slider = Slider(start=0, end=num_subsegs, value=0, step=1,\n                     title=\"Sub-segment number\", callback=callback)\ncallback.args[\"slider\"] = time_slider\n                            \nlayout = column(time_slider,  pl)\nshow(layout)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# The power spectrum for one testing segment\n\nDividing the first testing segment (150000 time steps) into 10 sub-segments."},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', \n                         index_col='seg_id', \n                         dtype={\"time_to_failure\": np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(submission)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seg_id = submission.index[np.random.randint(len(submission))]\nseg_test = pd.read_csv('../input/test/' + seg_id + '.csv')\nseg_test['time_to_failure'] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# a random testing sample\nseg_spectrums = np.array([])\nwindow = 10\nsampling_step = 10\nnum_subsegs = 10\nlen_subsegs = int(rows/num_subsegs)\n\nfor k in tqdm_notebook(range(num_subsegs)):\n    sub_seg = seg_test.iloc[k*len_subsegs : k*len_subsegs+len_subsegs]\n    _, frequencies, seg_spectrum = get_segment_spectrum(sub_seg, \n                                                        max_frequency=max_frequency,\n                                                        step=sampling_step, window=window)\n    seg_spectrums = np.append(seg_spectrums, seg_spectrum)\n    \nseg_spectrums = seg_spectrums.reshape(num_subsegs,-1)\nnormalized_seg_spectrums = seg_spectrums/(np.sum(seg_spectrums, axis=1)[:,np.newaxis])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"normalized_seg_spectrums.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"frequencies.size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"average_seg_spectrum1 = pd.Series(np.mean(normalized_seg_spectrums[0],axis=0)) \n# average all 10 sub-segments for first sample\naverage_seg_spectrum1 = average_seg_spectrum.rolling(window=window).mean()\naverage_seg_spectrum1.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pl = figure(plot_width=750, plot_height=400,\n            title=f\"Normalized average spectrum with rolling window {window} for segment {seg_id}\",\n            tooltips=TOOLTIPS)\npl.line(frequencies, average_seg_spectrum1)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency\"\npl.yaxis.axis_label = \"Normalized Power\"\nshow(pl)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## plot the subsegment power spectrum\nsubseg_intervals = np.arange(0, 150000+0.01, len_subsegs)\n\npl = figure(plot_width=750, plot_height=400, \n            title= f\"Testing Seg No.1: normalized spectrum for {0} <= time step <= {len_subsegs}\",\n            y_range=(0, 1e-1), tooltips=TOOLTIPS)\n\ndata_dict = dict()\nfor i in range(num_subsegs):    \n    data_dict[str(i)] = normalized_seg_spectrums[i]\n\nall_data = ColumnDataSource(data=data_dict)\nsource = ColumnDataSource(data=dict(frequencies=frequencies, \n                                    seg_spectrum=normalized_seg_spectrums[0]))\n\npl.line('frequencies', 'seg_spectrum', source=source)\npl.xaxis[0].formatter = PrintfTickFormatter(format=\"%d\")\npl.xaxis.axis_label = \"Frequency\"\npl.yaxis.axis_label = \"Normalized Power\"\n\n                            \ncallback = CustomJS(args=dict(source=source, \n                              all_data=all_data,plot=pl,\n                              subseg_intervals=subseg_intervals,\n                              len_subsegs=len_subsegs), \n                    code=\"\"\"\n    var data = source.data;\n    var interval = slider.value;\n    var seg_spectrum = all_data.data[String(interval)];    \n    var y = data['seg_spectrum']\n    for (var i = 0; i < y.length; i++) {\n        y[i] = seg_spectrum[i];\n    }\n    var _seg_time = subseg_intervals[interval]\n    console.log(String(_seg_time));\n    console.log(String(_seg_time+len_subsegs));\n    plot.title.text = String(\"Testing Seg No.1: Normalized spectrum for \"+(_seg_time).toFixed(0)) + \" <= time step <= \" + String((_seg_time+len_subsegs).toFixed(0));\n    source.change.emit();\n\"\"\")\n\ntime_slider = Slider(start=0, end=num_subsegs, value=0, step=1,\n                     title=\"Sub-segment number\", callback=callback)\ncallback.args[\"slider\"] = time_slider\n                            \nlayout = column(time_slider,  pl)\nshow(layout)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}