Python IIR filter to Pure Data using biquad~.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d191504e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import IPython.display as ipd\n",
    "import librosa as lib\n",
    "import subprocess\n",
    "import soundfile as sf\n",
    "import librosa.display\n",
    "\n",
    "%matplotlib inline\n",
    "%config IPCompleter.greedy=True\n",
    "\n",
    "sr=48000 # Defining the sample rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b58d74b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def apply_filter_python(filt, input_file, output_file):\n",
    "    \"\"\"Applies filter using lfilter2() in SciPy\"\"\"\n",
    "    audio = lib.load(input_file, sr=sr)[0] # Loads the audio file to which the filter will be applied\n",
    "    python_filtered = signal.lfilter(b=filt[0], a=filt[1], x=audio) # Applies the filter\n",
    "    sf.write(python_output_file, python_filtered, sr) # Writes the filtered signal to the specific output audio file\n",
    "    \n",
    "    print(\"Python filtering complete. Filtered output saved to:\", output_file)\n",
    "    \n",
    "    return python_filtered # Returns the filtered audio for visualisation\n",
    "\n",
    "def apply_filter_pd(filt, input_file, output_file):\n",
    "    \"\"\"Breaks down the IIR filter into second order sections, passes these to Pure Data to be used in biquad~\n",
    "    and then applies the filter in Pure data from the specified input file to the specified output file.\"\"\"\n",
    "\n",
    "    # Breaks down IIR filter into second-order section coefficients\n",
    "    sos = signal.tf2sos(filt[0], filt[1], pairing=\"nearest\")\n",
    "    # Scales coefficients by dividing all by a0 (a0 is usually 1 anyway, but just in case) and inverting a1 and a2\n",
    "    sos_scaled = __pd_scale_sos_coeffs(sos)\n",
    "\n",
    "    command = __pd_build_command(sos_scaled, input_file, output_file) # Builds the command line command\n",
    "\n",
    "    subprocess.run(command) # Actually runs the command\n",
    "    print(\"Pure Data filtering complete. Filtered output saved to:\", output_file) # Only prints once Pd is closed\n",
    "    \n",
    "    pd_filtered = lib.load(output_file, sr=sr)[0] # Loads the filtered audio file generated by Pure Data\n",
    "    \n",
    "    return pd_filtered # Returns the filtered audio for visualisation\n",
    "\n",
    "def __pd_build_command(sos_scaled, input_file, output_file):\n",
    "    \"\"\"Builds the command used to open Pure Data and apply the filter\"\"\"\n",
    "    \n",
    "    pd_executable = \"C:/Program Files/Pd/bin/pd.exe\" # <- will need changing depending on Pure Data install location!\n",
    "    pd_patch = \"tf2sos_biquad.pd\"\n",
    "    \n",
    "    section_strings = [] # Initialised empty list which will contains strings of coefficients\n",
    "    sends = ' -send \"' # Initialises string containing all sends\n",
    "    coeff_order = [4,5,0,1,2] # List of new order for coefficients i.e., fb1/a1, fb2/a2, ff1/b0, ff2/b1, ff3/b2 as biquad~ in PD needs them\n",
    "\n",
    "    for sec_i, section in enumerate(sos_scaled):\n",
    "        coeff_string = \"\"\n",
    "        for i in coeff_order: # iterates through section in order defined by coeff_order list - this should put them in the right order for Pd\n",
    "            coeff_string = coeff_string + str(section[i]) + \" \"\n",
    "        section_string = \"; sec\" + str(sec_i + 1) + \" \" + coeff_string[:-1] # Appends together Pd message for each section\n",
    "        sends = sends + section_string # Builds up sends string from section_string\n",
    "\n",
    "    # Finishes sends string by adding input and output filenames\n",
    "    sends = sends + \"; input_file \" + input_file + \"; output_file \" + output_file + '\"'\n",
    "\n",
    "    # Assembles the command\n",
    "    command = pd_executable + \" -open \" + pd_patch + \" -r \" + str(sr) + sends #+ \" -nogui -batch\" # \n",
    "    # command = pd_executable + \" -open \" + pd_patch + \" -r \" + str(sr) + sends + \" -nogui -batch\" <- runs Pure Data in batch mode w/o the GUI\n",
    "        \n",
    "    return command\n",
    "        \n",
    "def __pd_scale_sos_coeffs(sos):\n",
    "    \"\"\"Formats coefficients by dividing every coefficient by the a0 coefficient\"\"\"\n",
    "    sos_scaled = np.array([]) # Initialises empty array for scaled sos (every coefficient divided by its section's a0)\n",
    "    for section in sos:\n",
    "        section_scaled = np.array([])\n",
    "        for i, coeff in enumerate(section): # iterates through coefficients in section\n",
    "            coeff_scaled = coeff/section[3] # Divides coefficient by index 3, i.e., the a0 coefficient\n",
    "            if i >= 4:\n",
    "                coeff_scaled = -coeff_scaled # Inverts a1 and a2 coefficients to -a1 and -a2\n",
    "            section_scaled = np.append(section_scaled, coeff_scaled) # Appends scaled coeff to scaled section\n",
    "        sos_scaled = np.append(sos_scaled, section_scaled) # Appends scaled section to scaled sos\n",
    "\n",
    "    # Reshapes sos_scaled array so the input and output have the same shape i.e., sections of 6 coefficients\n",
    "    sos_scaled = np.reshape(sos_scaled, newshape=(-1, 6))\n",
    "\n",
    "    return sos_scaled\n",
    "    \n",
    "def spectrogram(s, title):\n",
    "    \"\"\"Generates a spectrogram of a given signal.\"\"\"\n",
    "    # Generates data\n",
    "    S = lib.stft(s)\n",
    "    S_dB = lib.amplitude_to_db(abs(S), ref=np.max)\n",
    "\n",
    "    # Creates figure\n",
    "    fig, ax = plt.subplots(figsize=(12,4))\n",
    "\n",
    "    # Adds data to figure\n",
    "    spectro = lib.display.specshow(S_dB, sr=sr, x_axis=\"time\", y_axis=\"log\")\n",
    "\n",
    "    # Figure setup\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(\"Time (s)\")\n",
    "    ax.set_ylabel(\"Frequency (Hz)\")\n",
    "    ax.set_ylim(0, sr*0.5)\n",
    "\n",
    "    # Adds colourbar to the spectrogram\n",
    "    plt.grid(False)\n",
    "    plt.colorbar(spectro, label=\"Gain (dB)\")\n",
    "\n",
    "    # Displays figure\n",
    "    img_save_loc = \"./spectrograms/\" + title + \"_spectrogram.png\"\n",
    "    plt.savefig(img_save_loc)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "843cb55c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Designing an IIR filter using SciPy signal\n",
    "filt = signal.iirdesign(wp=0.05, ws=0.07, gpass=1, gstop=20, ftype=\"cheby2\")\n",
    "\n",
    "input_file = \"./test_audio/to_be_filtered.wav\" # Defines input file location\n",
    "python_output_file = \"./test_audio/filtered_python.wav\" # Defines Python output file location\n",
    "pd_output_file = \"./test_audio/filtered_pd.wav\" # Defines Pure Data output file location\n",
    "\n",
    "# Before filtering\n",
    "before_filtering = lib.load(input_file, sr=sr)[0] # Loads unfiltered audio file\n",
    "spectrogram(before_filtering, \"Audio before filtering\") # Displays spectrogram of unfilitered signal\n",
    "display(ipd.Audio(before_filtering, rate=sr)) # Allows playback of unfiltered signal\n",
    "\n",
    "# Python filtering\n",
    "python_filtered = apply_filter_python(filt, input_file, python_output_file) # Applies filter to test file using SciPy\n",
    "spectrogram(python_filtered, \"Audio filtered by Python\") # Displays spectrogram of the filtered signal\n",
    "display(ipd.Audio(python_filtered, rate=sr)) # Allows playback of the filtered audio\n",
    "\n",
    "# Pure Data filtering\n",
    "pd_filtered = apply_filter_pd(filt, input_file, pd_output_file) # Applies filter to test file using Pure Data\n",
    "spectrogram(pd_filtered, \"Audio filtered by Pure Data\") # Displays spectrogram of the filtered signal\n",
    "display(ipd.Audio(python_filtered, rate=sr)) # Allows playback of the filtered audio"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('base')",
   "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "33e85536a753e4289d1e3ca1904fea47091c61e030b0b2f85895695a2354c09d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}