{"cells":[{"metadata":{"_uuid":"c1b8b3283f6bf77de1ced2733040943ff7041ab7"},"cell_type":"markdown","source":"# Analyze Power Line Signal Like a Physicist\nWe read and explore the data, especially the labels. Then, we focus on the frequency domain. Digital filtering with infinitely fast roll-off (sharp cutoff) is demonstrated. \nStatistician likes to use moving average and other smoothing techniques, which are basically low-pass filters with very slow roll-off. Engineers prefer more realistic filters with finite roll-off, because they have to implement filter in the real world. That is why scipy.signal provides an array of different filters. If you are a physicist who doesn't care about real world, why don't we just filter with infinitely fast roll-off?  \nBefore we begin, I would like to thank https://www.kaggle.com/xhlulu/exploring-signal-processing-with-scipy and the host: https://www.kaggle.com/sohier/reading-the-data-with-python"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nfrom scipy import fftpack, signal\nfrom matplotlib import pyplot as plt\nimport pandas as pd\nimport seaborn as sns\nimport pyarrow.parquet as pq","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"242ba94e7ea14877e8f6df052ef74800d0fa691f"},"cell_type":"markdown","source":"## Read and Explore\nRead signals (variable x). "},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"scrolled":true},"cell_type":"code","source":"xs = pq.read_table('../input/train.parquet', columns=[str(i) for i in range(999)]).to_pandas()\n# xs = pq.read_table('../input/train.parquet').to_pandas()\nprint((xs.shape))\nxs.head(2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b3627b19a6d8416e15e2fa7ec5af00071c80990c"},"cell_type":"raw","source":"Read labels (variable y). "},{"metadata":{"trusted":true,"_uuid":"17892210ab795f1ca0d2fb0f516f849f1a094ef1"},"cell_type":"code","source":"train_meta = pd.read_csv('../input/metadata_train.csv')\nprint(train_meta.shape)\ntrain_meta.head(6)\ntrain_meta_good = train_meta[train_meta.target == 0]\ntrain_meta_error = train_meta[train_meta.target == 1]\nid_error = train_meta_error.groupby('id_measurement')['target'].count()\nid_error.head(6)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ae3e758795b955bf50488de41adc73598824ddf9"},"cell_type":"markdown","source":"Examine how many phases with one ID were labelled problematic. For 80.4% faulty lines, three phases were all labelled faulty, while one-faulty-phase and two-faulty-phase lines contribute about 10% and 10%, respectively.  "},{"metadata":{"trusted":true,"_uuid":"ba2bbe3ca5eec3d9a4530578d36157bdd81e7e89"},"cell_type":"code","source":"id_error_c = id_error.astype('category')\nprint(id_error_c.value_counts())\nprint(id_error_c.value_counts() / id_error_c.value_counts().sum())\nid_error_c.value_counts().plot(kind='bar')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"16a4bdba3ce8895218df5fb5d2d97ddb1c2774ab"},"cell_type":"code","source":"period = 0.02\ntime_step = 0.02 / 800000.\ntime_vec = np.arange(0, 0.02, time_step)\nf_sampling = 1 / time_step\nprint(f'Sampling Frequency = {f_sampling / 1e6} MHz')\n# print (str(50* 800000 /1e6) + ' MHz')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c3ba4c5160e99c6bebc83ead6f4dc2026451e1d1"},"cell_type":"code","source":"# Fetch one signal from xs\nidx = 1\nsig = xs.iloc[:, idx]\nidx_error = 3\nsig_error = xs.iloc[:, idx_error]\nprint(sig.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a1d699380d0ca33f26176c38525c4185a013a0ed"},"cell_type":"code","source":"# https://www.scipy-lectures.org/intro/scipy/auto_examples/plot_fftpack.html\n# The FFT of the signal\nsig_fft = fftpack.fft(sig)\n# And the power (sig_fft is of complex dtype)\npower = np.abs(sig_fft)\n# The corresponding frequencies\nsample_freq = fftpack.fftfreq(sig.size, d=time_step)\n\n# Find the peak frequency: we can focus on only the positive frequencies\npos_mask = np.where(sample_freq >= 0)\nfreqs = sample_freq[pos_mask]\npeak_freq = freqs[power[pos_mask].argmax()]\n\nplt.figure(figsize=(6, 5))\n# plt.plot(sample_freq[pos_mask], power[pos_mask])\nplt.semilogy(sample_freq[pos_mask], power[pos_mask])\nplt.ylim([1e-0, 1e8])\nplt.xlabel('Frequency [Hz]')\nplt.ylabel('Power [A.U./Hz]')\n\n\n# Check that it does indeed correspond to the frequency that we generate\n# the signal with\nnp.allclose(peak_freq, 1./period)\n\n# An inner plot to show the peak frequency\naxes = plt.axes([0.55, 0.6, 0.3, 0.2])\nplt.title('Peak frequency')\nplt.plot(freqs[:8], power[:8])\nplt.setp(axes, yticks=[])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"699dca70d7baacbaa5b3098733d8b1b8ca0b49ce"},"cell_type":"markdown","source":"The inset figure above shows that peak frequency is 50 Hz, as we expected. Note that the frequency step in fft spectrum is 50 Hz, limited by the total duration of the signal. \nNext we use a slightly different method (signal.periodogram) and plot the power spectrum in log-log scale, which make more sense. The unit in the y axis is different by a fixed factor. It is OK. "},{"metadata":{"trusted":true,"_uuid":"7451f176b95f81d96fb712f905602eac79ee69b5"},"cell_type":"code","source":"def plot_ps(sig, f_sampling, label='sig', style='loglog'):\n    f, Pxx_den = signal.periodogram(sig, f_sampling)\n    if style == 'semilogy':\n        plt.semilogy(f, Pxx_den, label=label)\n    else:\n        plt.loglog(f, Pxx_den, label=label)\n    plt.ylim([1e-9, 1e2])\n    plt.xlabel('frequency [Hz]')\n    plt.ylabel('PSD [A.U./Hz]')\n    \nplot_ps(sig, f_sampling, 'Good sig')\nplot_ps(sig_error, f_sampling, 'Bad sig')\nplt.legend(loc='best')\nplt.show()\n# The horizontal line at 50 Hz is artificial, as the log scale in the x axis cannot show 0 Hz. ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7a9ce0e714a5d1f394d899b7a886150e5182bbb5"},"cell_type":"code","source":"def bandpassfilter(spec, sample_freq, lowcut, highcut):\n    # a digital bandpass filter with a infinite roll off. \n    # Note that we will keep the frequency point right at low cut-off and high cut-off frequencies. \n    spec1 = spec.copy()\n    spec1[np.abs(sample_freq) < lowcut] = 0\n    spec1[np.abs(sample_freq) > highcut] = 0\n    filtered_sig = fftpack.ifft(spec1)\n    return filtered_sig","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"65e732396d126f4ab3592e07e8a1188d79439bfa"},"cell_type":"markdown","source":"## Digital filtering\nThe peak_freq should be 50 Hz. \nWe demonstrated differnt low-pass, high-pass, and band-pass filtered signals. You can see 10-1000 Hz can capture a lot of low frequency features. \n"},{"metadata":{"trusted":true,"_uuid":"fdca023ca1fb5a4b772d5ac23d14c457df014577"},"cell_type":"code","source":"# We demonstrated differnt low-pass filtered signals. You can see 10-1000 Hz can capture a lot of low frequency features. \nlowcut, highcut = 10, 100\nfiltered_sig0 = bandpassfilter(sig_fft,sample_freq, lowcut, highcut)\nlowcut, highcut = 10, 300\nfiltered_sig1 = bandpassfilter(sig_fft,sample_freq, lowcut, highcut)\nlowcut, highcut = 10, 1000\nfiltered_sig2 = bandpassfilter(sig_fft,sample_freq, lowcut, highcut)\n\nplt.figure(figsize=(6, 5))\nplt.plot(time_vec, sig, label='Original signal')\nplt.plot(time_vec, filtered_sig0, linewidth=3, label='10-100 Hz')\nplt.plot(time_vec, filtered_sig1, linewidth=3, label='10-300 Hz')\nplt.plot(time_vec, filtered_sig2, linewidth=3, label='10-1000 Hz')\nplt.xlabel('Time [s]')\nplt.ylabel('Amplitude')\nplt.legend(loc='best')\n\n# We also demonstrate a band-pass filtered and a high-pass filtered signals. \nlowcut, highcut = 1000, 1e6\nfiltered_sig3 = bandpassfilter(sig_fft,sample_freq, lowcut, highcut)\nlowcut, highcut = 1000, 40e6\nfiltered_sig4 = bandpassfilter(sig_fft,sample_freq, lowcut, highcut)\nplt.figure(figsize=(6, 5))\nplt.plot(time_vec, sig, label='Original signal')\nplt.plot(time_vec, filtered_sig4, linewidth=3, label='Above 1 kHz')\nplt.plot(time_vec, filtered_sig3, linewidth=3, label='1 kHz-1 MHz')\nplt.xlabel('Time [s]')\nplt.ylabel('Amplitude')\nplt.legend(loc='best')\n","execution_count":null,"outputs":[]}],"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}