{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4b1b36d9",
   "metadata": {},
   "source": [
    "# Wave turbulence in gongs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cc6b20c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import signal as sig\n",
    "from scipy import fft as fft\n",
    "#import scipy as sc\n",
    "import scipy.io.wavfile as wav\n",
    "%matplotlib inline\n",
    "#%matplotlib notebook  \n",
    "#for jupiter notebook\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9796e975",
   "metadata": {},
   "source": [
    "# Acoustic measurement\n",
    "\n",
    "A microphone is available to record the sound emitted by the gong. Use the Audacity sofware to record the sound. Save the data in .wav file format that python can read using the scipy library.\n",
    "\n",
    "A quasi-stationary forcing can be achieved by hitting reapeatedly the gong with constant impact strength.\n",
    "\n",
    "#### 1 Perform a series of experiments with different forcing strengths (and possibly different gongs). Document precisely for dataset. Make long records."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4feddd1",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d669b204",
   "metadata": {},
   "source": [
    "#### 2. load one audio record (use the function `wav.read` from `scipy` to read .wav files). There are two channels but we will use only one (whichever). Plot the signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "259a631a",
   "metadata": {},
   "outputs": [],
   "source": [
    "rate, data = wav.read('backup2.wav')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2693e0e4",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a08054cd",
   "metadata": {},
   "source": [
    "#### 3. Stationary regime\n",
    "Compute the Fourier spectrum of the data, in the quasi-stationary regime, using the Welch method (`welch` function in Scipy, read the documentation of this function). Plot the data (choose between linear of logarithmic scales). Compare the different intensities and the different gongs). Optimize the values of the different parameters of the `welch`method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d047ddd2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "fb952b26",
   "metadata": {},
   "source": [
    "#### 4. Transitory regime\n",
    "To study the transitory part of the signal (at the beginning), one needs to compute the evolution of the frequency content as a function of time. The adapted tool for this is the spectrogram. It consists in computing the Fourier transform of successive short time wondows of the signal. \n",
    "\n",
    "Compute the spectrogram of one of your signals (use the `spectrogram` function from `scipy`, read the doc). Show the result in a figure. Optimise the values of the parameters to make the figure more readable, especially at the beginning of the relevant sound. To make the computation faster, it is relevant to restrict the computation to the most interesting part of the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8dc58f6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dd51378",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "39e38a06",
   "metadata": {},
   "source": [
    "#### Compare the different experiments and different gongs. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f3da2a7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "04ffdb30",
   "metadata": {},
   "source": [
    "# Spatial measurement: time-resolved profilometry \n",
    "\n",
    "In order to understand more the acoustic observation, we would like to have a spatial measurement of vibration. For this we will use the Fourier Transform Profilometry technique described by the article by Maurel et al. \"Experimental and theoretical inspection of the phase-to-height relation in Fourier transform profilometry\" Applied Optics 2009. Ask your professor to explain you the principle.\n",
    "\n",
    "#### 1. Optimise the setup to perform the measurement on the biggest gong. Record a few images with static gong and with a homogeneous image.\n",
    "\n",
    "#### 2. Parameters of the method\n",
    "Identify the parameters of the method that you have to measure to be able to process the images to obtain the deformation of the gong. Give the values. Explain how you choose the frame rate of the camera."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a940fd98",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "78d8a704",
   "metadata": {},
   "source": [
    "#### 3. Define a protocol to be able to record the motion in the stationary regime. Record the movie. Give the settings of the camera. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0538addb",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3216cfb4",
   "metadata": {},
   "source": [
    "#### 4. Ask your professor to show you the various steps of the image processing. The actual processing takes time and will be perform later on"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a67484c",
   "metadata": {},
   "source": [
    "#### 5. Load the data from the profilometry. You can use the `netCDF4` library.  The signal is a 3 dimensional matrix $h(x,y,t)$. At some point, you will have to concatenate the data from several files to have a long enough duration of the loaded data (use the ` concatenate` function from `numpy`). We are rather interested in the normal velocity signal rather than the deformation. Compute the velocity by simple finite difference. Show a snapshot of the velocity at a given time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebf5e575",
   "metadata": {},
   "outputs": [],
   "source": [
    "import netCDF4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45a97690",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for this part\n",
    "# Nb files to concatenate\n",
    "Nf=16\n",
    "# index of first file\n",
    "N1=1\n",
    "# number of windows for Welch method\n",
    "Nw=4;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8841d6bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "fname=\"imt\"+f\"{N1}\"+\".nc\"\n",
    "print(fname)\n",
    "nc = netCDF4.Dataset(fname)\n",
    "h=nc.variables['h']\n",
    "#the dimensions are flipped compared to matlab: t,x,y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ca57f4d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8ee33698",
   "metadata": {},
   "source": [
    "#### We want to compute the Fourier spectrum in both $(\\vec k,\\omega)$ space. \n",
    "For this we will use the Welch method and compute  $$S(\\vec k,\\omega)=\\langle |\\hat v(\\vec k,\\omega)|^2\\rangle$$ The Fourier transform in space will be performed over all available space but the Fourier transform in time will be performed over a window of finite duration $T$ that you have to choose. $T$ must be long enough to get a good resolution in frequency. The average $\\langle \\cdot \\rangle$ is then performed over successive time windows.\n",
    "\n",
    "As the signal is not periodic neither in space and time, a Hanning (or other) window must be used in all dimensions (space and time). The 3D Hanning window can be computed from 1D Hanning windows simply by $$h_{3D}(x,y,t)=h_x(x)h_y(y)h_t(t)$$ where $h_i$ is the Hanning window in the $i$ dimension. The spectrum is then computed by performing the Fourier transform of $v(x,y,t)h_{3D}(x,y,t)$.\n",
    "\n",
    "#### 6. Compute the 3D hanning window for your choice of $T$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a226978d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for this part\n",
    "# Nb files to concatenate\n",
    "Nf=8\n",
    "# index of first file\n",
    "N1=1\n",
    "# number of windows for Welch method\n",
    "Nw=25;\n",
    "\n",
    "Lx=65;\n",
    "Ly=50;\n",
    "\n",
    "hx=np.hanning(Lx)\n",
    "hy=np.hanning(Ly)\n",
    "ht=np.hanning(Nf*500-1)\n",
    "h2D=np.zeros((Lx,Ly))\n",
    "for i in range(0,Lx):\n",
    "    h2D[i,:]=hy*hx[i]\n",
    "h3D=np.zeros((Nf*500-1,Lx,Ly))\n",
    "for i in range(0,Nf*500-1):\n",
    "    h3D[i,:,:]=h2D*ht[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "803e4ee3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "827f2a1d",
   "metadata": {},
   "source": [
    "#### 7. By loading successive time windows of the signal (at least 4), compute the spectrum $S(\\vec k,\\omega)$.  For a faster computation, you should use the `fft` function from `scipy` with the option `workers=-1` or  use `workers=c` were `c` is the number of logical processors of your computer (for example `fft.fftn(v*h3D,workers=-1`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "863b35ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "85e04065",
   "metadata": {},
   "source": [
    "#### 8. Plot a slice of the spectrum (which is a 3D matrix) to show that the vibration is indeed made of waves. The dispersion relation of the bending waves of a flat plate is $\\omega=Ck^2$. Is it consistent with your observation ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64b8caff",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "beff771a",
   "metadata": {},
   "source": [
    "#### 9. From the 3D spectrum, compute the frequency only spectrum $S(\\omega)$. Plot it in a figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c55d7328",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfb89ef5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "190e6267",
   "metadata": {},
   "source": [
    "The spectrum is made of a large collection of peaks. To understand the reason why the spectrum is not continuous as expected from the weak turbulence theory, we would like to see the deformation associated to each one of these peaks. For this we can compute $\\langle|v(x,y,\\omega)|.^2\\rangle$ using the Welch method. In this case, a Hanning function can be used only in the time dimension.\n",
    "\n",
    "#### 10. Why should one use the modulus square and not simply $\\langle v(x,y,\\omega)\\rangle$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "415b7368",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0ddd29ad",
   "metadata": {},
   "source": [
    "#### 11. Compute the relevant Hanning window and $\\langle|v(x,y,\\omega)|.^2\\rangle$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b723bf18",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27e6e100",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ce4614a1",
   "metadata": {},
   "source": [
    "#### 12. For the lowest frequency peaks in the spectrum, show the result. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0f94dd7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e99de3e7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7896a673",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}