Skip to content

Commit

Permalink
experiments/piano_heterodyne: Optimize synth_harmonics so full analys…
Browse files Browse the repository at this point in the history 
…is takes 2 min not 15 min.
  • Loading branch information
@dpwe
dpwe committed Feb 3, 2025
1 parent db6e1a6 commit d522358
Showing 1 changed file with 46 additions and 26 deletions.
72 changes: 46 additions & 26 deletions experiments/piano_heterodyne.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,22 @@
"print(replicate_last(np.array([[1,2,3,4],[5,6,7,8]]), 5))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "292a19ac-6fce-477e-8a69-ca8d40fa238d",
"metadata": {},
"outputs": [],
"source": [
"w = np.hanning(256)\n",
"x = np.random.rand(1000)\n",
"print(np.convolve(w, x, 'same').shape)\n",
"print(scipy.signal.convolve(x, w, 'same').shape)\n",
"# \"same\" means \"same as 1st arg\" for scipy, but \"same as longer arg\" for numpy.\n",
"#convolve = np.convolve\n",
"convolve = scipy.signal.convolve"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -89,7 +105,7 @@
" interpolation_window = np.hstack([np.linspace(0, 1, 1 + interp)[1:], np.linspace(1, 0, 1 + interp)[1:]])\n",
" for i in range(x_in_flat.shape[0]):\n",
" x_atrous = np.reshape(np.hstack([x_in_flat[[i], :].T, np.zeros((n_cols, interp - 1))]), (n_cols * interp,))\n",
" x_out_flat[i] = np.convolve(x_atrous, interpolation_window, 'same')\n",
" x_out_flat[i] = convolve(x_atrous, interpolation_window, 'same')\n",
" if truncate:\n",
" if interp > 1: # Skip for degenerate interp == 1; indexing by [:, :-(1 - 1)] deletes the entire array.\n",
" x_out_flat = x_out_flat[:, :-(interp - 1)]\n",
Expand All @@ -108,7 +124,7 @@
" interpolation_window = np.expand_dims(np.hstack([np.linspace(0, 1, 1 + interp)[1:], np.linspace(1, 0, 1 + interp)[1:]]), 0)\n",
" for i in range(x_in_flat.shape[-1]):\n",
" #x_atrous = np.reshape(np.hstack([x_in_flat[[i], :].T, np.zeros((n_cols, interp - 1))]), (n_cols * interp,))\n",
" #x_out_flat[i] = np.convolve(x_atrous, interpolation_window, 'same')\n",
" #x_out_flat[i] = convolve(x_atrous, interpolation_window, 'same')\n",
" x_out_flat[..., i * interp : (i + 2) * interp] += x_in_flat[..., i] * interpolation_window\n",
" x_out_flat = x_out_flat[..., (interp - 1):-1]\n",
" if truncate:\n",
Expand Down Expand Up @@ -213,15 +229,17 @@
" assert len(freq) == len(mag)\n",
" # Ignore the final freq value, which tends to go crazy. Replicate the preceding one.\n",
" #freq[-1] = freq[-2]\n",
" freq = slinterp(freq, interp_factor, truncate=truncate)\n",
" if interp_factor > 1:\n",
" freq = slinterp(freq, interp_factor, truncate=truncate)\n",
" carrier = np.cos(np.cumsum(2 * np.pi / sample_rate * freq))\n",
" if carrier is None:\n",
" # freq must be a scalar\n",
" #freq = freq * np.ones(1 + (len(mag) - 1) * interp_factor)\n",
" carrier = np.cos(2 * np.pi * freq / sample_rate * np.arange(truncate + (len(mag) - truncate) * interp_factor))\n",
"\n",
" # Can now expand mag to match.\n",
" mag = slinterp(mag, interp_factor, truncate=truncate)\n",
" if interp_factor > 1:\n",
" mag = slinterp(mag, interp_factor, truncate=truncate)\n",
" # freq is cycles per sec. We want radians per sample\n",
" if mag_is_db:\n",
" mag = db_to_lin(mag)\n",
Expand Down Expand Up @@ -471,7 +489,7 @@
" n_bins = np.sum(freqs < f_max)\n",
" X_spec = 20 * np.log10(np.abs(long_fft[:n_bins]))\n",
" # Local average\n",
" smoo_X_spec = np.convolve(np.ones(local_smooth_points) / local_smooth_points, X_spec, 'same')\n",
" smoo_X_spec = convolve(X_spec, np.ones(local_smooth_points) / local_smooth_points, 'same')\n",
" unsmoo_X_spec = X_spec - smoo_X_spec\n",
"\n",
" #peaks = scipy.signal.argrelextrema(np.maximum(np.max(unsmoo_X_spec) - 30, unsmoo_X_spec), np.greater)[0]\n",
Expand Down Expand Up @@ -785,29 +803,14 @@
"`heterodyne_extract` extracts sample-by-sample magnitude and frequency contours for a single harmonic in the neighborhood of a specified input frequency (specified by a fundamental and a real-valued harmonic number, so the actual center frequency is simply the product of the two). This is done by 'heterodyne extraction' - frequency-shifting the signal by multiplying it by the conjugate complex sinusoid at the center frequency (which shifts the center frequency down to 0 Hz), then low-pass filtering with a hanning window whose length is `smooth_cycles` times the fundamental period corresponding to `fundamental_freq`. This window has zeros at the fundamental, so modulation at the fundamental rate is completely canceled. This gives smooth contours for true harmonics of a fundamental."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "292a19ac-6fce-477e-8a69-ca8d40fa238d",
"metadata": {},
"outputs": [],
"source": [
"w = np.hanning(256)\n",
"x = np.random.rand(1000)\n",
"print(np.convolve(w, x, 'same').shape)\n",
"print(scipy.signal.convolve(x, w, 'same').shape)\n",
"# \"same\" means \"same as 1st arg\" for scipy, but \"same as longer arg\" for numpy."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5278f393-a5b0-48ee-8b23-74293225a815",
"metadata": {},
"outputs": [],
"source": [
"#convolve = np.convolve\n",
"convolve = scipy.signal.convolve\n",
"#convolve = scipy.signal.convolve\n",
"\n",
"def heterodyne_extract(waveform, fundamental_freq, harmonic_number=1.0, sr=44100, smooth_cycles=2.0):\n",
" \"\"\"Return a full-sample-rate amplitude and frequency envelope near the specified frequency.\"\"\"\n",
Expand Down Expand Up @@ -1118,6 +1121,19 @@
"Audio(data=(audio).T, rate=sr)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "89829014-3f7e-4f39-8242-7d5d2a31b043",
"metadata": {},
"outputs": [],
"source": [
"f = avg_freqs[0]\n",
"ms = mags[0]\n",
"%timeit a = synth_harmonic(f, ms, 1)\n",
"plt.plot(synth_harmonic(f, ms[::2], 2))"
]
},
{
"cell_type": "markdown",
"id": "cf340f1a-993c-4447-b993-ef65bdf74627",
Expand Down Expand Up @@ -1227,7 +1243,7 @@
"trim_end_samps = 1024\n",
"\n",
"hnum = 1\n",
"mag_filt = lin_to_db(np.convolve(filt_win, db_to_lin(mags[hnum]), 'same'))\n",
"mag_filt = lin_to_db(convolve(db_to_lin(mags[hnum]), filt_win, 'same'))\n",
"mag_filt_trim = mag_filt[trim_start_samps:-trim_end_samps]\n",
"max_times = scipy.signal.argrelmax(mag_filt_trim)[0]\n",
"max_vals = mag_filt_trim[max_times]\n",
Expand Down Expand Up @@ -1279,7 +1295,8 @@
"# harms_params is a list of harmonic descriptions\n",
"# each harmonic description is [const_freq, [(delta_time_ms, lin_mag), (delta_time_ms, lin_mag), ...]]\n",
"\n",
"def make_harms_params(mags, avg_freqs, max_harmonic=None, sr=44100, mag_floor=-100.0, terminal_slope=20.0, adaptive_timing=True, verbose=False):\n",
"def make_harms_params(mags, avg_freqs, max_harmonic=None, sr=44100, mag_floor=-100.0, terminal_slope=20.0, \n",
" filt_len=31, adaptive_timing=True, verbose=False):\n",
" if not max_harmonic:\n",
" max_harmonic = mags.shape[0]\n",
" \n",
Expand All @@ -1290,11 +1307,14 @@
" #trim_ends_samps = int(round(trim_ends_sec * sr))\n",
" trim_start_samps = 0 # 128\n",
" trim_end_samps = 1024\n",
" \n",
"\n",
" #filt_len = 31\n",
" filt_win = np.hanning(filt_len)/np.sum(np.hanning(filt_len))\n",
"\n",
" harms_params = []\n",
"\n",
" for hnum in range(max_harmonic):\n",
" mag_filt = lin_to_db(np.convolve(filt_win, db_to_lin(mags[hnum]), 'same'))\n",
" mag_filt = lin_to_db(convolve(db_to_lin(mags[hnum]), filt_win, 'same'))\n",
" mag_filt_trim = mag_filt[trim_start_samps : -trim_end_samps]\n",
" #mag_filt_trim = mag_filt\n",
" mmag = np.maximum(mag_floor, mag_filt_trim)\n",
Expand Down Expand Up @@ -1404,7 +1424,7 @@
" h_decay = np.zeros(num_harmonics)\n",
" h_nums = range(num_harmonics)\n",
" for h in h_nums:\n",
" mag_filt = lin_to_db(np.convolve(filt_win, db_to_lin(mags[h]), 'same'))\n",
" mag_filt = lin_to_db(convolve(db_to_lin(mags[h]), filt_win, 'same'))\n",
" mag_filt_trim = mag_filt[trim_end_samps:-trim_end_samps]\n",
" t = np.arange(len(mag_filt_trim)) / sr\n",
" linfit = lin_fit(mag_filt_trim)\n",
Expand Down

0 comments on commit d522358

Please sign in to comment.