Skip to content

Commit

Permalink
Merge pull request #293 from shorepine/heterodyne
Browse files Browse the repository at this point in the history 
Heterodyne analysis optimization
  • Loading branch information
@dpwe
dpwe authored Feb 3, 2025
2 parents 54b05fa + d522358 commit a811a26
Show file tree
Hide file tree
Showing 4 changed files with 93 additions and 43 deletions.
2 changes: 1 addition & 1 deletion amy_headers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -354,7 +354,7 @@ def make_piano_patch():
amy.send(chorus=0) # Piano sounds weird with chorus on
amy.send(osc=0, wave=amy.INTERP_PARTIALS, patch=0)
amy.send(osc=20, wave=amy.PARTIAL)
return 21
return 25 # We now use up to 24 partials per voice + 1 control osc.

def make_patches(filename):
def nothing(message):
Expand Down
92 changes: 64 additions & 28 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 Expand Up @@ -1703,8 +1723,13 @@
"\n",
"#audio, ep = synth_params([harms_params[i] for i in 1 + np.arange(10)])\n",
"#filename = 'Piano.ff.C4.wav'\n",
"filename = 'Piano.ff.C1.wav'\n",
"h_audio, ep = synth_params(harms_dict[filename][:40])\n",
"filename = 'Piano.ff.C3.wav'\n",
"#hd = harms_dict[filename][:40]\n",
"# How does it sound if we skip some of the higher harmonics?\n",
"# In particular, above the 16th harmonic, skip harmonics that are 2.n.f0 and 3.n.f0. This lets us get higher, for brighter sounds, \n",
"# and the missing harmonics aren't particularly noticeable in those regions.\n",
"hd = [harms_dict[filename][i] for i in [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 22, 24, 28, 30, 34, 36]]\n",
"h_audio, ep = synth_params(hd)\n",
"#h_audio, ep = synth_params(harms_params)\n",
"\n",
"#f0 = harms_params[0][0]\n",
Expand All @@ -1726,6 +1751,17 @@
"# (start_pt + np.arange(npts)) / sr, old_audio[start_pt : start_pt + npts])\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad4e2638-a2e3-48ac-a0af-e26085bb9646",
"metadata": {},
"outputs": [],
"source": [
"print(filename, len(audio), sr)\n",
"Audio(data=audio, rate=sr)"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down
40 changes: 27 additions & 13 deletions src/interp_partials.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,16 @@ typedef struct {

#define MAX_NUM_MAGNITUDES 24

#define MAX_NUM_HARMONICS 20
#define MAX_NUM_HARMONICS 40

// Map to drop out some higher harmonics, namely the 2x and 3x overtones above 16th harmonic
typedef uint8_t bool;
const bool use_this_partial_map[MAX_NUM_HARMONICS] = {
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1-10
1, 1, 1, 1, 1, 1, 1, 0, 1, 0, // 11-20
0, 0, 1, 0, 1, 0, 0, 0, 1, 0, // 21-30
1, 0, 0, 0, 1, 0, 1, 0, 0, 0, // 31-40
};

void _cumulate_scaled_harmonic_params(float *harm_param, int harmonic_index, float alpha, const interp_partials_voice_t *partials_voice) {
int num_bps = partials_voice->num_sample_times_ms;
Expand Down Expand Up @@ -131,24 +140,29 @@ void interp_partials_note_on(uint16_t osc) {
// harmonic_base_index_pl_vl, pitch_alpha, vel_alpha,
// alpha_pl_vl, alpha_pl_vh, alpha_ph_vl, alpha_ph_vh);
for (int h = 0; h < num_harmonics; ++h) {
for (int i = 0; i < MAX_NUM_MAGNITUDES + 1; ++i) harm_param[i] = 0;
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_pl_vl + h,
alpha_pl_vl, partials_voice);
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_pl_vh + h,
alpha_pl_vh, partials_voice);
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_ph_vl + h,
alpha_ph_vl, partials_voice);
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_ph_vh + h,
alpha_ph_vh, partials_voice);
//fprintf(stderr, "harm %d freq %.2f bps %.3f %.3f %.3f %.3f\n", h, harm_param[0], harm_param[1], harm_param[2], harm_param[3], harm_param[4]);
_osc_on_with_harm_param(osc + 1 + h, harm_param, partials_voice);
if (use_this_partial_map[h]) {
for (int i = 0; i < MAX_NUM_MAGNITUDES + 1; ++i) harm_param[i] = 0;
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_pl_vl + h,
alpha_pl_vl, partials_voice);
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_pl_vh + h,
alpha_pl_vh, partials_voice);
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_ph_vl + h,
alpha_ph_vl, partials_voice);
_cumulate_scaled_harmonic_params(harm_param, harmonic_base_index_ph_vh + h,
alpha_ph_vh, partials_voice);
//fprintf(stderr, "harm %d freq %.2f bps %.3f %.3f %.3f %.3f\n", h, harm_param[0], harm_param[1], harm_param[2], harm_param[3], harm_param[4]);
++osc;
_osc_on_with_harm_param(osc, harm_param, partials_voice);
}
}
}

void interp_partials_note_off(uint16_t osc) {
//const interp_partials_voice_t *partials_voice = &interp_partials_map[synth[osc].patch % NUM_INTERP_PARTIALS_PATCHES];
//int num_oscs = partials_voice->num_harmonics[0]; // Assume first patch has the max #harmonics.
int num_oscs = MAX_NUM_HARMONICS;
int num_oscs = 0; //MAX_NUM_HARMONICS;
// Actual max num harmonics we may use is the number of 1s in the use_this_partial_map.
for (int i = 0; i < MAX_NUM_HARMONICS; ++i) num_oscs += use_this_partial_map[i];
for(uint16_t i = osc + 1; i < osc + 1 + num_oscs; i++) {
uint16_t o = i % AMY_OSCS;
AMY_UNSET(synth[o].note_on_clock);
Expand Down
2 changes: 1 addition & 1 deletion src/patches.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,6 +262,6 @@ static const char * const patch_commands[257] PROGMEM = {
/* 256: dpwe piano */ "k0Zv0w12p0Zv20w9Z",
};
const uint16_t patch_oscs[257] PROGMEM = {
5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,21,
5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,25,
};
#endif

0 comments on commit a811a26

Please sign in to comment.