{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Мустаев Айрат (@airat)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Обработка данных с АЦП ZetLab "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В рамках исследования по классификации отдельных типов объектов (человек, группа людей, животное) фиксируемых сигнализационными датчиками, таким как например: <b>точечный сейсмомагнитометрический датчик</b>, проводился эксперемент. В данной статье рассматривается алгоритм обработки эксперементальных (сырых) данных, получаемых с АЦП [(АЦП ЦАП Zet 230)](https://zetlab.com/shop/izmeritelnoe-oborudovanie/moduli-atsp-tsap/atsp-tsap-zet-230), для последуещего формирования признакового пространства и классификации. \n",
    "\n",
    "Под сырыми данными, понимаются синхронизированные по времени две реализации (сейсмический и магнитометрический сигналы) с частотой дискретизации Fs = 500 Гц. \n",
    "\n",
    "### Краткое описание природы данных\n",
    "\n",
    "#### Сейсмическая волна\n",
    "Сейсмическая волна - это волны, переносящие энергию упругих (механических) колебаний в горных породах. Источником сейсмической волны может быть землетрясение, взрыв, вибрация или <b>удар (в нашем случае проход объекта классификации)</b>. \n",
    "<p>Существует следующая классификация сейсмический волн:</p>\n",
    " * Объёмные волны - Объёмные волны проходят через недра Земли. Путь волн преломляется различной плотностью и жёсткостью подземных пород.\n",
    " * P-волны (первичные волны) — продольные, или компрессионные волны. Обычно их скорость в два раза быстрее S-волн, проходить они могут через любые материалы.\n",
    " * P- и S-волны в мантии и ядре.\n",
    " * Поверхностные волны несколько похожи на волны воды, но в отличие от них они путешествуют по земной поверхности. Их обычная скорость значительно ниже скорости волн тела. Из-за своей низкой частоты, времени действия и большой амплитуды они являются самыми разрушительными изо всех типов сейсмических волн.<img src=\"Image/Rayliegh_Love.jpg\" width=\"500\" height=\"500\" >\n",
    "\n",
    " \n",
    "#### Магнитометрический сигнал\n",
    "Магнитометрический сигнал - это наведеное ЭДС в катушке при переносе над плоскостью катушки ферромагнитоного предмета."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загрузим файлы сигналов:\n",
    " * для магнитометрических сигналов файл magnito.txt;\n",
    " * для сейсмических сигналов файл seismic.txt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_Magnito = open('file_signal/magnito.txt', 'r')\n",
    "data_m = file_Magnito.read().split(\"\\n\")\n",
    "file_Seismic = open('file_signal/seismic.txt', 'r')\n",
    "data_s = file_Seismic.read().split(\"\\n\")\n",
    "data_s = list(map(lambda signal: \n",
    "             list(map(lambda value: \n",
    "                 float(value), \n",
    "                     filter(lambda value2: value2 != '' and len(value2) > 0, \n",
    "                            signal.split(';')\n",
    "                           )\n",
    "                )), \n",
    "             data_s))\n",
    "\n",
    "data_m = list(map(lambda signal: \n",
    "             list(map(lambda value: \n",
    "                 float(value), \n",
    "                     filter(lambda value2: value2 != '' and len(value2) > 0, \n",
    "                            signal.split(';')\n",
    "                           )\n",
    "                )), \n",
    "             data_m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Отобразим реализацию сейсмического и магнитометрического сигнала"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для примера возьмем реализацию сигналов с индексом 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fs=500 # Частота дискретизации\n",
    "num_signal = 144\n",
    "magnito_signal = data_m[num_signal][144:5004]\n",
    "seismic_signal = data_s[num_signal][144:5004]\n",
    "time_x = np.linspace(0,len(seismic_signal)/Fs, len(seismic_signal))\n",
    "plt.figure(figsize = (15, 6))\n",
    "plt.title('Совмещенные по времени реализации сейсмического и магнитометрического сигналов')\n",
    "plt.plot(time_x,magnito_signal, 'r')\n",
    "plt.plot(time_x,seismic_signal, )\n",
    "plt.legend(['Магнитометрический сигнал', 'Сейсмический сигнал'])\n",
    "plt.xlabel('Время, с')\n",
    "plt.ylabel('Амплитуда')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Полученные сигналы сильно зашумлены, а для магнитометрического сигнала и вовсе не удается его рассмотреть."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Рассмотрим спектры сигналов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для построения спекторов воспользуемся библиотекой scipy и функцией fft для быстрого преобразования Фурье"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "from scipy import fft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Подготовим функцию построения, визуализации спектра сигнала и самого сигнала"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CreatSpectrum(signal,Fs):\n",
    "    n = len(signal) # Длительность сигнала\n",
    "    k = arange(n)\n",
    "    T = n/Fs\n",
    "    frq = k/T \n",
    "    frq = frq[range(int(n/2))] # Диапазон частот\n",
    "    \n",
    "    Y = fft(signal)/n # Вычисление быстрого преобразования Фурье и его нормализация\n",
    "    Y = np.abs(Y[range(int(n/2))])\n",
    "    #print(np.max(Y))\n",
    "    #print (list(Y).index(np.max(Y)))\n",
    "    #print (frq[list(Y).index(np.max(Y))])\n",
    "    PlotSpectrum(frq, Y)\n",
    "    \n",
    "def PlotSpectrum(frq, Y):\n",
    "    \n",
    "    plot(frq,abs(Y),'r')\n",
    "    title('Спектр сигнала')\n",
    "    xlabel('Частота (Гц)')\n",
    "    ylabel('|Y(freq)|')\n",
    "    plt.grid()\n",
    "    \n",
    "def PlotSignal(time_x,signal, name):\n",
    "    \n",
    "    plot(time_x,signal)\n",
    "    title(name+' сигнал')\n",
    "    xlabel('Время, с')\n",
    "    ylabel('Амплитуда')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Отобразим магнитометрический сигнал и его спектр, и тоже самое сделаем с сейсмическим сигналом"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (15, 8))\n",
    "subplot(2,2,1)\n",
    "PlotSignal(time_x,magnito_signal, 'Магнитометрический')\n",
    "subplot(2,2,2)\n",
    "CreatSpectrum(magnito_signal,Fs)\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (15, 8))\n",
    "subplot(2,2,1)\n",
    "PlotSignal(time_x,seismic_signal, 'Сейсмический')\n",
    "subplot(2,2,2)\n",
    "CreatSpectrum(seismic_signal,Fs)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "При рассмотрении спектров реализаций сильно выделяется составляющая равна 50 Гц и ее гармоники (100, 150, 200, 250 Гц), 50 Гц помеховая частота по питанию, которую требуется подавить."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Фильтрация сигналов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмортим четыре типа электронных фильтров:\n",
    " * Фильтр Баттерворта - он проектируется так, чтобы его амплитудно-частотная характеристика была <b>максимально гладкой</b> на частотах полосы пропускания.\n",
    " * Фильтр Чебышёва - он проектируется так, чтобы амплитудно-частотная характеристика имела <b>более крутой спад </b> на частотах полосы пропускания.\n",
    " * Эллиптический фильтр (Фильтр Кауэра) — характерной особенностью которого являются <b>пульсации</b> амплитудно-частотной характеристики как в полосе пропускания, так и полосе подавления. \n",
    " \n",
    "Построим АЧХ данных фильтров:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * Фильтр Баттерворта"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def butter_plot(cutOff, order = 5):\n",
    "    leg =[]\n",
    "    b, a = signal.butter(order, cutOff, 'low', analog=True)\n",
    "    w, h = signal.freqs(b, a)\n",
    "    plt.semilogx(w, 20 * np.log10(abs(h)), linewidth=3)\n",
    "    plt.title('Фильтр Баттерворта')\n",
    "    plt.xlabel('Частота, Гц')\n",
    "    plt.ylabel('Амплитуда, дБ')\n",
    "    leg.append('фильтр Баттерворта')\n",
    "    plt.axvline(cutOff, color='red')\n",
    "    leg.append('частота среза = ' + str(cutOff)+' Гц')\n",
    "    plt.legend(leg)\n",
    "    plt.grid(which='both', axis='both')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * Фильтры Чебышёва 1 и 2 рода"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cheby1_plot(cutOff, rp, order = 5):\n",
    "    leg =[]\n",
    "    for rp_i in rp:\n",
    "        b, a = signal.cheby1(order, rp_i, cutOff, 'low', analog=True)\n",
    "        w, h = signal.freqs(b, a)\n",
    "        plt.semilogx(w, 20 * np.log10(abs(h)), linewidth=3)\n",
    "        plt.title('Фильтра Чебышёва 1 рода')\n",
    "        plt.xlabel('Частота, Гц')\n",
    "        plt.ylabel('Амплитуда, дБ')\n",
    "        leg.append('Пульсация частоты = '+str(rp_i)+' дБ')\n",
    "    plt.axvline(cutOff, color='red')\n",
    "    leg.append('частота среза = ' + str(cutOff)+' Гц')\n",
    "    plt.legend(leg)\n",
    "    plt.grid(which='both', axis='both')\n",
    "\n",
    "def cheby2_plot(cutOff, rp, order = 5):\n",
    "    leg =[]\n",
    "    for rp_i in rp:\n",
    "        b, a = signal.cheby2(order, rp_i, cutOff, 'low', analog=True)\n",
    "        w, h = signal.freqs(b, a)\n",
    "        plt.semilogx(w, 20 * np.log10(abs(h)), linewidth=3)\n",
    "        plt.title('Фильтр Чебышёва 2 рода')\n",
    "        plt.xlabel('Частота, Гц')\n",
    "        plt.ylabel('Амплитуда, дБ')\n",
    "        leg.append('Пульсация частоты = '+str(rp_i)+' дБ')\n",
    "    plt.axvline(cutOff, color='red')\n",
    "    leg.append('частота среза = ' + str(cutOff)+' Гц')\n",
    "    plt.legend(leg)\n",
    "    plt.grid(which='both', axis='both')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ellip_plot(cutOff, rp, rc, order = 5):\n",
    "    leg =[]\n",
    "    for rp_i in rp:\n",
    "        for rc_i in rc:\n",
    "            b, a = signal.ellip(order, rp_i, rc_i, cutOff, 'low', analog=True)\n",
    "            w, h = signal.freqs(b, a)\n",
    "            plt.semilogx(w, 20 * np.log10(abs(h)), linewidth=3)\n",
    "            plt.title('Фильтр Эллиптический')\n",
    "            plt.xlabel('Частота, Гц')\n",
    "            plt.ylabel('Амплитуда, дБ')\n",
    "            leg.append('Пульсация = '+str(rp_i)+' дБ' + ', затухание = '+str(rc_i)+' дБ')\n",
    "    plt.axvline(cutOff, color='red')\n",
    "    leg.append('частота среза = ' + str(cutOff)+' Гц')\n",
    "    plt.legend(leg)\n",
    "    plt.grid(which='both', axis='both')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Визуализируем АЧХ подготовленных фильтров"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fCut = 100 # Частота среза\n",
    "plt.figure(figsize = (15, 12))\n",
    "plt.subplot(2,2,1)\n",
    "cheby1_plot(fCut, [5, 25, 45], order = 10)\n",
    "plt.subplot(2,2,2)\n",
    "cheby2_plot(fCut, [5, 20, 40], order = 4)\n",
    "plt.subplot(2,2,3)\n",
    "butter_plot(fCut, order = 10)\n",
    "plt.subplot(2,2,4)\n",
    "ellip_plot(fCut, [5], [20, 50], order = 10)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для примера рассмотрим четыре типа фильтров Баттерворта:\n",
    " * Фильтр нижних частот - пропускает частоты ниже частоты среза и не пропускает частоты выше частоты среза.\n",
    " * Фильтр верхних частот - пропускает частоты выше частоты среза и не пропускает частоты ниже частоты среза.\n",
    " * Полосовый фильтр - получает на вход две частоты нижняя и верхняя частоты среза, который в пределах этих значений пропускает частоты.\n",
    " * Режекторный фильтр - получает на вход две частоты нижняя и верхняя частоты среза, который в пределах этих значений не пропускает частоты."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для фильтрации сигналов используется функция butter, которая позволяет создавать различные по типу и полосам пропускания фильтры, так например:\n",
    " * btype='lowpass' - означает ФНЧ.\n",
    " * btype='highpass' - означает ФВЧ.\n",
    " * btype='bandpass' - означает полосовый фильтр\n",
    " * btype='bandstop' - означает режекторный фильтр\n",
    " \n",
    "Сгенерируем функции для этих фильтров, и визуалиазируем их (будем использовать цифровой фильтр analog = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Фильтр нижних частот\n",
    "def butter_lowpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = cutoff / nyq\n",
    "    b, a = signal.butter(order, normal_cutoff, btype='low', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_lowpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_lowpass(cutoff, fs, order=order)\n",
    "    y = signal.lfilter(b, a, data)\n",
    "    return y\n",
    "def butter_plot_lowpass(cutOff, fs, order=5):\n",
    "    b, a = butter_lowpass(cutOff, fs, order = order)\n",
    "    w, h = signal.freqz(b, a)\n",
    "    #plt.plot(0.5*fs*w/np.pi, 20 * np.log10(abs(h)), 'b', linewidth=5)\n",
    "    plt.plot(0.5*fs*w/np.pi, abs(h), 'b', linewidth=5)\n",
    "    #plt.xlim(0, 0.5*fs)\n",
    "    plt.title(\"АЧХ ФНЧ\")\n",
    "    plt.fill_between(0.5*fs*w/np.pi, 0, np.abs(h), 0.5*fs*w/np.pi< cutOff, alpha = 0.07, color = 'r')\n",
    "    plt.xlabel('Частота, Гц')\n",
    "    plt.axvline(cutOff, color='r', linewidth=2)\n",
    "    plt.legend(['ФНЧ','Частота среза='+str(cutOff)+' Гц'])\n",
    "    plt.grid()\n",
    "butter_plot_lowpass(48, Fs, order = 4)\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Фильтр верхних частот\n",
    "def butter_highpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = cutoff / nyq\n",
    "    b, a = signal.butter(order, normal_cutoff, btype='highpass', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_highpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_highpass(cutoff, fs, order=order)\n",
    "    y = signal.lfilter(b, a, data)\n",
    "    return y\n",
    "def butter_plot_highpass(cutOff, fs, order=5):\n",
    "    b, a = butter_highpass(cutOff, fs, order = order)\n",
    "    w, h = signal.freqz(b, a, worN=8000)\n",
    "    plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b', linewidth=5)\n",
    "\n",
    "    #plt.xlim(0, 0.5*fs)\n",
    "    plt.title(\"АЧХ ФВЧ\")\n",
    "    plt.xlabel('Частота, Гц')\n",
    "    plt.axvline(cutOff, color='r', linewidth=2)\n",
    "    plt.fill_between(0.5*fs*w/np.pi, 0, np.abs(h), 0.5*fs*w/np.pi> cutOff, alpha = 0.07, color = 'r')\n",
    "    plt.legend(['ФВЧ','Частота среза='+str(cutOff)+' Гц'])\n",
    "    plt.grid()\n",
    "butter_plot_highpass(48, Fs)\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Полосовой фильтр\n",
    "def butter_bandpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = [cutoff[0]/nyq, cutoff[1]/nyq]\n",
    "    b, a = signal.butter(order, normal_cutoff, btype='bandpass', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_bandpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_bandpass(cutoff, fs, order=order)\n",
    "    y = signal.lfilter(b, a, data)\n",
    "    return y\n",
    "def butter_plot_bandpass(cutOff, fs, order=5):\n",
    "    b, a = butter_bandpass(cutOff, fs, order = order)\n",
    "    w, h = signal.freqz(b, a, worN=8000)\n",
    "    plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b', linewidth=5)\n",
    "\n",
    "    #plt.xlim(0, 0.5*fs)\n",
    "    plt.title(\"АЧХ полосового фильтра\")\n",
    "    plt.xlabel('Частота, Гц')\n",
    "    plt.axvline(cutOff[0], color='r', linewidth=2)\n",
    "    plt.axvline(cutOff[1], color='r', linewidth=2)\n",
    "    plt.legend(['Полосовой', 'Частота среза='+str(cutOff[0])+' Гц', 'Частота среза='+str(cutOff[1])+' Гц'])\n",
    "    plt.fill_between(0.5*fs*w/np.pi, 0, np.abs(h), 0.5*fs*w/np.pi> cutOff[0], alpha = 0.07, color = 'r')\n",
    "    plt.fill_between(0.5*fs*w/np.pi, 0, np.abs(h), 0.5*fs*w/np.pi< cutOff[1], alpha = 0.07, color = 'r')\n",
    "    plt.grid()\n",
    "butter_plot_bandpass([40, 60], Fs)\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Режекторный фильтр\n",
    "def butter_bandstop(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = [cutoff[0]/nyq, cutoff[1]/nyq]\n",
    "    b, a = signal.butter(order, normal_cutoff, btype='bandstop', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_bandstop_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_bandstop(cutoff, fs, order=order)\n",
    "    y = signal.lfilter(b, a, data)\n",
    "    return y\n",
    "def butter_plot_bandstop(cutOff, fs, order=5):\n",
    "    b, a = butter_bandstop(cutOff, fs, order = order)\n",
    "    w, h = signal.freqz(b, a, worN=8000)\n",
    "    plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b', linewidth=5)\n",
    "    plt.title(\"АЧХ режекторного фильтра\")\n",
    "    plt.xlabel('Частота, Гц')\n",
    "    plt.axvline(cutOff[0], color='r', linewidth=2)\n",
    "    plt.axvline(cutOff[1], color='r', linewidth=2)\n",
    "    plt.legend(['Режекторный', 'Частота среза='+str(cutOff[0])+' Гц', 'Частота среза='+str(cutOff[1])+' Гц'])\n",
    "    plt.fill_between(0.5*fs*w/np.pi, 0, np.abs(h), 0.5*fs*w/np.pi< cutOff[0], alpha = 0.07, color = 'r')\n",
    "    plt.fill_between(0.5*fs*w/np.pi, 0, np.abs(h), 0.5*fs*w/np.pi> cutOff[1], alpha = 0.07, color = 'r')\n",
    "    plt.grid()\n",
    "butter_plot_bandstop([40, 60], Fs)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Применим один раз ФНЧ и четыре раза режекторный фильтр для наших сигналов. \n",
    "\n",
    "Сначала применим фильтры к магнитометрическому сигналу"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sub_plot_signal_spectr(signal, name, Fs):\n",
    "    time_x = linspace(0,len(signal)/Fs, len(signal))\n",
    "    plt.figure(figsize = (15, 8))\n",
    "    subplot(2,2,1)\n",
    "    PlotSignal(time_x,signal, name)\n",
    "    subplot(2,2,2)\n",
    "    CreatSpectrum(signal,Fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_filter = int(Fs / 2 / 50)\n",
    "width_window = 10\n",
    "\n",
    "magnito_signal_filter = magnito_signal\n",
    "magnito_signal_filter = butter_lowpass_filter(magnito_signal_filter, 250-width_window, Fs, order = 5)\n",
    "sub_plot_signal_spectr(magnito_signal_filter, 'Магнитометрический', 500)\n",
    "\n",
    "for i in range(n_filter-1):\n",
    "    magnito_signal_filter = butter_bandstop_filter(magnito_signal_filter, [50*(i+1)-width_window, 50*(i+1)+width_window], Fs)\n",
    "    sub_plot_signal_spectr(magnito_signal_filter, 'Магнитометрический', 500)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Результаты уже лучше, но еще не достаточно хорошо. Вспомним, магнитометрический сигнал расположен в диапозоне от единиц до десятков Гц. Воспользуемся этим и применим ФНЧ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "width_window = 10\n",
    "\n",
    "magnito_signal_filter_hard = magnito_signal_filter\n",
    "magnito_signal_filter_hard = butter_lowpass_filter(magnito_signal_filter_hard, 25, Fs)\n",
    "sub_plot_signal_spectr(magnito_signal_filter_hard, 'Магнитометрический', 500)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "А теперь применим для сейсмического сигнала."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_filter = int(Fs / 2 / 50)\n",
    "width_window = 10\n",
    "\n",
    "seismic_signal_filter = seismic_signal\n",
    "seismic_signal_filter = butter_lowpass_filter(seismic_signal_filter, 250-width_window, Fs, order = 5)\n",
    "sub_plot_signal_spectr(seismic_signal_filter, 'Сейсмический', 500)\n",
    "\n",
    "for i in range(n_filter-1):\n",
    "    seismic_signal_filter = butter_bandstop_filter(seismic_signal_filter, [50*(i+1)-width_window, 50*(i+1)+width_window], Fs)\n",
    "    sub_plot_signal_spectr(seismic_signal_filter, 'Сейсмический', 500)\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "width_window = 10\n",
    "\n",
    "seismic_signal_filter_hard = seismic_signal_filter\n",
    "seismic_signal_filter_hard = butter_lowpass_filter(seismic_signal_filter_hard, 25, Fs)\n",
    "sub_plot_signal_spectr(seismic_signal_filter_hard, 'Магнитометрический', 500)\n",
    "show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Совместим две реализации"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_x = linspace(0,len(seismic_signal_filter_hard)/Fs, len(seismic_signal_filter_hard))\n",
    "plt.figure(figsize = (15, 6))\n",
    "plt.title('Совмещенные по времени реализации сейсмического и магнитометрического сигналов')\n",
    "plt.plot(time_x,seismic_signal_filter_hard, )\n",
    "plt.plot(time_x,magnito_signal_filter_hard, 'r')\n",
    "plt.legend(['Сейсмический сигнал', 'Магнитометрический сигнал'])\n",
    "plt.grid()\n",
    "show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}