{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"My First Kaggle Kernel with lot of enthusiasm! Eager to listen from all Kagglers!  All suggestions, improvements, ideas, questions and queries are well come.\nThe signal is sum of sinusoidal wave, other interfering waves, random disturbances and measurement noise. The sinusoidal wave is the indented signal and rest is deviations from various sources. One of the sources is fault on line. Using superimposition of wave theory we can subtract the sinusoidal wave from signal and get all deviation. From the difference in faulty and non-faulty residuals we can get patterns exhibited by fault.\nTo reach this goal, we divide it following sub tasks,\n1.\tSmooth the signal so we can fit sine wave\n2.\tFit the sine wave\n3.\tCalculated residual\n4.\tCompare some good and faulty signals in term of residuals.\n\nThe optimize library from scipy is used for sine wave fitting. The find_peak function from scipy.signal is used to identify peaks and Kmeans from sklearn.cluster is used for clustering those peaks. The smooth function denoises data and sine_func implements the sine function.\n"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"code","source":"import os\nimport pyarrow.parquet as pq\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom scipy import optimize\nfrom scipy.signal import find_peaks\nfrom sklearn.cluster import KMeans\n\n#The number of sample in every record\ndataPoints=800000\ndef smooth(y, box_pts):\n    box = np.ones(box_pts)/box_pts\n    y_smooth = np.convolve(y, box, mode='same')\n    return y_smooth\n\ndef sine_func(x, a, b, c,d):\n    return a * np.sin((b*x)+c)+d","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"107182add3bd9463a2207afb09873bd316cedaba"},"cell_type":"markdown","source":"Let us read a signal and plot it to see how it looks."},{"metadata":{"trusted":true,"_uuid":"a5344b34baf997f7b9b186a453f74aac6b0448ae"},"cell_type":"code","source":"t = pd.read_parquet('../input/train.parquet', engine='pyarrow', columns=['575'])  #8021 8031 8081  \nplt.figure(figsize=(24, 10))\nplt.plot(t[:],label='Raw series')\nplt.legend(loc='best')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c6f5a0f2c90c8f51b55d5a7aa8922dc35d57e566"},"cell_type":"markdown","source":"The fitting of sine wave is a challenging task, especially to noisy data. To smooth data the numpy's convolve method is used. The second parameter decided the convolution window.  The smoothed wave is clearly a sine wave. The sine wave is fitted to smoothed data."},{"metadata":{"trusted":true,"_uuid":"010197e1eef355ac7d86e133f76dcddca2e4b64d"},"cell_type":"code","source":"y_data = t[0:dataPoints].values.reshape(dataPoints)\nx_data = np.linspace(0,1, num=dataPoints).T.reshape(dataPoints)\ny_data_smooth = smooth(y_data, 10000)\nplt.figure(figsize=(24, 10))\nplt.plot( y_data,'r+')\nplt.plot(y_data_smooth,\"b-\",label='Smoothed Data')\nplt.legend(loc='best')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"855211ca037ae1ec9b731fac676f608b82c48bfd"},"cell_type":"markdown","source":"To fit sine wave, one need to provide initial guess values for parameters like amplitude, phase and shift. Other parameters are known from context we need to find out phase. I will discuss the sine equation and other parameters while fitting sine wave. Now we focus on phase. Phase is the position of a point in time on a waveform cycle.(Further reference https://en.wikipedia.org/wiki/Phase_(waves) ) We estimate the phase  by locating the positive peak. To identify peaks we use find_peaks method from scipy.signal library. We specify width to avoid sudden fluctuations and smaller peaks."},{"metadata":{"trusted":true,"_uuid":"8690eb0a9e20b805cc4f447a055c64ccc0f023c9"},"cell_type":"code","source":"peaks, _ = find_peaks(y_data_smooth,threshold =(0.00001,0.0001),width = 100)\nplt.plot(y_data_smooth)\nplt.plot(peaks, y_data_smooth[peaks], \"x\")\nplt.plot(np.zeros_like(y_data_smooth), \"--\", color=\"gray\",label='Peaks')\nplt.legend(loc='best')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ac00553aed9286b4ed9a8db862c2b3ec03796d12"},"cell_type":"code","source":"The more than two peaks are found. We cluster this peak into two clusters and take average \nfor finding their position on x-axis. Kmeans cluster with 2 numbers of clusters returns the cluster \nindex for every peak. We check which is peak is positive by comparing them and take average of \npeak positions for it. From this location of positive peak we calculated the phase with below reasoning.\n\nThe entire signel is capturing a full wave i.e 360 degress or 2 pi in 800000 time slices. \nSo the every displacement counts for 2*pi/800000. \nThe positive peack will appear at pi/2 for zero phase. and hence\nphase = (((800000- posivePeackLoc)+200000)% 800000)/800000*np.pi/2         ---phase formula\n\nThe sine wave has equation,\ny(t) = A\\sin(2 \\pi f t + \\varphi) = A\\sin(\\omega t + \\varphi)\nwhere:\n\nA = the amplitude, the peak deviation of the function from zero.\nf = the ordinary frequency, the number of oscillations (cycles) that occur each second of time.\nω = 2πf, the angular frequency, the rate of change of the function argument in \n        units of radians per second \n{\\displaystyle \\varphi } \\varphi  = the phase, specifies (in radians) where in its cycle \n        the oscillation is at t = 0.    When {\\displaystyle \\varphi } \\varphi  is non-zero, \n    the entire waveform appears to be shifted in time by the amount {\\displaystyle \\varphi } \n    \\varphi /ω seconds. A negative value represents a delay, and a positive value represents an advance.\n\nWe also add a displacement parameter so the fitted curve can move horizontally.\n\nThe intial guasses are below,\nA = 20 visual inspection  and problem description.\nf = 1/800000. There are 800000 samples per full cycle.\nω = 2πf = 4*np.pi/800000\nphase = (((800000- posivePeackLoc)+200000)% 800000)/800000*np.pi/2         ---phase formula \n\nThe fitted wave is plotted in blue over red raw signal. Good fit!","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6fb5dc278b27751ae6df6064189291b4619b50ba"},"cell_type":"code","source":"km = KMeans(n_clusters=2)\nkm.fit(peaks.reshape(-1,1))\nclu = km.predict(peaks.reshape(-1,1))\n\nif ( np.mean(y_data_smooth[peaks[clu==1]]) > np.mean(y_data_smooth[peaks[clu==0]])) :\n        posivePeackLoc = np.mean(peaks[clu==1])\nelse:\n        posivePeackLoc = np.mean(peaks[clu==0])\n\ny_data = t[0:dataPoints].values.reshape(dataPoints)\nx_data = np.linspace(0,1, num=dataPoints).T.reshape(dataPoints)\ny_data[y_data>45]=45\ny_data[y_data<-45]=-45\nphase = (((800000- posivePeackLoc)+200000)% 800000)/800000*np.pi/2\nparams, params_covariance = optimize.curve_fit(sine_func, x_data, y_data,\n                                               p0=[20, 4*np.pi/800000, phase,5])\n\nplt.figure(figsize=(24, 10))\nplt.plot( y_data,'r+')\nplt.plot(x_data*dataPoints, sine_func(x_data, params[0], params[1],params[2],params[3]),\n         label='Fitted Sine Wave')\nplt.legend(loc='best')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a29616ef79bdde2e0461bd4559ffb8de58ba62b0"},"cell_type":"markdown","source":"Now getting residul is very easy. Just substract the fitted sine wave from raw data. Below is residual for non fault case. (signal id 575)"},{"metadata":{"trusted":true,"_uuid":"1312df627bae51a49436ee09cd682542c744604e"},"cell_type":"code","source":"residual=y_data-sine_func(x_data, params[0], params[1],params[2],params[3])\nplt.figure(figsize=(24, 10))\nplt.plot(x_data*dataPoints, y_data-sine_func(x_data, params[0], params[1],params[2],params[3]),\n         label='Residual')\nplt.legend(loc='best')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b711d53dc32732be692a778c2e20cb0e3591c748"},"cell_type":"markdown","source":"[Imgur](https://i.imgur.com/X1PnOr3.png)"},{"metadata":{"_uuid":"1347eea7e273fe4dd20951b86683604d7bb2de1c"},"cell_type":"markdown","source":"These lines are disturbances. The circular patterns are due to fittment error at peak. Please compare this image with the Residual graph for faulty case. (signal id 7919) shown below,"},{"metadata":{"_uuid":"216cfde359a61df4fbe97e1d4ad16086f9bd600c"},"cell_type":"markdown","source":"[Imgur](https://i.imgur.com/ERYuetn.png)"},{"metadata":{"_uuid":"06de72993e4ca9087cb48c41d9b3cc57fff1368e"},"cell_type":"markdown","source":"The two images indicate 1) the frequencies of vertical lines is more in faulty lines.\n                                            2) The width of these lines is more in faulty lines.\n                                            3) The height of these lines is more in faulty lines.\n                                            4) The fitment error is dominated in non-faulty lines.\n   These study suggests that we can use some image processing, CNN, LSTM to use visual clues. Also we can extract has disturbance frequencies, time length, density for feature selection."},{"metadata":{"_uuid":"7abbaed528786be1190d7c02b5802918ef5722a3"},"cell_type":"markdown","source":""}],"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}