{"cells":[{"metadata":{"_uuid":"fb4817148314369eb9077574693fe790ce926813"},"cell_type":"markdown","source":"# Rock Music! Lomb-Scargle Periodograms for HF Noise Spectral Analysis\n---\n\n\nAs this competition will hinge largely on signal analysis, I thought I would make a feature engineering contribution that tries to isolate the main constituent frequencies amid the noise. \n\nFrequency analysis is typically performed with the [Fast Fourier Transform](http://mathworld.wolfram.com/FastFourierTransform.html). This breaks down a complex waveform into its sinusoidal constituents, mapping them into freqency space with their respective amplitudes:\n\n![](https://cdn.iopscience.com/images/books/978-1-627-05419-5/live/bk978-1-627-05419-5ch8f1_online.jpg)\n\nFourier transforms act on continuous periodic functions. Since real-life data is discretised, numerical solutions have to be found via methods like DFFT (Discrete Fast Fourier Transform). However, calculating this requires a time axis with regular intervals - that is, the time difference between successive amplitude measurements is constant. Many real-life situations, along with this competition, require frequency analysis on data where the time between measurements is unequal. What's the solution? The fantasically named **Lomb-Scargle Periodogram**:\n\n![](https://static.packt-cdn.com/products/9781785282287/graphics/B04223_06_11.jpg)\n\nIf this is new to you, please take a moment to say it again out loud. Great stuff! LSP is one approach among many for least-squares spectral analysis problems, and I'm sure other kagglers will find better solutions. But it's early days in the competition so hopefully this will help you get started! For a comprehensive primer on LSP, I strongly recommend [reading this article (PDF)](https://arxiv.org/pdf/1703.09824.pdf).\n\nThe test data doesn't include the time axis since that is our target, so when adding features you'll have to assume the time difference between measurements is even. With train.csv we can't make that assumption, which is why I tried this approach.\n\nFor speed, I'll be running this kernel on the data for the first earthquake period in train.csv only. I encourage other readers to expand on it for their own work - you might stumble on something critical to your final performance! As always, the first step is to read in the data and have a look at our variables:"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing\nimport matplotlib.pyplot as plt\nimport gc\n\n#thanks to EEK: https://www.kaggle.com/theupgrade/creating-breaks-as-start-ends-of-bins\nttf_resets = np.array([5656575,50085879,104677357,138772454,187641821,\n                  218652631,245829586,307838918,338276288,375377849,\n                  419368881,461811624,495800226,528777116,585568145,\n                  621985674]) - 1\n\ntrain = pd.read_csv('../input/train.csv', nrows=ttf_resets[0])\n#rename columns - quicker to type!\ntrain.columns = ['signal', 'ttf']\nprint(train.describe())\ntrain.plot(kind='scatter',x='ttf',y='signal', s=1)\nplt.xlim(1.5, -0.05)\nplt.title('Acoustic signal/ttf')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b68eae8ef0c63aa456bf9caa7af1686ddc4cca96"},"cell_type":"markdown","source":"The amplitude spikes during certain events. Periodic trends are difficult to identify visually. Taking a look at the signal amplitude distribution alone:"},{"metadata":{"trusted":true,"_uuid":"e4a9f67db5c88f311f2a5c008703e8fffc0360f6","_kg_hide-input":true},"cell_type":"code","source":"train.signal.hist(bins=200)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f20c44a66004319a0bab177360452c1c72b94fb5"},"cell_type":"markdown","source":"A histogram of the entire signal data isn't particularly helpful since by nature, notable seismic activity manifests as high-amplitude burst phenomena. Removing the extreme outliers at TTF ~ 0.35s, we can get a clearer view of the majority of the data:"},{"metadata":{"trusted":true,"_uuid":"8d772916a3b389e4f2b03f5e58e319e9c1023fbd"},"cell_type":"code","source":"train.loc[(train.signal < 40) & (train.signal > -40), ].signal.hist(bins=200)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bde797bafe752ff7d2057b8c9c6df258f635e0b1"},"cell_type":"markdown","source":"An ordinary, symmetric distribution with a mean of around 4.56. I'm going to avoid mean normalisation since the general form of LSP doesn't require it, and a mean value != 0 is a feature in itself. Looking at the measurement times: "},{"metadata":{"trusted":true,"_uuid":"5b69b7b7536bfd99a359dd2351f00f4af965d8d4","_kg_hide-input":true},"cell_type":"code","source":"train.ttf.hist()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"87d6ba25dbe17b67ffd242229fa0b29429b9ee86"},"cell_type":"markdown","source":"As expected, we have a mostly - but not perfectly - even time distribution. Python's inbuilt rounding limits how accurately we can see the see the time difference between data points but we can inspect them with the `decimal` package:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"a3294dfb518719dbe95e9cf2874a70e05cec5622"},"cell_type":"code","source":"import decimal as dec\nt1 = dec.Decimal(train.iloc[0, 1])\nt2 = dec.Decimal(train.iloc[1, 1])\nt3 = dec.Decimal(train.iloc[2, 1])\nprint('First TTF: ', t1, '(s)')\nprint('Second TTF: ', t2, '(s)')\nprint('Third TTF: ', t3, '(s)')\n\nprint('First delta t: ', t1-t2, '(s)')\nprint('Second delta t: ', t2-t3, '(s)')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6b82222a4efe585b6a6587d1a59a2dbfaa2fb1b5"},"cell_type":"markdown","source":"Delta-t for these adjacent readings is approximately 1.1E-9, meaning the [Nyquist Frequency](http://mathworld.wolfram.com/NyquistFrequency.html) - the upper bound on the frequency we can reliably detect - would be about 450MHz. For audio data, this is an extremely high frequency.\n\nHigh-frequency waves are readily attenuated by rock, and actual seismological waves generally lie in the range of 0.1-1Hz since they travel further and hence have more predictive and forensic value. This experiment was at a much smaller scale however, so analysis of higher frequencies may yield valuable insights. So what is the *minimum* frequency we can reliably measure?\n\nThis poses something of a conundrum. The data is time-dependant, and ideally we'd like to have a number of LSP frequency spectra to examine as TTF approaches 0. This way, any changes in the signal profile over time can be examined - this is a prediction challenge after all. But the more samples we take, but smaller the timespan they will cover, and hence the lower bound for our frequency range gets higher. For this kernel I'm naively choosing a range of 10000 rows to get 566 different readings for the first earthquake, equivalent to around 0.0021 seconds per segment. For this the minimum frequency we could expect to differentiate would be approximately 1/0.021s or ~ 475Hz.\n\nFor future kernels I'll try a sliding-window approach with a larger number of samples to increase the frequency resolution and examine if low-frequencies have any notable trends in the training data. Be warned if doing this yourself as this is a very memory intensive process - running LSP on more than a couple of millions rows at once will give you the dreaded Memory Error! To start let's analyse the first sample alone:"},{"metadata":{"trusted":true,"_uuid":"a069bf4602d4d711dd5e90692443c8c63925433a"},"cell_type":"code","source":"from astropy.stats import LombScargle\n\nt = train.ttf.values\ny = train.signal.values\n\n#astropy implementation works directly on frequency, not angular frequency\n#autopower() calculates a frequency grid automatically based on mean time separation\nfrequency, power = LombScargle(t[0:10000], y[0:10000]).autopower(nyquist_factor=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"45356e5541d188f042ca2d65dbb47a637e438f29"},"cell_type":"code","source":"plt.plot(frequency, power) \nplt.title('Noise LSP Frequency-Amplitude Spectrum')\nplt.xlabel('Frequency (Hz)')\nplt.ylabel('Amplitude')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"db09ce2bb061ef28601d3e890ea0c5269b33efc7"},"cell_type":"markdown","source":"Overall the frequencies are highly smeared, characteristic of noisy data collected over time. A number of peaks are present but the amplitude is very low, far lower than anything one would consider to be statistically important. This is to be expected since this segment of the data doesn't display any notable acoustic activity. Taking a closer look at the initial peaks:"},{"metadata":{"trusted":true,"_uuid":"9889be1f8267c5ff2c31347ddc6197af4ff4c6c9","_kg_hide-input":true},"cell_type":"code","source":"#power multiplied for graphical visibility\nfreq_df = pd.DataFrame({'freq': frequency.round(),\n                       'amp': (power*1e6).round(1)})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7a106bbf87e3c7003269aa03795eab6c3d6ef4c9","_kg_hide-input":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfreq_df.loc[freq_df.freq < 1000000, :].plot(kind='scatter',x='freq',y='amp', s=2)\nplt.title('Noise LSP Frequency-Amplitude Spectrum')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"70547f562ac2f0e695eb312fe2090dc3608079ec","_kg_hide-input":true},"cell_type":"code","source":"freq_df.loc[freq_df.freq < 150000, :].plot(kind='scatter',x='freq',y='amp', s=2)\nplt.title('Noise LSP Frequency-Amplitude Spectrum')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7a266299e6750b6d2f6421873cfc3d343a3ac281"},"cell_type":"markdown","source":"In frequency space, the acoustic noise displays sinusoidal trends along the 10kHz harmonics that are likely artefacts of the frequency grid created by `LombScargle`. Examining the frequencies sorted by amplitude:"},{"metadata":{"trusted":true,"_uuid":"4bf043931058887a910d0262f085ffeee0e7cda7","_kg_hide-input":true},"cell_type":"code","source":"freq_df.sort_values('amp', ascending=False).head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"121fe274221963e13eb6d26a20f223bba6aadf85"},"cell_type":"markdown","source":"From the first 10 results, values that lie on the same frequency 'hump' can be seen with values close to eachother, which would yield unhelpful results. The initial method I'll be using to examine this data is as follows:\n\n* Run LSP for each segment of 10000 rows up to TTF==0\n* Remove frequencies with an associated period smaller than the time window of these rows\n* Remove frequencies higher than a heuristically chosen value (optional)\n* Extract the 5 frequencies with the highest amplitudes **and** a minium separation of 0.2MHz and add them to train\n* Plot how these frequencies change in the build up to the quake\n* Examine how the extracted frequencies could be used to engineer new features \n\nThis approach has its flaws. Notably, [as examined here](https://arxiv.org/pdf/1703.09824.pdf), taking the highest-amplitude frequency from LSP is a naive means of interpreting the data. The frequencies are highly spread out, so to remove values which are clustered on the same peak, adjacent frequency recordings have to be at least 0.2MHz apart. This corresponds to the approximate distance between 'humps'. Ideally we'd like to obtain one frequency for each 'hump' in the frequency spectrum, as are clearly visible below:"},{"metadata":{"trusted":true,"_uuid":"b12df61f006e0fb44164f428a9502f30012ee2b5","_kg_hide-input":true},"cell_type":"code","source":"freq_df.loc[(freq_df.freq > 125000) & (freq_df.freq < 1000000), :].plot(kind='scatter',x='freq',y='amp', s=2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"40cafef47d66ff5a399d3f659c0a243736e9084a"},"cell_type":"markdown","source":"The quick and dirty function `arr_dist()` will run across an array and keep only values whose separation is greater than 0.2MHz until it has 5 values. The frequencies recorded here will be placed in 5 new columns in train. If expanding this process to other sections of the data, you may want to decrease the number of new columns since this particular segment of the data is rather small, and is less intensive to work with than most of the earthquakes.\n\nThis function is rather long but quite straightforward and easily modifiable. If you have any questions, I'll be happy to help in the comments."},{"metadata":{"trusted":true,"_uuid":"458fd08ae272f404c6602c6ed03713d6e84dece8"},"cell_type":"code","source":"def arr_dist(arr, sep, n=5):\n    output = []\n    for x in arr:\n        keep=True\n        for y in output:\n            if abs(y-x)<sep:\n                keep=False\n                break\n        if(keep):\n            output.append(x)\n            if len(output)==n:\n                return(np.asarray(output))\n\nROWS_PER_SEGMENT = 10000\nTIME_PER_SEGMENT = train.iloc[0, 1] - train.iloc[ROWS_PER_SEGMENT, 1]\nMIN_FREQ = round(1/TIME_PER_SEGMENT)\nMAX_FREQ = 1e8\nFREQ_SEP = 0.2e6\n\ndef LSP_freq(df, signal_col='signal', time_col='ttf',\n             nrows=ROWS_PER_SEGMENT, min_freq=MIN_FREQ, max_freq=MAX_FREQ):\n    print('Lomb-Scargle Periodogram analysis commencing.')\n    print('Minimum detection frequency: {}Hz'.format(MIN_FREQ))\n    print('Manual maximum frequency cutoff: {}Hz'.format(MAX_FREQ))\n    print('Number of segments: ', round(len(df)/ROWS_PER_SEGMENT))\n    #initialise empty arrays for frequency outputs to be concatenated to DataFrame \n    freq_1 = np.zeros(len(df))\n    freq_2 = np.zeros(len(df))\n    freq_3 = np.zeros(len(df))\n    freq_4 = np.zeros(len(df))\n    freq_5 = np.zeros(len(df))\n    segment_num = np.zeros(len(df))\n    #loop through input DataFrame in chunks of length=nrows\n    init_id = 0\n    segment_id =1\n    while init_id < len(df):\n        if segment_id==1:\n            print('Processing segment {:d}...'.format(segment_id))\n        if segment_id%25==0:\n            print('Processing segment {:d}...'.format(segment_id))\n        end_id = min(init_id + nrows, len(df))\n        ids = range(init_id, end_id)\n        df_chunk = df.iloc[ids]\n        #np arrays of amplitude and time columns\n        signal = df_chunk[signal_col].values\n        ttf = df_chunk[time_col].values\n        #clear memory\n        del df_chunk\n        gc.collect()\n        #calulate Lomb-Scargle periodograms for spectral analysis\n        freq, amp = LombScargle(ttf, signal).autopower(nyquist_factor=2)\n        freq_df = pd.DataFrame({'freq': freq.round(),\n                               'amp': amp})\n        #obtain frequencies sorted by highest amplitude as np.array\n        top_freqs = freq_df.loc[(freq_df.freq > min_freq) & (freq_df.freq < max_freq)].sort_values('amp', ascending=False).freq.values\n        del freq_df, freq, amp\n        gc.collect()\n        #obtain top 5 values that do not lie within 1kHz of eachother\n        top_freqs = arr_dist(top_freqs, sep=FREQ_SEP)\n        #sort principal frequencies from highest to lowest\n        top_freqs = -np.sort(-top_freqs)\n        #update main frequency component arrays\n        freq_1[ids] = top_freqs[0]\n        freq_2[ids] = top_freqs[1]\n        freq_3[ids] = top_freqs[2]\n        freq_4[ids] = top_freqs[3]\n        freq_5[ids] = top_freqs[4]\n        segment_num[ids] = segment_id\n        del top_freqs\n        init_id += nrows\n        segment_id += 1\n    print('...Done. Adding main component frequencies as DataFrame columns...')\n    df['Freq_1'] = freq_1\n    df['Freq_2'] = freq_2\n    df['Freq_3'] = freq_3\n    df['Freq_4'] = freq_4\n    df['Freq_5'] = freq_5\n    df['Segment'] = segment_num\n    df['Freq_MinMax'] = df['Freq_1'] - df['Freq_5']\n    print('...Done.')\n    \nLSP_freq(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"54600a1a24795fcb135918cafa2cd4ac6a43ad53"},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"80460ca45b49bdb5cd5b252e785752245e0be598"},"cell_type":"markdown","source":"Now we can examine the noise component frequencies over time, smoothing them via a rolling-mean:"},{"metadata":{"trusted":true,"_uuid":"85a2fcbb192c4075d0cbbd2390c3c7f6e13e4894"},"cell_type":"code","source":"cols = ['Freq_1', 'Freq_2','Freq_3','Freq_4','Freq_5']\nfor x in cols:\n    train[x + '_RM'] = train[x].rolling(window=100000,center=False).mean()\ncols_rm = []\nfor x in range(0, 5):\n    cols_rm.append(cols[x] + '_RM')\n\nax = plt.gca()\ntrain.plot(kind='line',x='ttf',y=cols_rm ,ax=ax, figsize=(12, 6))\nplt.xlim(1.5, -0.05)\nplt.title('5 Main Noise Component Frequencies (LSP, min separation 200kHz)')\nplt.xlabel('TTF(s)')\nplt.ylabel('Frequency (Hz)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"379ed9d07a41399c90e37a754e680945cbcfc7e6"},"cell_type":"markdown","source":"There's a lot to examine here, and certainly many parameters to change that may lead to more informative results. The main frequencies spike during the earthquake event, but this is not the only time they do so. This data is independent of the signal amplitude, and seems to be sensitive to smaller jumps in amplitude that also cause high-frequency noise. "},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"9925848c432c782a375e359a39b68497be409e35"},"cell_type":"code","source":"train['FREQ_MEAN'] = train[cols].mean(axis=1)\ntrain['FREQ_MEAN_RM'] = train['FREQ_MEAN'].rolling(window=250000, center=False).mean()\ntrain.plot(kind='line',x='ttf',y='FREQ_MEAN_RM', figsize=(12, 6))\nplt.xlim(1.5, -0.05)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"07ab3b47d506544bab755e502a8970a041733ca8"},"cell_type":"code","source":"train['FREQ_RATIO'] = train['Freq_1']/train['Freq_2']\ntrain['FREQ_RATIO_RM'] = train['FREQ_RATIO'].rolling(window=500000, center=False).mean()\ntrain.plot(kind='line',x='ttf',y='FREQ_RATIO_RM', figsize=(12, 6))\nplt.xlim(1.5, -0.05)\nplt.title('Ratio of Highest Main Frequency/Second Highest: Rolling-Mean')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b36f4266dbdec31a2ae00afe2674a314473f8a19"},"cell_type":"markdown","source":"Determining trends visually is difficult, which is no surprise since the data by nature is both stochastic and noisy. Also, will identified trends be reapeated in the data for other earthquakes?\n\nTo isolate instances where the frequency detected by LSP is statistically viable, a typical approach is a false-alarm filter on the results, from [Scargle's original 1982 paper](https://www.researchgate.net/profile/Thomas_Ruf3/publication/245535651_The_Lomb-Scargle_Periodogram_in_Biological_Rhythm_Research_Analysis_of_Incomplete_and_Unequally_Spaced_Time-Series/links/02e7e51d72fa498ab1000000.pdf): ![](https://i.imgur.com/O0ByWdL.png)  \n\nThis approach has not been successful during my work, and has been [identified as an inappropriate method](https://github.com/astropy/astropy/issues/7618) for assessing the viability of identified frequencies. In any case, the nature of this dataset limits the occurance of statistically notable frequencies. Instead I'll be using an amplitude cut-off which displays the occurances of frequencies that, while not having a p-value < 0.05, still stand out from the background noise."},{"metadata":{"trusted":true,"_uuid":"e5569f7828a2c43552626f136ab12cf6a7c43a44"},"cell_type":"code","source":"#false-alarm filter - unused in this kernel\ndef LSP_filter(p, M):\n    threshold = -np.log(1-((1-p)**(1/M)))\n    return(threshold)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"80823fbfd5e772fcf2c60cd543c37e024c441b55"},"cell_type":"code","source":"ROWS_PER_SEGMENT = 1000\nTIME_PER_SEGMENT = train.iloc[0, 1] - train.iloc[ROWS_PER_SEGMENT, 1]\nMIN_FREQ = round(1/TIME_PER_SEGMENT)\nMAX_FREQ = 1e10\nTHRESHOLD = 0.1\n\ndef LSP_freq_filtered(df, signal_col='signal', time_col='ttf',\n             nrows=ROWS_PER_SEGMENT, min_freq=MIN_FREQ, max_freq=MAX_FREQ):\n    print('Lomb-Scargle Periodogram analysis commencing.')\n    print('Minimum detection frequency: {}Hz'.format(MIN_FREQ))\n    print('Manual maximum frequency cutoff: {}Hz'.format(MAX_FREQ))\n    #initialise empty arrays for frequency outputs to be concatenated to DataFrame \n    freq_1 = np.zeros(len(df))\n    amps_1 = np.zeros(len(df))\n    segment_num = np.zeros(len(df))\n    #loop through input DataFrame in chunks of length=nrows\n    init_id = 0\n    segment_id =1\n    while init_id < len(df):\n        if segment_id==1:\n            print('Processing segment {:d}...'.format(segment_id))\n        if segment_id%500==0:\n            print('Processing segment {:d}...'.format(segment_id))\n        end_id = min(init_id + nrows, len(df))\n        ids = range(init_id, end_id)\n        df_chunk = df.iloc[ids]\n        #np arrays of amplitude and time columns\n        signal = df_chunk[signal_col].values\n        ttf = df_chunk[time_col].values\n        #clear memory\n        del df_chunk\n        gc.collect()\n        #calulate Lomb-Scargle periodograms for spectral analysis\n        freq, amp = LombScargle(ttf, signal).autopower(nyquist_factor=2)\n        freq_df = pd.DataFrame({'freq': freq.round(),\n                               'amp': amp})\n        freq_df = freq_df.loc[freq_df['amp'] >= THRESHOLD] \n        #obtain frequencies sorted by highest amplitude as np.array\n        top_freqs = freq_df.loc[(freq_df.freq > min_freq) & (freq_df.freq < max_freq)].sort_values('amp', ascending=False).freq.values\n        amps = freq_df.loc[(freq_df.freq > min_freq) & (freq_df.freq < max_freq)].sort_values('amp', ascending=False).amp.values\n        del freq_df, freq, amp, signal, ttf\n        gc.collect()\n        #update main frequency component arrays\n        try:\n            freq_1[ids] = top_freqs[0]\n            amps_1[ids] = amps[0]\n        except:\n            freq_1[ids] = 0\n        segment_num[ids] = segment_id\n        del top_freqs\n        init_id += nrows\n        segment_id += 1\n    print('...Done. Adding main component frequencies as DataFrame columns...')\n    df['Freq_periodic'] = freq_1\n    df['Amps_periodic'] = amps_1\n    df['Segment'] = segment_num\n    print('...Done.')\n    \nLSP_freq_filtered(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5ef0515b17d79430108fa515abecdf6e54ba8269"},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6f9b32ae3da6fd112a40cdf7930f9a1a9657ddf9","_kg_hide-input":true},"cell_type":"code","source":"train['Freq_RM'] = train['Freq_periodic'].rolling(window=200000,center=False).std()\ntrain.plot(kind='scatter',x='ttf',y='Freq_RM', figsize=(12, 6), s=1)\nplt.xlim(1.5, -0.05)\nplt.axvline(0.32)\nplt.text(0.3,2e8,'Main Quake')\nplt.title('Lomb-Scargle Main Component Frequency Standard Deviation (100MHz): Rolling-Mean ')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"7f66f0244d771cb6c790050779afe56978efc941"},"cell_type":"code","source":"train['Amp_RM'] = train['Amps_periodic'].rolling(window=200000,center=False).mean()\ntrain.plot(kind='scatter',x='ttf',y='Amp_RM', figsize=(12, 6), s=1)\nplt.xlim(1.5, -0.05)\nplt.title('Lomb-Scargle Main Component Frequency Amplitude: Rolling-Mean')\nplt.axvline(0.32)\nplt.text(0.3,0.01,'Main Quake')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fd55528770098528033fdd4cc70fc34f1dea2050"},"cell_type":"markdown","source":"These are potentially significant results for predicting a TTF < 0.2 seconds. Following an earthquake event, the statistical randomness of the noise increases, so the maximum amplitude for the main frequency detected by LSP drops below the mean. Meanwhile, the standard deviation of the main frequency isolated by LSP decreases. Physical reasons for this may include a decrease in the acoustic energy, and an increase in entropy following the rapid release of tension in a system. I'll leave that to the researchers to explain! \n\nThere are other features that also decrease as TTF approaches 0, so nothing in this kernel will be a magic bullet. But there may still be some useful features that can be engineered from LSP frequency analysis, and I hope this kernel helps my fellow competitors. Even if these conclusions aren't helpful to you, there is much for you to build on.Good luck! "}],"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}