From 45dfb5debd4a0d0bc79ad8b6691576e8462d1396 Mon Sep 17 00:00:00 2001 From: David Doukhan Date: Mon, 7 Feb 2011 13:09:27 +0000 Subject: cw_binaural~ code and examples! Makefile still missing svn path=/trunk/externals/ddoukhan/cw_binaural~/; revision=14855 --- src/hrtfcont.cpp | 256 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 256 insertions(+) create mode 100644 src/hrtfcont.cpp (limited to 'src/hrtfcont.cpp') 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 . +*/ + + +#include +#include +#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; +} -- cgit v1.2.1