FFT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FFT windowing\n",
    "---\n",
    "\n",
    "\n",
    "This notebook visualizes the FFT windows with different windowing functions. To start working with this notebook, you need to download the [IEEE dataset](https://ieee-dataport.org/open-access/experimental-database-detecting-and-diagnosing-rotor-broken-bar-three-phase-induction), unpack it and place in the folder named *IEEE* above the current folder, containing the notebook. Also, you should have the libraries from *requirements.txt* installed in the virtual environment, as described in the README file.\n",
    "\n",
    "The plots are made with interactive module to enable the inspection of the window cycle number, motor state and load on the current signal frequency spectrum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Libraries import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import h5py\n",
    "import scipy.fft\n",
    "from ipywidgets import interact\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "data loading"
    ]
   },
   "outputs": [],
   "source": [
    "# array of data[class[0BRB-4BRB]][load[0.125-1]]\n",
    "# data[0][0] - 0BRB, 12.5% load\n",
    "# data[3][2] - 3BRB, 37.5% load etc.\n",
    "exp_data=[[j for j in range(8)] for i in range(5)]\n",
    "load_levels=['torque05','torque10','torque15','torque20','torque25','torque30','torque35','torque40']\n",
    "states=['rs','r1b','r2b','r3b','r4b']\n",
    "i=0\n",
    "F_sampling=50000 #Hz, current value sampling\n",
    "for state in states:\n",
    "    j=0\n",
    "    for level in load_levels:\n",
    "        mat=h5py.File('../IEEE/struct_'+state+'_R1.mat','r')\n",
    "        dataset_refname=mat[state][level]['Ia'][0]\n",
    "        data=mat[dataset_refname[0]]\n",
    "        data=np.ravel(data)\n",
    "        data=data[120000:] # discards the transient part of the signals\n",
    "        d={'Current, A':data}\n",
    "        exp_data[i][j]=pd.DataFrame(data=d,index=np.ravel([k/F_sampling for k in range (len(data))]))\n",
    "        exp_data[i][j].index.name='Time, s'\n",
    "        j+=1\n",
    "    mat.close()\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Shape of the data: {np.shape(exp_data)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_data[0][0].head() #array of data[class[0BRB-3BRB]][load[0-1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Current signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "@interact(brb=widgets.IntSlider(description='BRB',style=style,min=0,max=4,step=1,value=0),\n",
    "          load=widgets.FloatSlider(description='Load level, %',style=style,min=12.5,max=100,step=12.5,value=12.5))\n",
    "def current_signal(brb,load):\n",
    "    load=int(load/12.5-1)\n",
    "    plt.rcParams['figure.dpi']=300\n",
    "    fig, ax=plt.subplots(figsize=(16,7))\n",
    "    plt.plot(exp_data[brb][load]['Current, A'][:0.06],linewidth=1)\n",
    "    ax.set_xlabel('Time, s', fontsize=18)\n",
    "    ax.set_ylabel('Current, A', fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply windowing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rectangular"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "@interact(brb=widgets.IntSlider(description='BRB',style=style,min=0,max=4,step=1,value=0),\n",
    "          load=widgets.FloatSlider(description='Load level, %',style=style,min=12.5,max=100,step=12.5,value=12.5),\n",
    "          window_length=widgets.IntSlider(description='Window length',style=style,min=1,max=200,step=1,value=40))\n",
    "def fft_hamming(brb,load,window_length):\n",
    "    #signal info\n",
    "    F_sampling=50000 #Hz\n",
    "    resolution=1/F_sampling\n",
    "    F_grid=60 #Hz, grid frequency\n",
    "    load=int(load/12.5-1)\n",
    "    data=exp_data[brb][load][:window_length/F_grid].copy()\n",
    "    N=np.shape(data)[0]\n",
    "    ywf_rect=scipy.fft.rfft(np.ravel(data['Current, A']))\n",
    "    yfreq=scipy.fft.fftfreq(N,d=resolution)[:N//2]\n",
    "    fig,ax=plt.subplots(figsize=(16,7))\n",
    "    plt.semilogy(yfreq[1:N//2],2/N*np.abs(ywf_rect[1:N//2]));\n",
    "    ax.set_xlim([0,120])\n",
    "    ax.set_xlabel('Frequency, Hz',fontsize=18)\n",
    "    ax.set_ylabel('Amplitude, dB',fontsize=18)\n",
    "    ax.set_title('Rectangular window',fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Blackman"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "@interact(brb=widgets.IntSlider(description='BRB',style=style,min=0,max=4,step=1,value=0),\n",
    "          load=widgets.FloatSlider(description='Load level, %',style=style,min=12.5,max=100,step=12.5,value=12.5),\n",
    "          window_length=widgets.IntSlider(description='Window length',style=style,min=1,max=200,step=1,value=40))\n",
    "def fft_blackman(brb,load,window_length):\n",
    "    #signal info\n",
    "    plt.rcParams['figure.dpi']=300\n",
    "    F_sampling=50000 #Hz\n",
    "    resolution=1/F_sampling\n",
    "    F_s=60 #supply frequency\n",
    "    load=int(load/12.5-1)\n",
    "    data=exp_data[brb][load][:window_length/F_s].copy()\n",
    "    window=scipy.signal.windows.blackman(len(data))\n",
    "    #signal length\n",
    "    N=np.shape(data)[0]\n",
    "    ywf=scipy.fft.fft(np.ravel(data['Current, A'])*window)\n",
    "    yfreq=scipy.fft.fftfreq(N,d=resolution)[:N//2]\n",
    "    fig,ax=plt.subplots(figsize=(16,7));\n",
    "    plt.semilogy(yfreq[1:N//2],2/N*np.abs(ywf[1:N//2]));\n",
    "    ax.set_xlim([0,120])\n",
    "    ax.set_xlabel('Frequency, Hz', fontsize=18)\n",
    "    ax.set_ylabel('Amplitude, dB', fontsize=18)\n",
    "    ax.set_title('Blackman window',fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hanning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "@interact(brb=widgets.IntSlider(description='BRB',style=style,min=0,max=4,step=1,value=0),\n",
    "          load=widgets.FloatSlider(description='Load level, %',style=style,min=12.5,max=100,step=12.5,value=12.5),\n",
    "          window_length=widgets.IntSlider(description='Window length',style=style,min=1,max=200,step=1,value=40))\n",
    "def fft_hanning(brb,load,window_length):\n",
    "    #signal info\n",
    "    F_sampling=50000 #Hz\n",
    "    resolution=1/F_sampling\n",
    "    F_s=60 #grid frequency\n",
    "    load=int(load/12.5-1)\n",
    "    data=exp_data[brb][load][:window_length/F_s].copy()\n",
    "    N=np.shape(data)[0]\n",
    "    window_hanning=np.hanning(len(data['Current, A']))\n",
    "    ywf_hanning=scipy.fft.rfft(np.ravel(data['Current, A'])*window_hanning)\n",
    "    yfreq=scipy.fft.fftfreq(N,d=resolution)[:N//2]\n",
    "    fig,ax=plt.subplots(figsize=(16,7))\n",
    "    plt.semilogy(yfreq[1:N//2],2/N*np.abs(ywf_hanning[1:N//2]))\n",
    "    ax.set_xlim([0,120])\n",
    "    ax.set_xlabel('Frequency, Hz',fontsize=18)\n",
    "    ax.set_ylabel('Amplitude, dB',fontsize=18)\n",
    "    ax.set_title(f'Hanning window',fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hamming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "@interact(brb=widgets.IntSlider(description='BRB',style=style,min=0,max=4,step=1,value=0),\n",
    "          load=widgets.FloatSlider(description='Load level, %',style=style,min=12.5,max=100,step=12.5,value=12.5),\n",
    "          window_length=widgets.IntSlider(description='Window length',style=style,min=1,max=200,step=1,value=40))\n",
    "def fft_hamming(brb,load,window_length):\n",
    "    #signal info\n",
    "    F_sampling=50000 #Hz\n",
    "    resolution=1/F_sampling\n",
    "    F_s=60 #Hz, grid frequency\n",
    "    load=int(load/12.5-1)\n",
    "    data=exp_data[brb][load][:window_length/F_s].copy()\n",
    "    N=np.shape(data)[0]\n",
    "    window_hamming=np.hamming(len(data['Current, A']))\n",
    "    ywf_hamming=scipy.fft.rfft(np.ravel(data['Current, A'])*window_hamming)\n",
    "    yfreq=scipy.fft.fftfreq(N,d=resolution)[:N//2]\n",
    "    fig,ax=plt.subplots(figsize=(16,7))\n",
    "    plt.semilogy(yfreq[1:N//2],2/N*np.abs(ywf_hamming[1:N//2]))\n",
    "    ax.set_xlim([0,120])\n",
    "    ax.set_xlabel('Frequency, Hz',fontsize=18)\n",
    "    ax.set_ylabel('Amplitude, dB',fontsize=18)\n",
    "    ax.set_title('Hamming window',fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Flattop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "@interact(brb=widgets.IntSlider(description='BRB',style=style,min=0,max=4,step=1,value=0),\n",
    "          load=widgets.FloatSlider(description='Load level, %',style=style,min=12.5,max=100,step=12.5,value=12.5),\n",
    "          window_length=widgets.IntSlider(description='Window length',style=style,min=1,max=200,step=1,value=40))\n",
    "def fft_hamming(brb,load,window_length):\n",
    "    #signal info\n",
    "    F_sampling=50000 #Hz\n",
    "    resolution=1/F_sampling\n",
    "    F_s=60 #Hz, grid frequency\n",
    "    load=int(load/12.5-1)\n",
    "    data=exp_data[brb][load][:window_length/F_s].copy()\n",
    "    N=np.shape(data)[0]\n",
    "    window_flattop=scipy.signal.windows.flattop(len(data['Current, A']))\n",
    "    ywf_flattop=scipy.fft.rfft(np.ravel(data['Current, A'])*window_flattop)\n",
    "    yfreq=scipy.fft.fftfreq(N,d=resolution)[:N//2]\n",
    "    fig,ax=plt.subplots(figsize=(16,7))\n",
    "    plt.semilogy(yfreq[1:N//2],2/N*np.abs(ywf_flattop[1:N//2]))\n",
    "    ax.set_xlim([0,120])\n",
    "    ax.set_xlabel('Frequency, Hz',fontsize=18)\n",
    "    ax.set_ylabel('Amplitude, dB',fontsize=18)\n",
    "    ax.set_title('Flattop window',fontsize=20)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.10.6"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}