{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "%reset -f\n",
    "\n",
    "# import various packages\n",
    "\n",
    "import numpy as np                    # matrices & etc\n",
    "from matplotlib import pyplot as plt  # general plotting \n",
    "from matplotlib import patches        # plotting shapes\n",
    "from datetime import date             # date and time arithmetic\n",
    "from math import exp, pi, sin, cos, tan, atan, atan2, sqrt, floor, log10, nan   # math functions\n",
    "import scipy.linalg as la             # linear algebra functions\n",
    "import os                             # operating system\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "import scipy.signal as sg\n",
    "import scipy.stats as st\n",
    "import scipy.sparse.linalg as las\n",
    "from scipy import sparse\n",
    "import matplotlib\n",
    "import pydub as pd\n",
    "import array\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def gda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = gda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if( N==0 ):\n",
    "            continue; # skip vector of zero length\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = gda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0:1] = v;  # patch 20230418 was #w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=gda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=gda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def gda_cvec1(v):\n",
    "    if isinstance(v, int) or isinstance(v, np.int32):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"gda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"gda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = np.array(v); # patch v -> np.array(v)\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"gda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"gda_cvec: %s not supported\" % type(v));\n",
    "\n",
    "        \n",
    "# gda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "def gda_draw(*argv):\n",
    "    DOCOLOR=True;\n",
    "    if( DOCOLOR ):\n",
    "        bwcmap = matplotlib.colormaps['jet'];\n",
    "    else:\n",
    "        bw = np.zeros((256,4));\n",
    "        v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "        bw[:,0] = v;\n",
    "        bw[:,1] = v;\n",
    "        bw[:,2] = v;\n",
    "        bw[:,3] = np.ones(256);\n",
    "        bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "# bandpass filter, used in seismological example, but hand\n",
    "# in a variety of settings involving time series\n",
    "def gda_chebyshevfilt(d, Dt, flow, fhigh):\n",
    "    # (dout,u,v)=gda_chebyshevfilt(d, Dt, flow, fhigh);\n",
    "    # chebyshev IIR bandpass filter\n",
    "    # d - input array of data\n",
    "    # Dt - sampling interval\n",
    "    # flow - low pass frequency, Hz\n",
    "    # fhigh - high pass frequency, Hz\n",
    "    # dout - output array of data\n",
    "    # u - the numerator filter\n",
    "    # v - the denominator filter\n",
    "    # these filters can be used again using dout=filter(u,v,din);\n",
    "\n",
    "    # make sure input timeseries is a column vector\n",
    "    s = np.shape(d);\n",
    "    N = s[0];\n",
    "    if(N==1):\n",
    "        dd = np.zeros((N,1));\n",
    "        dd[:,0] = d;\n",
    "    else:\n",
    "        dd=d;\n",
    "        \n",
    "    # sampling rate\n",
    "    rate=1/Dt;\n",
    "\n",
    "    # ripple parameter, set to ten percent\n",
    "    ripple=0.1;  \n",
    "\n",
    "    # normalise frequency\n",
    "    fl=2.0*flow/rate;\n",
    "    fh=2.0*fhigh/rate;\n",
    "\n",
    "    # center frequency \n",
    "    cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );\n",
    "\n",
    "    # bandwidth\n",
    "    bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );\n",
    "\n",
    "    # ripple parameter factor\n",
    "    rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));\n",
    "    a = 0.5*(rpf-1.0/rpf);\n",
    "    b = 0.5*(rpf+1.0/rpf);\n",
    "\n",
    "    u=np.zeros((5,1));\n",
    "    v=np.zeros((5,1));\n",
    "    theta = 3*pi/4;\n",
    "    sr = a * cos(theta);\n",
    "    si = b * sin(theta);\n",
    "    es = sqrt(sr*sr+si*si);\n",
    "    tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;\n",
    "    v[0,0] = 1.0;\n",
    "    v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;\n",
    "    v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;\n",
    "    v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    tmp = 4.0*es*es*bw*bw/tmp;\n",
    "    u[0,0] = tmp;\n",
    "    u[1,0] = 0.0;\n",
    "    u[2,0] = -2.0*tmp;\n",
    "    u[3,0] = 0.0;\n",
    "    u[4,0] = tmp;\n",
    "\n",
    "    dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());\n",
    "    return (gda_cvec(dout),gda_cvec(u),gda_cvec(v));\n",
    "\n",
    "# function to perform the multiplication FT (F v)\n",
    "def gdaFTFmul(v):\n",
    "    # this function is used by bicongugate gradient to\n",
    "    # solve the generalized least squares problem Fm=F.\n",
    "    # Note that \"F\" must be called \"gdaFsparse\".\n",
    "    global gdaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    # weirdly, \"*\" multiplies sparse matrices just fine\n",
    "    temp = gdaFsparse*vv;\n",
    "    return gdaFsparse.transpose()*temp;\n",
    "\n",
    "# function to perform the multiplication FT (G v)\n",
    "def gda_FTFmul(v):\n",
    "    # this function is used by bicongugate gradient to\n",
    "    # solve the generalized least squares problem Fm=F.\n",
    "    # Note that \"F\" must be called \"gdaFsparse\".\n",
    "    global gdaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    # weirdly, \"*\" multiplies sparse matrices just fine\n",
    "    temp = gdaFsparse*vv;\n",
    "    return gdaFsparse.transpose()*temp;\n",
    "\n",
    "\n",
    "# function to perform the multiplication FT f\n",
    "def gda_FTFrhs(f):\n",
    "    # this function is used by bicongugate gradient to\n",
    "    # solve the generalized least squares problem Fm=F.\n",
    "    # Note that \"F\" must be called \"gdaFsparse\".\n",
    "    global gdaFsparse;\n",
    "    return gdaFsparse.transpose()*f;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved office_data_Dt.txt with Dt 2.2676e-05\n",
      "Saved office_data_data.txt with 114688 samples\n",
      "Done\n"
     ]
    }
   ],
   "source": [
    "myfileroot = 'office_data';\n",
    "# myfileroot = 'plate_data';\n",
    "# myfileroot = 'whistle_data';\n",
    "# myfileroot = 'whistle_data';\n",
    "# myfileroot = 'teapot_data';\n",
    "# myfileroot = 'penny_data';\n",
    "\n",
    "myfileextension = 'm4a';\n",
    "\n",
    "myfileaudio=myfileroot+'.'+myfileextension;\n",
    "myfileDt=myfileroot+'_Dt.txt';\n",
    "myfiledata=myfileroot+'_data.txt';\n",
    "\n",
    "info = pd.utils.mediainfo(myfileaudio);\n",
    "rate = info['sample_rate'];\n",
    "Dt = 1.0 / float(rate);\n",
    "\n",
    "# write data file\n",
    "# concatenate time, discharge into a matrix C\n",
    "# save matrix to file of tab-separated values\n",
    "np.savetxt(myfileDt,gda_cvec(Dt), delimiter=\"\\t\");\n",
    "print(\"Saved %s with Dt %.4e\" % (myfileDt,Dt) );\n",
    "\n",
    "audio = pd.AudioSegment.from_file(myfileaudio,myfileextension);\n",
    "intdata = array.array('H',bytearray(audio.raw_data));\n",
    "audio = np.array(intdata,dtype=np.int16);\n",
    "N, = np.shape(audio);\n",
    "np.savetxt(myfiledata,audio.astype(float),delimiter=\"\\t\");\n",
    "print(\"Saved %s with %d samples\" % (myfiledata,N) );\n",
    "\n",
    "print('Done');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file office_data.m4a with rate 44100 and sampling 2.2676e-05\n",
      "file office_data.m4a with 114688 data\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 38140 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34032 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29724 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38388 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# reads a m4a file file made by the \"Voice Recorder\" app on my Galaxy cell phone\n",
    "# beats me how veratile it is.  Needed to install pydum with a\n",
    "# conda install -c conda-forge pydub\n",
    "# and then import pydub as pd and import array\n",
    "\n",
    "# myfile = 'tone_data.m4a';\n",
    "myfile = 'office_data.m4a';\n",
    "# myfile = 'whistle_data.m4a';\n",
    "# myfile = 'penny_data.m4a';\n",
    "# myfile = 'teapot_data.m4a';\n",
    "# myfile = 'plate_data.m4a';\n",
    "\n",
    "\n",
    "info = pd.utils.mediainfo(myfile);\n",
    "rate = info['sample_rate'];\n",
    "Dt = 1.0 / float(rate);\n",
    "print(\"file %s with rate %s and sampling %.4e\" % (myfile,rate,Dt) );\n",
    "\n",
    "audio = pd.AudioSegment.from_file(myfile,'m4a');\n",
    "intdata = array.array('H',bytearray(audio.raw_data));\n",
    "audio = np.array(intdata,dtype=np.int16);\n",
    "dobs = gda_cvec( audio.astype(float) );\n",
    "N, i = np.shape(dobs);\n",
    "if( N%2 ):\n",
    "    N=N-1;\n",
    "    dobs = dobs[0:N,0:1];\n",
    "    print(\"truncating to even, N=%d\" % (N) );\n",
    "print(\"file %s with %d data\" % (myfile,N) );\n",
    "\n",
    "t = gda_cvec( Dt*np.linspace(0,N-1,N) );\n",
    "tmin = t[0,0];\n",
    "tmax = t[N-1,0];\n",
    "dabsmax = np.max(np.abs(dobs));\n",
    "\n",
    "N1 = 0;\n",
    "N2 = N;\n",
    "fig1 = plt.figure();\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [t[N1,0], t[N2-1,0], -1.1*dabsmax, 1.1*dabsmax] );\n",
    "plt.plot(t[N1:N2,0:1],dobs[N1:N2,0:1],'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('d');\n",
    "plt.show();\n",
    "\n",
    "# generic frequency up\n",
    "No2 = floor(N/2);  # integer version of N/2\n",
    "fmax=1/(2.0*Dt);   # nyquist frequency \n",
    "Df=fmax/No2;       # frequency smaping\n",
    "Nf=No2+1;          # number on non-negative frequencies\n",
    "f=gda_cvec(Df*np.concatenate(( np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2*pi*Df;        # angular frequeny increment\n",
    "w=2*pi*f;          # full (pos and neg) ana\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=gda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos;    # non-negative angular frequencies\n",
    "\n",
    "# remove mean\n",
    "dobs = dobs - np.mean(dobs);\n",
    "\n",
    "#  Hamming window function\n",
    "W = gda_cvec( 0.54 - 0.46*np.cos(2*pi*np.linspace(0,N-1,N)/(N-1)) );\n",
    "dtapered = np.multiply( W, dobs );\n",
    "\n",
    "fig2 = plt.figure();\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [t[N1,0], t[N2-1,0], -1.1*dabsmax, 1.1*dabsmax] );\n",
    "plt.plot(t[N1:N2,0:1],dtapered[N1:N2,0:1],'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('d');\n",
    "plt.show();\n",
    "\n",
    "#  compute fourier coefficients. Note normalization of Dt\n",
    "dbar = Dt*np.fft.fft(dtapered, axis=0); # fourier transform\n",
    "dbarpos = dbar[0:Nf,0:1]; # non-negative frequencies\n",
    "T = N*Dt;  # time series length\n",
    "s2 = (2/T) * np.real(np.multiply(dbarpos.conj(), dbarpos));\n",
    "\n",
    "fplotmin=0.0;\n",
    "fplotmax=fmax;\n",
    "# linear in frequency plot\n",
    "fig3 = plt.figure(figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([fplotmin, fplotmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in frequency');\n",
    "plt.show();\n",
    "\n",
    "fplotmin=0.0;\n",
    "fplotmax=120.0;\n",
    "# linear in frequency plot\n",
    "fig4 = plt.figure(figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([fplotmin, fplotmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in frequency');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file office_data_Dt.txt with sampling 2.2676e-05\n",
      "file office_data_data.txt with 114688 data\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 38140 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34032 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29724 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38388 characters]",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# reads a m4a file file made by the \"Voice Recorder\" app on my Galaxy cell phone\n",
    "# beats me how veratile it is.  Needed to install pydum with a\n",
    "# conda install -c conda-forge pydub\n",
    "# and then import pydub as pd and import array\n",
    "\n",
    "myfileroot = 'office_data';\n",
    "# myfileroot = 'plate_data';\n",
    "# myfileroot = 'whistle_data';\n",
    "# myfileroot = 'whistle_data';\n",
    "# myfileroot = 'teapot_data';\n",
    "# myfileroot = 'penny_data';\n",
    ";\n",
    "myfileextension = 'm4a';\n",
    "\n",
    "myfileaudio=myfileroot+'.'+myfileextension;\n",
    "myfileDt=myfileroot+'_Dt.txt';\n",
    "myfiledata=myfileroot+'_data.txt';\n",
    "\n",
    "# get data from file\n",
    "Dt = np.genfromtxt(myfileDt, delimiter='\\t');\n",
    "print(\"file %s with sampling %.4e\" % (myfileDt,Dt) );\n",
    "\n",
    "audio = np.genfromtxt(myfiledata, delimiter='\\t');\n",
    "dobs = gda_cvec( audio.astype(float) );\n",
    "N, i = np.shape(dobs);\n",
    "if( N%2 ):\n",
    "    N=N-1;\n",
    "    dobs = dobs[0:N,0:1];\n",
    "    print(\"truncating to even, N=%d\" % (N) );\n",
    "print(\"file %s with %d data\" % (myfiledata,N) );\n",
    "\n",
    "t = gda_cvec( Dt*np.linspace(0,N-1,N) );\n",
    "tmin = t[0,0];\n",
    "tmax = t[N-1,0];\n",
    "dabsmax = np.max(np.abs(dobs));\n",
    "\n",
    "N1 = 0;\n",
    "N2 = N;\n",
    "fig1 = plt.figure();\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [t[N1,0], t[N2-1,0], -1.1*dabsmax, 1.1*dabsmax] );\n",
    "plt.plot(t[N1:N2,0:1],dobs[N1:N2,0:1],'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('d');\n",
    "plt.show();\n",
    "\n",
    "# generic frequency up\n",
    "No2 = floor(N/2);  # integer version of N/2\n",
    "fmax=1/(2.0*Dt);   # nyquist frequency \n",
    "Df=fmax/No2;       # frequency smaping\n",
    "Nf=No2+1;          # number on non-negative frequencies\n",
    "f=gda_cvec(Df*np.concatenate(( np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2*pi*Df;        # angular frequeny increment\n",
    "w=2*pi*f;          # full (pos and neg) ana\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=gda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos;    # non-negative angular frequencies\n",
    "\n",
    "# remove mean\n",
    "dobs = dobs - np.mean(dobs);\n",
    "\n",
    "#  Hamming window function\n",
    "W = gda_cvec( 0.54 - 0.46*np.cos(2*pi*np.linspace(0,N-1,N)/(N-1)) );\n",
    "dtapered = np.multiply( W, dobs );\n",
    "\n",
    "fig2 = plt.figure();\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [t[N1,0], t[N2-1,0], -1.1*dabsmax, 1.1*dabsmax] );\n",
    "plt.plot(t[N1:N2,0:1],dtapered[N1:N2,0:1],'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('d');\n",
    "plt.show();\n",
    "\n",
    "#  compute fourier coefficients. Note normalization of Dt\n",
    "dbar = Dt*np.fft.fft(dtapered, axis=0); # fourier transform\n",
    "dbarpos = dbar[0:Nf,0:1]; # non-negative frequencies\n",
    "T = N*Dt;  # time series length\n",
    "s2 = (2/T) * np.real(np.multiply(dbarpos.conj(), dbarpos));\n",
    "\n",
    "fplotmin=0.0;\n",
    "fplotmax=fmax;\n",
    "# linear in frequency plot\n",
    "fig3 = plt.figure(figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([fplotmin, fplotmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in frequency');\n",
    "plt.show();\n",
    "\n",
    "fplotmin=0.0;\n",
    "fplotmax=120.0;\n",
    "# linear in frequency plot\n",
    "fig4 = plt.figure(figsize=(10,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([fplotmin, fplotmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in frequency');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
