notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import soundfile as sf\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import random\n",
    "import sklearn\n",
    "import librosa\n",
    "import librosa.display\n",
    "\n",
    "from scipy.io import wavfile\n",
    "from scipy import signal\n",
    "from scipy.fft import fft, fftfreq\n",
    "from sklearn.decomposition import FastICA\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (15,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_animal = \"Data Sets/Bird Song/\"\n",
    "dir_music  = \"Data Sets/GTZAN/\"\n",
    "dir_speech = \"Data Sets/TIMIT/\"\n",
    "dir_envir  = \"Data Sets/TUT Acoustic Scenes\"\n",
    "dir_synths = \"Data Sets/Synthesizers\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_chosen = dir_speech\n",
    "\n",
    "# For music, bird song, and acoustic: Fs = 22050, 44100, and 44100\n",
    "# For speech: Fs = 16000\n",
    "num_samples = 800\n",
    "num_sounds  = 256 #len(os.listdir(dir_chosen)) \n",
    "\n",
    "X = np.zeros((num_sounds, num_samples))\n",
    "\n",
    "l = []\n",
    "\n",
    "while (len(l) < num_sounds):\n",
    "    temp = random.randrange(0, len(os.listdir(dir_chosen)), 1)\n",
    "    \n",
    "    if (temp in l):\n",
    "        continue\n",
    "    \n",
    "    else:\n",
    "        l.append(temp)\n",
    "\n",
    "index = 0\n",
    "for idx, file in enumerate(os.listdir(dir_chosen)):\n",
    "    filename = os.fsdecode(file)\n",
    "    \n",
    "    if (idx not in l):\n",
    "        continue\n",
    "    \n",
    "    if filename.endswith(\".wav\") or filename.endswith(\".WAV\"): \n",
    "        Fs, data = wavfile.read(dir_chosen + \"/\" + filename)\n",
    "                                            \n",
    "        if (len(data.shape) > 1):\n",
    "            data = (data[:, 0] + data[:, 1]) / 2\n",
    "        \n",
    "        start = random.randrange(0, len(data) - num_samples, 1)\n",
    "        \n",
    "        hp_filt  = signal.butter(16, Fs / num_samples, btype='highpass', output='sos', fs=Fs)\n",
    "        filtdata = signal.sosfilt(hp_filt, data)\n",
    "        \n",
    "        data = data[start:start + num_samples]\n",
    "        \n",
    "        if (index == 0):\n",
    "            plt.plot(data)\n",
    "            plt.title(\"Original Signal Segement\")\n",
    "            plt.show()\n",
    "        \n",
    "        if (index == 0):\n",
    "            plt.plot(filtdata[start:start + num_samples])\n",
    "            plt.title(\"Signal After High-Pass Filter\")\n",
    "            plt.show()\n",
    "        \n",
    "        X[index, :] = filtdata[start:start + num_samples]\n",
    "        \n",
    "        index = index + 1\n",
    "        \n",
    "    else:\n",
    "        continue\n",
    "        \n",
    "print(\"All signals high-pass filtered at frequency {}. {} samples of length {} randomly\"\n",
    "      \" selected for study.\".format(Fs / num_samples, num_sounds, num_samples))\n",
    "\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ica(x, samplerate, sources):\n",
    "    \n",
    "    tolerance = 1e-20\n",
    "    iterations = 2000\n",
    "    \n",
    "    # Data points\n",
    "    n = np.size(x, 1)\n",
    "\n",
    "    # Initialize and normalize W\n",
    "    w = np.random.rand(sources, np.size(x, 0))\n",
    "    w = w/np.sqrt(np.sum(np.power(w, 2), axis=1, keepdims=True))\n",
    "\n",
    "    # Centering and whitening\n",
    "    print('Whitening data..')\n",
    "    \n",
    "    # Rows\n",
    "    x_mean = x.mean(1, keepdims=True)\n",
    "    x_centered = x-x_mean\n",
    "    \n",
    "    # Whitening\n",
    "    r_cov = np.cov(x_centered)\n",
    "    u_cov, s_cov, v_cov = np.linalg.svd(r_cov, hermitian=True)\n",
    "    v = u_cov @ (np.diag(1/np.sqrt(s_cov)) @ u_cov.T)\n",
    "    z = v @ x_centered\n",
    "\n",
    "    # ICA loop\n",
    "    print('ICA loop..')\n",
    "    \n",
    "    for i in range(iterations):\n",
    "        # Previous W\n",
    "        wlast = w.copy()\n",
    "        \n",
    "        # Current source estimate\n",
    "        y = w @ z\n",
    "\n",
    "        # Gradient function\n",
    "        g = y*np.exp(-0.5*np.power(y, 2))\n",
    "        dg = (1-np.power(y, 2))*np.exp(-0.5*np.power(y, 2))\n",
    "\n",
    "        # Update W\n",
    "        w = (g @ z.transpose())/n - dg.mean(axis=1, keepdims=True)*w\n",
    "        \n",
    "        # Normalize W\n",
    "        w = w/np.sqrt(np.sum(np.power(w, 2), axis=1, keepdims=True))\n",
    "        \n",
    "        # Decorrelate W\n",
    "        u_w, s_w, v_w = np.linalg.svd(w)\n",
    "        w = u_w @ np.diag(1/s_w) @ u_w.T @ w\n",
    "\n",
    "        # Check convergence\n",
    "        delta = np.max(1-np.abs(np.sum(w*wlast, axis=1)))\n",
    "        \n",
    "        if np.mod(i, 50) == 0:\n",
    "            print(i, '=>', delta)\n",
    "            \n",
    "        if delta < tolerance:\n",
    "            y = w @ z\n",
    "            print(i, '=>', delta)\n",
    "            print('Converged.')\n",
    "            return y, w, v, x_mean\n",
    "\n",
    "    # Output\n",
    "    y = w @ z\n",
    "    print('Did not converge.')\n",
    "    return y, w, v, x_mean\n",
    "\n",
    "\n",
    "# Call ica function\n",
    "ICs, w, v, mu = ica(X, Fs, sources=128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ICs.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (15,2)\n",
    "\n",
    "filename = \"Results/Synthesizers/ICs/IC_\"\n",
    "\n",
    "for x in range(128):\n",
    "    plt.plot(ICs[x, :])\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "    #plt.savefig(filename + str(x), dpi=200)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (15,2)\n",
    "\n",
    "filename = \"Results/Synthesizers/Specs/Spec_\"\n",
    "\n",
    "centers = []\n",
    "downs   = []\n",
    "ups     = []\n",
    "\n",
    "for idx in range(128):\n",
    "\n",
    "    f, t, Sxx = signal.spectrogram(ICs[idx, :], Fs, nperseg=16)\n",
    "    plt.pcolormesh(t, f, Sxx, shading='gouraud')\n",
    "    plt.ylabel('Frequency (Hz)')\n",
    "    plt.xlabel('Time (s)')\n",
    "    #plt.savefig(filename + str(idx), dpi=200)\n",
    "    plt.show()\n",
    "\n",
    "    ps = np.square(np.abs(np.fft.fft(ICs[idx, :])))[0:128]\n",
    "    \n",
    "    center_flag = False\n",
    "    down_flag   = False\n",
    "    up_flag     = False\n",
    "    \n",
    "    summ = 0\n",
    "\n",
    "    for x in range(len(ps)):\n",
    "        temp = ps[x]\n",
    "        summ += temp\n",
    "\n",
    "        if (summ > ps.sum() / 2 and center_flag == False):\n",
    "            center = x * Fs / (len(ps) * 2)\n",
    "            center_flag = True\n",
    "\n",
    "        if (summ > ps.sum() / 4 and down_flag == False):\n",
    "            down = x * Fs / (len(ps) * 2)\n",
    "            down_flag = True\n",
    "\n",
    "        if (summ > ps.sum() * 3 / 4 and up_flag == False):\n",
    "            up = x * Fs / (len(ps) * 2)\n",
    "            up_flag = True\n",
    "    \n",
    "    centers.append(center)\n",
    "    downs.append(down)\n",
    "    ups.append(up)\n",
    "    \n",
    "ups   = np.asarray(ups)\n",
    "downs = np.asarray(downs)\n",
    "centers    = np.asarray(centers)\n",
    "bandwidths = ups - downs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (10,10)\n",
    "\n",
    "plt.plot(centers[:], bandwidths[:], '*')\n",
    "plt.xlabel('Center Frequency (Hz)')\n",
    "plt.ylabel('Width (Hz)')\n",
    "plt.title('Synthesizers Sound Filter Characteristics')\n",
    "#plt.savefig(\"Results/Synthesizers/ch_filter\", dpi=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}