{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e208aa51-4e26-4487-b80d-df83fdc00277",
   "metadata": {},
   "source": [
    "# make_piano_examples\n",
    "\n",
    "Quick notebook to compose sound examples for my planned blog post about piano synthesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d83c04c3-309e-4364-a3cc-eb3d31ee2977",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import time\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.io.wavfile as wav\n",
    "import scipy.signal\n",
    "\n",
    "from IPython.display import Audio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ec3aed7-be52-486d-974c-b927df746c7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spectrogram display, copied from piano_heterodyne.\n",
    "\n",
    "from scipy import fft\n",
    "\n",
    "def dB(x):\n",
    "    return 20 * np.log10(np.abs(x))\n",
    "\n",
    "def spectrogram(x, fft_len=1024, win_len=None, hop_len=None, window_fn=np.hanning, sr=1, f_max=None, title=None):\n",
    "    if win_len is None:\n",
    "        win_len = fft_len\n",
    "    if hop_len is None:\n",
    "        hop_len = win_len // 2\n",
    "    # Window.\n",
    "    prepad_len = (fft_len - win_len) // 2\n",
    "    win = np.hstack([np.zeros(prepad_len), window_fn(win_len), np.zeros(fft_len - win_len - prepad_len)])\n",
    "    # Frame.\n",
    "    frame_indices = np.arange(0, len(x) - win_len, hop_len)[:, np.newaxis] + np.arange(win_len)[np.newaxis, :]\n",
    "    x_chunks_windowed = x[frame_indices] * win[np.newaxis, :]\n",
    "    # Transform.\n",
    "    stft_mag_db = dB(fft.fft(x_chunks_windowed, n=fft_len)[:, :(fft_len // 2 + 1)])\n",
    "    num_frames, num_bins = stft_mag_db.shape\n",
    "    t_base = np.arange(num_frames) * hop_len / sr\n",
    "    f_base = np.arange(num_bins) * sr / fft_len\n",
    "    plt.imshow(stft_mag_db.T, aspect='auto', origin='lower', extent=[t_base[0], t_base[-1], f_base[0], f_base[-1]])\n",
    "    plt.clim(np.max(stft_mag_db) + [-80, 0])\n",
    "    #plt.ylim([0, f_max / sr * fft_len])\n",
    "    if f_max:\n",
    "        plt.ylim([0, f_max])\n",
    "    #x_times = np.arange(0, len(x) / sr, 0.25)\n",
    "    #y_freqs = np.arange(0, 10000, 1000)\n",
    "    #plt.xticks(x_times * sr / hop_len, x_times)\n",
    "    #plt.yticks(y_freqs / sr * fft_len, y_freqs)\n",
    "    plt.colorbar(label='level / dB')\n",
    "    plt.xlabel('time / sec')\n",
    "    plt.ylabel('freq / kHz')\n",
    "    if title:\n",
    "        plt.title(title)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d509e7f7-6f86-456b-a157-f2016154b793",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in the original piano samples, copied from piano_heterodyne\n",
    "\n",
    "def wavread(filename):\n",
    "  \"\"\"Read in audio data from a wav file.  Return d, sr.\"\"\"\n",
    "  # Read in wav file.\n",
    "  file_handle = open(filename, 'rb')\n",
    "  samplerate, wave_data = wav.read(file_handle)\n",
    "  # Normalize short ints to floats in range [-1..1).\n",
    "  data = np.asarray(wave_data, dtype=np.float32) / 32768.0\n",
    "  return data, samplerate\n",
    "\n",
    "def wavwrite(data, samplerate, filename):\n",
    "  \"\"\"Write a waveform to a WAV file.\"\"\"\n",
    "  wav.write(filename, samplerate, (32768.0 * data).astype(np.int16))\n",
    "\n",
    "def read_and_trim(filename, duration=2.0, rel_threshold=0.015, abs_threshold=0.0006, sr=44100, channel=0, do_plot=False, pre_time=0.002):\n",
    "    d, file_sr = wavread(filename)\n",
    "    #print(\"d.shape=\", d.shape)\n",
    "    if len(d.shape) > 1:  # Stereo file\n",
    "        if channel is None:\n",
    "            d = np.mean(d, axis=-1)\n",
    "        else:\n",
    "            d = d[:, channel]\n",
    "    assert file_sr == sr\n",
    "    d = d - np.mean(d)\n",
    "    threshold = np.maximum(np.max(np.abs(d)) * rel_threshold, abs_threshold)\n",
    "    # Tight window at start?\n",
    "    drop_initial_samps = np.maximum(0, -int(round(pre_time * sr)) + np.min(np.flatnonzero(np.abs(d) > threshold)))\n",
    "    d_return = d[int(round(drop_initial_samps)) + np.arange(int(round(duration * sr)))]\n",
    "    if do_plot:\n",
    "        t = np.arange(len(d)) / sr\n",
    "        plt.figure(figsize=(12, 4))\n",
    "        plt.subplot(121)\n",
    "        plt.plot(t[:50000], d[:50000])\n",
    "        plt.ylim(threshold * np.array([-1, 1]))\n",
    "        plt.subplot(122)\n",
    "        plt.plot(t[:len(d_return)], d_return)\n",
    "        plt.xlim([0, 0.1])\n",
    "        print(\"drop initial\", drop_initial_samps / sr)\n",
    "    return d_return\n",
    "\n",
    "def piano_filename(note='C', octave=4, strike='ff'):\n",
    "    #filename = 'Piano.ff.D4.wav'\n",
    "    directory = '/Users/dpwe/Downloads/uiowa-piano'\n",
    "    #octave = 2\n",
    "    #note = 'D'\n",
    "    #octave = 4\n",
    "    #note = 'Eb'\n",
    "    #strike = 'ff'  # ['pp', 'mf', 'ff']:\n",
    "    filename = os.path.join(directory, '.'.join(['Piano', strike, note + str(octave), 'wav']))\n",
    "    return filename\n",
    "\n",
    "filename = piano_filename('D', 4)\n",
    "print(filename)\n",
    "waveform = read_and_trim(filename)\n",
    "_, sr = wavread(filename)\n",
    "Audio(data=waveform.T, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9a44a35-b787-433d-aaee-4f9525d5d178",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Give it a 50ms fade-out\n",
    "def fade_out(waveform, fade_end_sec=0.5, fade_duration_sec=0.05, sample_rate=44100):\n",
    "    \"\"\"Apply a fade-out to a waveform.\"\"\"\n",
    "    fade_end_samples = int(round(fade_end_sec * sample_rate))\n",
    "    fade_duration_samples = int(round(fade_duration_sec * sample_rate))\n",
    "    fade_start_samples = fade_end_samples - fade_duration_samples\n",
    "    mask = np.ones_like(waveform)\n",
    "    mask[fade_start_samples : fade_start_samples + fade_duration_samples] = np.linspace(1, 0, fade_duration_samples)\n",
    "    mask[fade_start_samples + fade_duration_samples:] = 0\n",
    "    return mask * waveform\n",
    "\n",
    "fade_waveform = fade_out(waveform)\n",
    "#plt.plot(fade_waveform)\n",
    "Audio(data=fade_waveform.T, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbcfeb1a-8d63-448b-9588-502a8a19340a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compose the 3s (time, note, vel) sequence from piano_examples.py\n",
    "\n",
    "def emplace(main_waveform, waveform_to_add, start_time_sec, sample_rate=44100):\n",
    "    \"\"\"Add-in a waveform to a background at a given time, truncate as needed.\"\"\"\n",
    "    start_time_samples = int(round(start_time_sec * sample_rate))\n",
    "    length_to_use = min(len(waveform_to_add), len(main_waveform) - start_time_samples)\n",
    "    main_waveform[start_time_samples : start_time_samples + length_to_use] += waveform_to_add[:length_to_use]\n",
    "    return main_waveform\n",
    "\n",
    "sr = 44100\n",
    "waveform = np.zeros(int(round(3.3 * sr)))\n",
    "\n",
    "waveform = emplace(waveform, fade_out(read_and_trim(piano_filename('D', 4, 'pp')), 0.4), 0.05)\n",
    "waveform = emplace(waveform, fade_out(read_and_trim(piano_filename('D', 4, 'mf')), 0.4), 0.45)\n",
    "waveform = emplace(waveform, fade_out(read_and_trim(piano_filename('D', 4, 'ff')), 0.65), 0.85)\n",
    "\n",
    "waveform = emplace(waveform, fade_out(read_and_trim(piano_filename('D', 2, 'mf')), 1.8), 1.5)\n",
    "waveform = emplace(waveform, fade_out(read_and_trim(piano_filename('D', 6, 'ff')), 1.2), 2.1)\n",
    "\n",
    "Audio(data=waveform.T, rate=sr)                    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cdb36a0-50fa-46c0-8be0-308cf549eb13",
   "metadata": {},
   "outputs": [],
   "source": [
    "wavwrite(waveform, sr, '../sounds/piano_example_original.wav')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dad756fb-08e0-4653-8fa4-c7d495328f2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '../sounds/piano_example_original.wav'\n",
    "waveform, sr = wavread(filename)\n",
    "plt.figure(figsize=(16, 4))\n",
    "spectrogram(waveform, sr=sr, f_max=8000, title='UIowa Samples')\n",
    "Audio(data=waveform, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79ab2480-bd38-4b53-b95e-83dd7c85dc45",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '../sounds//piano_example_juno_patch_7.wav'\n",
    "waveform, sr = wavread(filename)\n",
    "plt.figure(figsize=(16, 4))\n",
    "spectrogram(waveform, sr=sr, f_max=8000, title='Juno-60 piano')\n",
    "Audio(data=waveform, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1436640c-56ca-4049-8ff2-d645ea39abe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '../sounds//piano_example_dx7_patch_137.wav'\n",
    "waveform, sr = wavread(filename)\n",
    "plt.figure(figsize=(16, 4))\n",
    "spectrogram(waveform, sr=sr, f_max=8000, title='DX7 piano')\n",
    "Audio(data=waveform, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27ad0bd1-7b2d-45ab-9902-591cfd0d4691",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '../sounds/piano_example_additive_fixed.wav'\n",
    "waveform, sr = wavread(filename)\n",
    "plt.figure(figsize=(16, 4))\n",
    "spectrogram(waveform, sr=sr, f_max=8000, title='Fixed-envelope additive piano')\n",
    "Audio(data=waveform, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06ff335b-4aba-437e-837b-ccbde61498ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '../sounds//piano_example_additive_interpolated.wav'\n",
    "waveform, sr = wavread(filename)\n",
    "plt.figure(figsize=(16, 4))\n",
    "spectrogram(waveform, sr=sr, f_max=8000, title='Interpolated additive piano')\n",
    "Audio(data=waveform, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0ee56a1-8391-42c8-b1ac-c9c193a6d02f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeat original again\n",
    "filename = '../sounds//piano_example_original.wav'\n",
    "waveform, sr = wavread(filename)\n",
    "plt.figure(figsize=(16, 4))\n",
    "spectrogram(waveform, sr=sr, f_max=8000, title='UIowa Samples')\n",
    "Audio(data=waveform, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5825360-1d12-4694-b40d-b051406c7395",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot FFT to illustrate two-sided spectrum\n",
    "from scipy import fft\n",
    "\n",
    "def plotspec(d, n_fft=4096, sr=44100, f_max=6000, title=None, db_range=80):\n",
    "    X = fft.fft(d[:n_fft] * np.hanning(n_fft))\n",
    "    # FFT shift to put zero in the middle.\n",
    "    fft_keep_points = int(round(f_max / sr * n_fft))\n",
    "    X = np.hstack([X[-fft_keep_points:], X[:fft_keep_points]])\n",
    "    freqs = np.arange(-fft_keep_points, fft_keep_points) * sr / n_fft / 1000\n",
    "    db_X = dB(X)\n",
    "    plt.plot(freqs, db_X)\n",
    "    db_max = 10 * np.ceil(np.max(db_X) / 10)\n",
    "    plt.ylim([db_max - db_range, db_max])\n",
    "    plt.xlabel('freq / kHz')\n",
    "    plt.ylabel('level / dB')\n",
    "    plt.grid()\n",
    "    if title:\n",
    "        plt.title(title)\n",
    "\n",
    "d = read_and_trim(piano_filename('D', 4, 'ff'))\n",
    "# Skip first 0.2 s to exclude \"thump\"\n",
    "plt.figure(figsize=(16, 4))\n",
    "plotspec(d[int(round(0.2 * sr)):], title='D4.ff spectrum', f_max=3000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05d60937-faa2-4c92-b6be-9e7d088fdbfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Demodulate by multiplying by complex exponential\n",
    "f0 = 293.665  # Nominal frequency of D4\n",
    "d_demod = d * np.exp(-1j * 2 * np.pi * f0 * np.arange(len(d))/ sr)\n",
    "plotspec(d_demod[int(round(0.2 * sr)):] , title='D4.ff shifted down by f0', f_max=3000)\n",
    "f0_on_2 = f0 / 2000\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.plot([-f0_on_2, f0_on_2, f0_on_2, -f0_on_2, -f0_on_2], [38, 38, -38, -38, 38], '--r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f65bcd09-1940-40db-ba04-ef1b01ff726d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Low-pass demodulated waveform to recover harmonic envelope.\n",
    "\n",
    "def demod_and_smooth(d, f0, harmonic=1, sr=44100):\n",
    "    d_demod = d * np.exp(-1j * 2 * np.pi * f0 * harmonic * np.arange(len(d))/ sr)\n",
    "    fundamental_period_samples = int(round(sr / f0))\n",
    "    d_demod_smoo = np.convolve(d_demod, np.hanning(2 * fundamental_period_samples), 'same')\n",
    "    return d_demod_smoo\n",
    "\n",
    "t = np.arange(len(d)) / sr\n",
    "plt.figure(figsize=(16, 6))\n",
    "\n",
    "plt.subplot(211)\n",
    "plt.plot(t, d)\n",
    "plt.ylabel('amplitude')\n",
    "plt.title('D4.ff waveform')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(t, dB(demod_and_smooth(d, f0)))\n",
    "plt.plot(t, dB(demod_and_smooth(d, f0, harmonic=2)))\n",
    "plt.plot(t, dB(demod_and_smooth(d, f0, harmonic=3)))\n",
    "plt.plot(t, dB(demod_and_smooth(d, f0, harmonic=4)))\n",
    "plt.ylim([-15, 25])\n",
    "plt.title('Extracted harmonic envelopes')\n",
    "plt.xlabel('time / sec')\n",
    "plt.ylabel('level / dB')\n",
    "plt.legend(['%d - %.2f Hz' % (n, n * f0) for n in range(1, 5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df559105-e346-45eb-bf66-7ee59c6da5b3",
   "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.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
