aboutsummaryrefslogtreecommitdiff
path: root/src/hrtfcont.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/hrtfcont.cpp')
-rw-r--r--src/hrtfcont.cpp256
1 files changed, 256 insertions, 0 deletions
diff --git a/src/hrtfcont.cpp b/src/hrtfcont.cpp
new file mode 100644
index 0000000..7e018f0
--- /dev/null
+++ b/src/hrtfcont.cpp
@@ -0,0 +1,256 @@
+/*
+ cw_binaural~: a binaural synthesis external for pure data
+ by David Doukhan - david.doukhan@gmail.com - http://www.limsi.fr/Individu/doukhan
+ and Anne Sedes - sedes.anne@gmail.com
+ Copyright (C) 2009-2011 David Doukhan and Anne Sedes
+
+ For more details, see CW_binaural~, a binaural synthesis external for Pure Data
+ David Doukhan and Anne Sedes, PDCON09
+
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <iostream>
+#include <math.h>
+#include "hrtfcont.hpp"
+#include "logstring.hpp"
+
+#define DEG2RAD (M_PIl/180)
+#define RAD2DEG (180/M_PIl)
+#define EPS 0.0000000000001
+
+HrtfCont::HrtfCont(const ir_key& k):
+ _vert_pol_coords(k.vertical_polar_coords) {}
+
+// set elevation in [-90, +90], and azimuth in [-180, 180[
+void HrtfCont::normalize_vertpolar_coords(float& azimuth, float& elevation) const
+{
+ // convert any elevation into a usefull value
+ // the value used should lie between -90 and +90
+ elevation = fmod(elevation, 360);
+ if (elevation < 0)
+ elevation += 360;
+
+ if (elevation > 270)
+ elevation -= 360;
+ else if (elevation > 90)
+ {
+ elevation = 180 - elevation;
+ azimuth += 180;
+ }
+
+ // azimuth values should be in [0,360[
+ azimuth = fmod(azimuth, 360);
+ if (azimuth < -180)
+ azimuth += 360;
+ else if (azimuth >= 180)
+ azimuth -= 360;
+}
+
+
+// return the distance in degree between 2 normalized angles
+// ie both angles are in [0, 360[
+float HrtfCont::angular_dist(float a1, float a2) const
+{
+ float d = fabs(a1-a2);
+ return d <= 180 ? d : 360 - d;
+}
+
+
+// iangle1, iangle2
+void HrtfCont::add_a2_candidates(interp_cdts& ic, float a1_key, float a2, float weight)
+{
+ // get impulse responses corresponding to the current elevation
+ const angle2_cont& a2_candidates = _m[a1_key];
+ //slog << "interp for fixed angle1=" << a1_key << ", optimal angle2="<< a2<< ", weight=" << weight << endl;
+ // there is only one angle2 measure at index angle1
+ if (a2_candidates.size() == 1)
+ return ic.add(a1_key, a2_candidates.begin()->first, weight);
+
+ // iterator the the 1st element of key bigger or equal
+ angle2_cit supeq_a2it = a2_candidates.lower_bound(a2);
+
+ // the requested angle2 corresponds exactly to an available measure
+ if (supeq_a2it != a2_candidates.end() && a2 == supeq_a2it->first)
+ return ic.add(a1_key, supeq_a2it->first, weight);
+
+ // index candidates for angle2
+ float a2_cand1, a2_cand2;
+ if (supeq_a2it == a2_candidates.end() || supeq_a2it == a2_candidates.begin())
+ {
+ //slog << "angle2 btwn 2 extremes" << (supeq_a2it == a2_candidates.end()) << (supeq_a2it == a2_candidates.begin()) << endl;
+ // requested angle2 is bigger than all available measures
+ // or smaller than all available measures
+ // interpolation between extreme measures can be done
+ // since angle2 is in [-180, 180[
+ a2_cand1 = a2_candidates.begin()->first;
+ supeq_a2it = a2_candidates.end();
+ advance(supeq_a2it, -1);
+ a2_cand2 = supeq_a2it->first;
+ }
+ else
+ {
+ // general case
+ //slog << "angle2 interp general case" << endl;
+ a2_cand1 = supeq_a2it->first;
+ advance(supeq_a2it, -1);
+ a2_cand2 = supeq_a2it->first;
+ }
+ // angular distance between the 2 candidates
+ float a2_cands_dist = angular_dist(a2_cand1, a2_cand2);
+ // weight of the second candidate
+ float a2_cand2_weight = angular_dist(a2, a2_cand1) / a2_cands_dist;
+ //slog << "dist between a2 candidates" << a2_cands_dist << endl;
+ //slog << "angular_dist(a2, a2_cand1)" << angular_dist(a2, a2_cand1) << endl;
+
+ // add 2 candidates
+ ic.add(a1_key, a2_cand1, weight * (1-a2_cand2_weight));
+ ic.add(a1_key, a2_cand2, weight * a2_cand2_weight);
+}
+
+
+void HrtfCont::set_candidates(interp_cdts& ic, float az, float el)
+{
+ // assumption: ic.size == 0, and map size != 0, not checked for RT issues :-(
+
+ // index angles in the hrtf map
+ float iangle1, iangle2;
+
+ //slog << "set candidates az:" << az << ", el:" << el << endl;
+ // get the angular indexes expressed in the coordianates of the HRTF db to use
+ if (_vert_pol_coords)
+ {
+ // Database using vertical polar coordinates
+ normalize_vertpolar_coords(az, el);
+ iangle1 = el;
+ iangle2 = az;
+ }
+ else
+ {
+ // Database using interaural polar coordiantes
+ vertpol2interaurpol(az,el);
+ iangle1 = az;
+ iangle2 = el;
+ }
+ //slog << "indexes " << iangle1 << ", " << iangle2 << endl;
+
+ // there is only one available value for first index
+ if (_m.size() == 1)
+ return add_a2_candidates(ic, _m.begin()->first, iangle2, 1);
+
+ // get the first index value >= iangle1
+ angle1_cit supeq_a1it = _m.lower_bound(iangle1);
+
+ // the requested angle1 corresponds exactly to an indexed element
+ if (supeq_a1it != _m.end() && iangle1 == supeq_a1it->first)
+ return add_a2_candidates(ic, iangle1, iangle2, 1);
+
+ if (supeq_a1it == _m.begin() || supeq_a1it == _m.end())
+ {
+ // requested angle1 is strictly smaller than
+ // the smallest index of the database
+ // or strictly bigger than the biggest index of the db
+ float range_without_measures;
+ float w2;
+ if (supeq_a1it == _m.end())
+ {
+ // requested index1 is bigger than available indexes
+ advance(supeq_a1it, -1);
+ range_without_measures = 2*(90 - supeq_a1it->first);
+ w2 = iangle1 - supeq_a1it->first;
+ }
+ else
+ {
+ // requested index1 is smaller than available indexes
+ range_without_measures = -2*(-90 - supeq_a1it->first);
+ w2 = supeq_a1it->first - iangle1;
+ }
+ w2 /= range_without_measures;
+ add_a2_candidates(ic, supeq_a1it->first, iangle2, 1-w2);
+ iangle2 = iangle2 > 0 ? iangle2 - 180 : iangle2 + 180;
+ add_a2_candidates(ic, supeq_a1it->first, iangle2, w2);
+ }
+ else
+ {
+ // general case: intepolation using 2 different indexed elements
+ float index1_candidate1, index1_candidate2;
+ index1_candidate2 = supeq_a1it->first;
+ advance(supeq_a1it, -1);
+ index1_candidate1 = supeq_a1it->first;
+ // distance between index1 candidates. NB: candidate 2 is bigger than candidate 1!
+ float dcands = index1_candidate2 - index1_candidate1;
+ add_a2_candidates(ic, index1_candidate2, iangle2, (iangle1-index1_candidate1)/dcands);
+ add_a2_candidates(ic, index1_candidate1, iangle2, (index1_candidate2-iangle1)/dcands);
+ }
+}
+
+
+// TODO: could be optimized using templates
+void HrtfCont::update_from_candidates(const interp_cdts& ic, float* left, float* right)
+{
+ //slog << "update from candidates " << left << right << endl;
+ const float *lc, *rc;
+ float wc;
+ ir_buffer* irb;
+
+ irb = &(_m[ic.angle_index1[0]][ic.angle_index2[0]]);
+ //slog << ic.angle_index1[0] << " " << ic.angle_index2[0] << " " << irb->fname << endl;
+ // return;
+
+ wc = ic.weight[0];
+ lc = irb->lbuf;
+ rc = irb->rbuf;
+ //slog << "update from " << _m[ic.angle_index1[0]][ic.angle_index2[0]].fname << " " << wc << endl;
+ //slog << lc << rc << left << right << endl;
+
+ for (size_t i = 0; i < _ir_length; ++i)
+ {
+ left[i] = wc * lc[i];
+ right[i] = wc * rc[i];
+ }
+
+ for (size_t icand = 1; icand < ic.size; ++icand)
+ {
+ irb = &(_m[ic.angle_index1[icand]][ic.angle_index2[icand]]);
+ wc = ic.weight[icand];
+ //slog << "update from " << _m[ic.angle_index1[icand]][ic.angle_index2[icand]].fname << " " << wc << endl;
+ lc = irb->lbuf;
+ rc = irb->rbuf;
+ for (size_t i = 0; i < _ir_length; ++i)
+ {
+ left[i] += wc * lc[i];
+ right[i] += wc * rc[i];
+ }
+ }
+ // slog << "endof candidate update" << endl;
+}
+
+
+// convert an azimuth/elevation couple expressed in
+// vertical polar coordinates to interaural polar coordinates
+void HrtfCont::vertpol2interaurpol(float& az, float& el) const
+{
+ const float raz = az * DEG2RAD;
+ const float rel = el * DEG2RAD;
+ const float cosaz = cos(raz);
+ const float sinaz = sin(raz);
+ const float cosel = cos(rel);
+ const float sinel = sin(rel);
+ const float cosazcosel = cosaz * cosel;
+
+ el = atan2(sinel, cosazcosel) * RAD2DEG;
+ az = atan2(cosel * sinaz, sqrt(cosazcosel * cosazcosel + sinel*sinel)) * RAD2DEG;
+}
generated by cgit v1.2.1 (git 2.18.0) at 2024-05-14 08:06:42 +0000