{"cells":[{"metadata":{"_uuid":"abda47372d10edee85b20464ec202ab41ac3c462"},"cell_type":"markdown","source":"## Short kernel for aligning waves based on their phases, which are obtained from fft.  "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nfrom scipy import fftpack # Fast Fourier Transform functions\n\nimport pyarrow.parquet as pq # for reading input data\n\n%matplotlib inline\nimport matplotlib.pyplot as plt # plotting\n\nplt.style.use('seaborn-whitegrid')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"**Load data**"},{"metadata":{"trusted":true,"_uuid":"7c1ebe67038973fc3c7a1ede8e07d1d5c98fb60b"},"cell_type":"code","source":"# load just three first signals (which belong to the same id measurement)\nn_signals_to_load = 3\nsignals = pq.read_pandas(\n    '../input/train.parquet', \n    columns=[str(i) for i in range(n_signals_to_load)]).to_pandas()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"521187e92ae6da2f1852e6ab870de5df916da332"},"cell_type":"markdown","source":"**Bookkeeping**"},{"metadata":{"trusted":true,"_uuid":"3a0f60818c2362258ee2d71d1b389a935395777a"},"cell_type":"code","source":"# sampling rate\nnum_samples = signals.shape[0] # 800,000 samples per signal\nperiod = 0.02 # over a 20ms period\nfs = num_samples / period # 40MHz sampling rate\n\n# time array support\nt = np.array([i / fs for i in range(num_samples)])\n\n# frequency vector fro FFT\nfreqs = fftpack.fftfreq(num_samples, d=1/fs)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b554080a49429f48a62a5e1bd9f84d20d46af9a"},"cell_type":"markdown","source":"# Outline\n\nThere are four main steps for each signal:\n1. Get FFT coefficients\n2. Find the coefficients with highest norm (should correspond to 50Hz)\n3. Find the phase of the complex coefficient\n4. Get the instant angular phase vector $\\omega_i$ of the main signal ($f_0 = 50$Hz), i.e. $\\omega_i = 2 \\pi t_i  f_0 + \\phi$, where $\\phi$ is found in step 3\n\nAfterwards, you only need to arbitrarily define a phase to align each signal (e.g. $\\frac{\\pi}{2}$)."},{"metadata":{"trusted":true,"_uuid":"b95bc56f3bd3f4450a31ec38fdc30e92ff542eca"},"cell_type":"code","source":"# get fft coeffs\ndef get_fft_coeffs(sig):\n    return fftpack.fft(sig)\n\n# get coeff with highest norm\ndef get_highest_coeff(fft_coeffs, freqs, verbose=True):\n    coeff_norms = np.abs(fft_coeffs) # get norms (fft coeffs are complex)\n    max_idx = np.argmax(coeff_norms)\n    max_coeff = fft_coeffs[max_idx] # get max coeff\n    max_freq = freqs[max_idx] # assess which is the dominant frequency\n    max_amp = (coeff_norms[max_idx] / num_samples) * 2 # times 2 because there are mirrored freqs\n    if verbose:\n        print('Dominant frequency is {:,.1f}Hz with amplitude of {:,.1f}\\n'.format(max_freq, max_amp))\n    \n    return max_coeff, max_amp, max_freq\n\n# get max coeff phase\ndef get_max_coeff_phase(max_coeff):\n    return np.angle(max_coeff)\n\n# construct the instant angular phase vector indexed by pi, i.e. ranges from 0 to 2\ndef get_instant_w(time_vector, f0, phase_shift):\n    w_vector = 2 * np.pi * time_vector * f0 + phase_shift\n    w_vector_norm = np.mod(w_vector / (2 * np.pi), 1) * 2 # range between cycle of 0-2 \n    return w_vector, w_vector_norm\n\n# find index of chosen phase to align\ndef get_align_idx(w_vector_norm, align_value=0.5):\n    candidates = np.where(np.isclose(w_vector_norm, align_value))\n    # since we are in discrete time, threre could be many values close to the desired one\n    # so let's take the one in the middle\n    return int(np.median(candidates))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"90185a604beb7bf8265fc05a2bcca2a4cd92a511"},"cell_type":"code","source":"fig = plt.figure(figsize=(16, 9))\nplot_number = 0\n\nfor signal_id in signals.columns:\n    # get samples\n    print('=== Signal {} ==='.format(signal_id))\n    sig = signals[signal_id]\n    \n    # fft\n    fft_coeffs = get_fft_coeffs(sig)\n    \n    # asses dominant frequency\n    max_coeff, amp, f0 = get_highest_coeff(fft_coeffs, freqs)\n    \n    # phase shift\n    ps = get_max_coeff_phase(max_coeff)\n    \n    # get angular phase vector\n    w, w_norm = get_instant_w(t, f0, ps)\n    \n    # generate dominant signal at f0\n    dominant_wave = amp * np.cos(w) # if np.sin(), then need to ajust by pi/2\n    (w)\n    \n    # plot signals\n    plot_number += 1\n    ax = fig.add_subplot(3, 2, plot_number)\n    \n    ax.plot(t * 1000, sig, label='Original') # original signal\n    ax.plot(t * 1000, dominant_wave, color='red', label='Wave at {:.0f}Hz'.format(f0)) # wave at f0\n    ax.legend()\n    ax.set_xlabel('time (ms)')\n    ax.set_ylabel('Amplitude')\n    ax.set_title('Signal {}'.format(signal_id))\n    \n    # plot phase\n    plot_number += 1\n    ax = fig.add_subplot(3, 2, plot_number)\n    \n    ax.plot(t * 1000, w_norm, label='phase') # instant phase\n    ax2 = ax.twinx() # secondary y\n    ax2.plot(t * 1000, dominant_wave, color='red', label='Wave at {:.0f}Hz'.format(f0)) # wave at f0\n    ax.legend()\n    ax.set_xlabel('time (ms)')\n    ax.set_ylabel('$\\omega_i (\\pi$ rad)')\n    ax2.set_ylabel('Wave amplitude')\n    ax.set_title('Instant angular phase of dominant wave')\n    \nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ff35d4f2d3949b70748995e4a1e1cb9fd5ce0a4e"},"cell_type":"markdown","source":"All amplitude should be the same, since it is the same triphasic circuit, but FFT cannot recover it perfectly due to noise."},{"metadata":{"_uuid":"cbde8e4ff3d57a7858a23af1a21caac86da313ba"},"cell_type":"markdown","source":"**Align waves**"},{"metadata":{"trusted":true,"_uuid":"428cf2ffb2235b21ad9ec822a890cb6b362f088a"},"cell_type":"code","source":"# align waves with np.roll()\nalign_phase = 0.5 # w_i = pi/2\n\n\nfig = plt.figure(figsize=(12, 9))\nplot_number = 0\n\nfor signal_id in signals.columns:\n    # get samples\n    sig = signals[signal_id]\n    \n    # fft\n    fft_coeffs = get_fft_coeffs(sig)\n    \n    # asses dominant frequency\n    max_coeff, amp, f0 = get_highest_coeff(fft_coeffs, freqs, verbose=False)\n    \n    # phase shift\n    ps = get_max_coeff_phase(max_coeff)\n    \n    # get angular phase vector\n    w, w_norm = get_instant_w(t, f0, ps)\n    \n    # generate dominant signal at f0\n    dominant_wave = amp * np.cos(w)\n    \n    # idx to roll\n    origin = get_align_idx(w_norm, align_value=align_phase)\n    \n    # roll signal and dominant wave\n    sig_rolled = np.roll(sig, num_samples - origin)\n    dominant_wave_rolled = np.roll(dominant_wave, num_samples - origin)\n    \n    # plot signals\n    plot_number += 1\n    ax = fig.add_subplot(3, 1, plot_number)\n    \n    ax.plot(t * 1000, sig_rolled, label='Rolled Original') # original signal\n    ax.plot(t * 1000, dominant_wave_rolled, color='red', label='Rolled Wave at {:.0f}Hz'.format(f0)) # wave at f0\n    ax.legend()\n    ax.set_xlabel('time (ms)')\n    ax.set_ylabel('Amplitude')\n    ax.set_title('Signal {} rolled'.format(signal_id))\n    \nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68eddafb521c7ce3859edf48d7392826527f00f3"},"cell_type":"markdown","source":"Signals aligned. It should work with any desired `align_phase`. Let's try with $\\pi / 4$"},{"metadata":{"trusted":true,"_uuid":"414b5d1fcab95db28514cd357da554a54086e0b2"},"cell_type":"code","source":"# align waves with np.roll()\nalign_phase = 0.25 # w_i = pi/4\n\n\nfig = plt.figure(figsize=(12, 9))\nplot_number = 0\n\nfor signal_id in signals.columns:\n    # get samples\n    sig = signals[signal_id]\n    \n    # fft\n    fft_coeffs = get_fft_coeffs(sig)\n    \n    # asses dominant frequency\n    max_coeff, amp, f0 = get_highest_coeff(fft_coeffs, freqs, verbose=False)\n    \n    # phase shift\n    ps = get_max_coeff_phase(max_coeff)\n    \n    # get angular phase vector\n    w, w_norm = get_instant_w(t, f0, ps)\n    \n    # generate dominant signal at f0\n    dominant_wave = amp * np.cos(w)\n    \n    # idx to roll\n    origin = get_align_idx(w_norm, align_value=align_phase)\n    \n    # roll signal and dominant wave\n    sig_rolled = np.roll(sig, num_samples - origin)\n    dominant_wave_rolled = np.roll(dominant_wave, num_samples - origin)\n    \n    # plot signals\n    plot_number += 1\n    ax = fig.add_subplot(3, 1, plot_number)\n    \n    ax.plot(t * 1000, sig_rolled, label='Rolled Original') # original signal\n    ax.plot(t * 1000, dominant_wave_rolled, color='red', label='Rolled Wave at {:.0f}Hz'.format(f0)) # wave at f0\n    ax.legend()\n    ax.set_xlabel('time (ms)')\n    ax.set_ylabel('Amplitude')\n    ax.set_title('Signal {} rolled'.format(signal_id))\n    \nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1edf9fea748926ece99e078219990d28f87ba3fd"},"cell_type":"markdown","source":"Worked!"},{"metadata":{"trusted":true,"_uuid":"62c6707f962f86853167271ffcb07e3ba7e0dd31"},"cell_type":"code","source":"","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}