头部跟踪系统设计在虚拟听觉的应用.docx

上 海 大 学 毕 业 论 文 （设 计）头部跟踪系统设计在虚拟听觉的应用学 生：许仕鹏 学 号：09123819 指导老师：朱梦尧 副教授专 业：电子信息工程 完成年月：2013年06月 目录摘要3abstract4第一章 绪论51.1 课题研究背景及意义51.2 虚拟声音与头部跟踪51.2.1 虚拟声音生成51.2.2 头部跟踪6第二章 头部跟踪系统设计72.1 头部跟踪作用与原理72.2 dcm互补滤波器设计82.2.1 dcm矩阵82.2.2 互补滤波器设计92.3 imu计算方程102.3.1 旋转的表示方式102.3.2 四元数112.3.3 由传感器数值与四元数132.3.4 四元数微分方程152.4 系统实现172.4.1 芯片选择172.4.2 硬件电路192.4.3 软件实现22第三章 3d声音渲染系统253.1 3d声音渲染原理253.2 光线追踪算法263.3 房间冲击响应计算313.3.1 hrtf数据库介绍313.3.2房间冲击响应计算方法333.4 软件实现363.4.1 系统构架363.4.2 多线程声音卷积播放37第四章 整体系统与效果404.1 整体系统404.2 效果说明41英文文献43英文文献中文翻译48总结52致谢53参考文献54附录56头部跟踪核心代码56imu.c563d声音渲染核心代码57hrtf.cpp57wav.c57primitive.cpp58main.cpp59有关图纸62头部跟踪系统原理图6264摘要虚拟听觉技术是继视频技术之后又一重要的人机接口。本课题实现了在虚拟环境中加入与虚拟场景相对应的三维虚拟声音，使收听者在虚拟环境中收听到的声信号逼近在真实环境声场中收听到的声信号，从而增强虚拟现实环境的逼真度和沉浸感。此设计中使用光线追踪的方式模拟声音的传播，通过与头相关函数数据库结合，计算出双耳房间冲击响应，使之与原始声音卷积即可得到与房间真实情况相仿的三维虚拟声音。为了使听者真实感提升，本课题将使用捷联惯导方法实现的头部跟踪系统与虚拟听觉相结合，可以给用户提供更逼真的交互式听觉体验。本设计实现了一个半实时的包含头部跟踪的三维声音渲染系统，效果良好，头部跟踪系统每秒计算姿态400次，每秒上传50帧数据给计算机，精度1-2，三维声音渲染部分可以连续播放声音，房间冲击响应每秒钟更新一次，声音方位感明显。关键词：头部跟踪，三维声音渲染，头相关函数， 双耳房间冲击响应abstractvirtual auditory is another important technology in human machine interface since video technologies. this design combined the 3d virtual sound rendering with a virtual scene, making the listeners hear acoustic signals that approximate to sound in a real scene, thereby enhancing the fidelity and immersive of the virtual reality environment. the system uses ray tracing to simulate sound propagation, with the head-related transfer function(hrtf) database from the mit kemar dummy head model the system generates the binaural room impulse responses(brir) of the virtual scene. after convolving the original sound with the brir, a virtual sound that is similar to the sound in a real environment can be generated. this design also combined the head tracking system that uses the inertial navigation method with the virtual auditory, which can provide users with a more realistic interactive listening experience. this design has implemented a half real time 3d sound rendering system with head tracking. the head tracking system calculates postures 400 times per second, uploads data 50 frames per second through bluetooth. the 3d sound rendering system can render continuously, updates binaural room impulse response once per second.key words: head tracking, 3d sound rendering, head-related transfer function (hrtf), binaural room impulse response (brir).第一章 绪论1.1 课题研究背景及意义虚拟现实环境不仅追求真实的视觉效果，还需有身临其境的听觉效果。虚拟听觉技术是继视频技术之后又一重要的人机接口。虚拟听觉技术就是在虚拟环境中加入与场景对应的三维虚拟声音，使收听者在虚拟环境中收听到的声信号逼近在真实环境声场中收听到的声信号，从而增强虚拟现实环境的逼真度和沉浸感。在移动平台上，音乐的再现往往是通过耳机来实现。而通过耳机的音频再现通常给人的感觉是不自然的，因为声源被定位在头的内部，如果使用3d虚拟声音技术，可以将声源位置移动到外部，给人更好的听觉感受1。心理声学的实验也表面，使用了虚拟声音技术，给以给人接近真实情况的感知。传统的虚拟听觉系统并没有将人头部的运动和姿态考虑进虚拟声音系统，这样限制了声音与人的交互，减少了真实性。虚拟声音中的前后混叠现象十分明显1。现实中，人们可以通过转动头部，感受声音新位置，从而准确判断声音的新位置，减弱前后混叠的影响。所以在虚拟声音中加入头部跟踪系统是必要的。使用了头部跟踪可以大大提高虚拟声源位置的精度。可以使前后混叠发生的平均几率由28%降低到7%1。本课题将头部跟踪系统与虚拟听觉相结合，可以给用户提供更逼真的交互式听觉体验。1.2 虚拟声音与头部跟踪1.2.1 虚拟声音生成虚拟听觉生成的主要算法为：头相关函数(hrtf)和基于光线追踪(ray tracing)的算法。头相关函数用人头录音方式记录了不同方向声源在人耳处的冲击响应，将声音和头相关函数卷积，就可以得到来自该方向的声音。基于光线追踪算法的虚拟听觉采用光线追踪来模拟声音的传播，通过计算可以得到声音从声源-经过场景的反射、衍射、折射-到听者的房间brir响应。由于声场的波长比场景中障碍物的尺寸大，所以使用这样的方法计算散射是合适的2。图1-1 hrtf数据以及计算出的 brir1.2.2 头部跟踪头部跟踪系统的应用十分广泛，它可以用来捕捉人头部的运动，近年来在游戏、虚拟现实领域发挥着越来越重要的角色。其实现主要有基于光学和基于捷联惯导算法两种方式。基于光学的方式比起捷联惯导方式的优势在于可以检测到头部的平动，而捷联惯导方式只可以检测到头部转动。不过光学方式需要高速红外摄像头，实现成本很大，同时实时性并不好。鉴于本系统用于虚拟听觉，对转动要求高，平动部分可以使用键盘控制，所以我们将使用捷联惯导算法来实现头部跟踪系统。为了对传感器数据滤波，本系统中使用的是dcm互补滤波3方法，虽然kalman滤波20方法可能会有更好的效果，不过考虑到在单片机中计算的时间复杂度，而且两者实际效果差距相似，所以我们使用互补滤波方法。图1-2 头部跟踪系统 人类听觉在水平面上可以对声源做非常准确的定位，因为水平面上不单单是依靠单耳来确定位置，而是依靠双耳的时间差itd、强度差ild来定位。而确定仰角时只能依靠单耳信息，也就是人头部、身体、外耳廓对声音频谱的滤波作用，所以是不准确的1。由于人类对水平方位感知更明显，在头部跟踪系统中我们着重考虑水平方向的精度。在水平面上，人类对声音角度感知的最小值是1左右，少于1的角度变换人类无法察觉。本系统的精度在水平方向是1-2，所以是符合虚拟声音对头部跟踪的要求的。第二章 头部跟踪系统设计2.1 头部跟踪作用与原理在虚拟现实系统中，如果可以知道用户头部的转动姿态，根据这些信息调整虚拟现实中场景声音的生产，使用户的获得更真实的感受。图2-1 头部跟踪在虚拟听觉中的作用如果没有头部跟踪，用户头部旋转了90度，虚拟声源的位置也会随之旋转90度；而如果加上头部跟踪系统，转动时可以将声源位置反向旋转同样角度，这样，声音对于用户依然保持在同一位置。图2-2 头部的转动头部跟踪系统的本质是通过传感器来检查头部在x, y, z三轴方向的转动角 roll ，pitch , yaw。使用加速度传感器可以检测到重力的方向，有了重力的方向，我们就可以用之计算出roll和pitch信息；使用地磁传感器可以检测出地磁场的方向，使用地磁信息，我们可以计算出yaw方向的数值。不过由于无论是重力计还是地磁计，它的响应速度都是较慢的，为了纠正他们的数据，还需要使用陀螺仪来检测角速度。使用角速度信息来加快姿态调整。传感器数据波动较大，如何更好地融合加速度计、地磁计和陀螺仪的数据计算出准确的姿态是亟待解决的问题，这就需要dcm互补滤波器。2.2 dcm互补滤波器设计传感器输入数据不稳定，为了得到稳定的结果，需要融合不同的传感器数据进行滤波。2.2.1 dcm矩阵首先，我们需要定义2个坐标系：地理坐标系：真实世界所在的坐标系。重力向量为（0,0,1g），地磁向量是（0,1,x）两个三维向量,对于地磁我们不关心其垂直分量。头部坐标系：上面有加计、陀螺、电子罗盘传感器，三个三维向量。图2-3 两个坐标系的转换方向余弦矩阵(direction cosine matrix)是表示从坐标系oxyz到坐标系oxyz变换的矩阵，使用它，我们可以实现地理坐标系与头部坐标系这两个坐标系间的转换。如上图所示，这是一个3维空间中的坐标系变化。在我们的系统中，oxyz是世界坐标系，oxyz是头部跟踪系统所在的局部坐标系。使用dcm的目的也就是将坐标从世界坐标系变换到头部坐标系，以便进行姿态计算。设x, y, z为世界坐标系中的3个单位向量：x=100 y=010 z=001 设x, y, z是局部坐标系中的3个单位向量：x=100 y=010 z=001在全局坐标中：xg=xxgxygxzg则：xxg=xgcos(x,xg)其中cos(x,xg)是x和xg的夹角。所以有：即为dcm矩阵，表示从世界坐标系到头部局部坐标系的变换。且有mg=dcmgmb, mb=dcmbmg。2.2.2 互补滤波器设计图2-4 互补滤波示意图为了准确融合加速度、地磁、角速度三个传感器的结果，我们需要使用互补滤波器。这样做有2个好处：1、高通滤波器可以抑制陀螺仪低频漂移2、低通滤波器可以很好的抑制加速度计的噪声也就是说，加计和地磁传感器在长久来看是很准确的，而陀螺仪在短时间内是很准确的，结合二者数据，可以得到较优的结果。伪代码：angle = (k)*(angle + gyro * dt) + (1-k)*(x_acc);上述伪代码即为互补滤波器融合加计和陀螺仪数值的简单例子。2.3 imu计算方程imu(inertial measurement unit, 惯性测量单元)是用于姿态解算的单元。姿态解算运算中需要用到大量的旋转操作，所以选择合适的旋转方式对系统的计算十分重要。2.3.1 旋转的表示方式常用的选择表示方法有：旋转矩阵、欧拉角、四元数。旋转矩阵、欧拉角、四元数主要用于:向量的旋转、坐标系之间的转换、方位的平滑插值计算。表2-1 各旋转表示方法比较任务/性质旋转矩阵欧拉角四元数在坐标系间(物体和惯性)旋转点能不能(必须转换到矩阵)不能(必须转换到矩阵)连接或增量旋转能,但经常比四元数慢,小心矩阵蠕变的情况不能能,比矩阵快插值基本上不能能,但可能遭遇万向锁或其他问题slerp提供了平滑插值易用程度难易难在内存或文件中存储9个数3个数4个数对给定方位的表达方式是否唯一是不是,对同一方位有无数多种方法不是,有两种方法,它们互相为共轭可能导致非法矩阵蠕变任意三个数都能构成合法的欧拉角可能会出现误差积累,从而产生非法的四元数欧拉角显然是最直接的方式，与人类对旋转的感知最为接近，但是欧拉角会出现万向锁这样的问题，即三轴系统中由于其中一轴旋转90度导致损失一轴的自由度，三轴系统的自由度退化为2轴。也就是说：当任意一轴旋转90的时候会导致该轴同其他轴重合，此时旋转被重合的轴无法产生任何效果，这就是万向锁(gimbal lock)。图2-5 万向锁作为旋转变换，矩阵是必要的，可以把四元数或欧拉角作为主数据，同时维护一个旋转矩阵，当四元数或欧拉角变化时更新矩阵。欧拉角虽然占用更少的内存(比四元数少1/4)，但是转换到矩阵的时间要比四元数更慢一些。平滑的插值只能用四元数完成。直接用欧拉角插值会产生奇怪的转动。综上考虑，我们选择四元数作为旋转的表示方式，因为四元数插值方便、没有万向锁问题、运算速度较快。2.3.2 四元数四元数是由哈密顿在1843年爱尔兰发现的。当时他正研究扩展复数到更高的维次（复数可视为平面上的点）。他不能做到三维空间的例子，但四维则造出四元数。根据哈密顿记述，他于10月16日跟他的妻子在都柏林的皇家运河（royal canal）上散步时突然想到i2=k2=j2=ijk=-1的方程解。之后哈密顿立刻将此方程刻在附近布鲁穆桥（brougham bridge，现称为金雀花桥 broom bridge）。这条方程放弃了交换律，是当时一个极端的想法（那时还未发展出向量和矩阵）4。四元数常见于控制论、信号处理、姿态控制，用来表示旋转和方位。这是由于四元数转换组合比很多矩阵转换组合在数字上更稳定。四元数一般定义如下：q=w+xi+yj+zk其中 w,x,y,z是实数。同时，有:i2=j2=k2=-1四元数也可以表示为：q=w,v其中v=(x,y,z)是矢量，w是标量，虽然v是矢量。当用一个四元数乘以一个向量时，实际上就是让该向量围绕着这个四元数所描述的旋转轴，转动这个四元数所描述的角度而得到的向量。四元数加法：四元数取模：四元数乘法：四元数共轭：四元数倒数：对于四元数，我们最重要的是用是用它来旋转一个向量，假设有一个用四元数表示的旋转为：，有一个向量 v=(vx,vy,vz)，则经过旋转后的结果v为：v=qvq*其中 q*为四元数q的共轭四元数，是转动角度。有了四元数，如果应用于变换，我们还需要将四元数转化成相对应的矩阵。把上述的qvq*表示为qv，则有上式在其次坐标下可以简化为：当n(q)=1时，上式可以进一步简化为：2.3.3 由传感器数值与四元数我们可以获取到地磁、加速度计和陀螺仪的数值，根据这些数值我们可以实时更新四元数方程，使四元数表示的方位与头部所在的姿态保持一致。首先，我们需要根据之前计算出的四元数来反算出前一时刻头部所在的姿态的重力方向，然后，计算此时传感器得到的重力信息和前一时刻的差，更新四元数。vx=2*(q1*q3-q0*q2)vy=2*(q0*q1+q2*q3) vz=q0*q0-q1*q1-q2*q2+q3*q3上述公式即为由前一状态四元数计算出重力参考方向。vx, vy, vz分别为重力在三轴的分量。这个公式来自从2.3.1中的：重力向量就是q的第三行。与重力一样，我们也需要计算地磁场的参考方向，并且依据地磁的方向来更新yaw方向的四元数。地磁参考方向计算如下：hx=2*mx*(0.5-q2q2-q3q3)+2*my*(q1q2-q0q3) + 2*mz*(q1q3+ q0q2); hy=2*mx*(q1q2+q0q3)+2*my*(0.5-q1q1-q3q3)+2*mz*(q2q3-q0q1); hz=2*mx*(q1q3-q0q2)+2*my*(q2q3+q0q1)+2*mz*(0.5-q1q1-q2q2); bx = sqrt(hx*hx) + (hy*hy);bz = hz; mx,my,mz是地磁传感器得到的数值的归一化表示，hx,hy,hz是计算出来的地磁参考方向。由于我们并不关心地磁垂直方向的量，所以把地磁的y轴和x轴做了融合，得到新的参考方向bx,bz。之后是由上一时刻四元数表示的方向计算地磁的方向：wx = 2*bx*(0.5 - q2q2 - q3q3) + 2*bz*(q1q3 - q0q2); wy = 2*bx*(q1q2 - q0q3) + 2*bz*(q0q1 + q2q3);wz = 2*bx*(q0q2 + q1q3) + 2*bz*(0.5 - q1q1 - q2q2); 上公式依旧来源于2.3.1中的有了重力和地磁的参考方向，我们可以计算出传感器得到的方向与场的方向间的差值，这种求差值本质上是进行了一次向量叉乘运算：e=av+mw其中，e是传感器和场方向的差，a是加速度计数值，v是重力场参考方向，m是地磁计数值，w是地磁场参考方向。将向量展开成分量，则有： ex = (ay*vz - az*vy)+(my*wz - mz*wy); ey = (az*vx - ax*vz)+(mz*wx - mx*wz); ez = (ax*vy - ay*vx)+(mx*wy - my*wx);最后，我们需要把陀螺仪的数值融合到误差e中，对e做一些pid控制，使结果更加平滑准确。 gx = gx + imu_kp*ex + exint; gy = gy + imu_kp*ey + eyint; gz = gz + imu_kp*ez + ezint;上述伪代码中的gx,gy,gz即是陀螺仪的三轴的数值，与上面的e相融合以后，就可以使用在四元数微分方程里来更新四元数了。2.3.4 四元数微分方程数值分析中，龙格库塔法（runge-kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔龙格和马丁威尔海姆库塔于1900年左右发明5。下面举一个具体例子来说明龙格库塔方法。令初值问题表述如：则，对于该问题的如果使用四阶龙格库塔方法由如下方程给出：其中下一时刻的输出yn+1是由当前值加上估算的斜率与时间差的乘积。估算的斜率是k1,k2,k3,k4斜率的加权平均，中点的斜率有更大的权值：上例中使用的是四阶方法，考虑到单片机的处理速度，我们使用的是1阶龙格库塔法。设有一个矩阵微分方程：x(t)=fxt,w(t)其一阶龙格库塔的解为：xt+t=xt+tfxt,w(t)对于四元数q，其微分方程为：q=12q*b矩阵形式为：则带入一阶龙格库塔的解为 qt+t=qt+t*12b*q(t)即：表示为伪代码为： iq0 = (-q1*gx - q2*gy - q3*gz)*halft; iq1 = (q0*gx + q2*gz - q3*gy)*halft; iq2 = (q0*gy - q1*gz + q3*gx)*halft; iq3 = (q0*gz + q1*gy - q2*gx)*halft; 其中，iq0,iq1,iq2,iq3分别为q的四个分量的增量，halft是上面公式中的12t。qnext=q+iq2.4 系统实现2.4.1 芯片选择实现头部跟踪设备时，我们综合考虑单片机的性能、价位、资料丰富程度等因数，选择了stm32系列的stm32f103c8t6单片机作为系统的处理器。stm32是基于arm cortex-m3内核的32位微处理器，最高到72mhz的主频，可以实现单周期乘法和硬件加速的除法，它以8位机的价格实现了32位机的性能。同时资料丰富，爱好者众多。系统中主要的传感器芯片主要有2片：mpu6050：3轴加速度计与3轴陀螺仪hmc5883l：3轴电子罗盘系统中还使用了蓝牙模块hc-06与pc通信。图2-6 mpu6050mpu-6000为全球首例整合性6轴运动处理组件，相较于多组件方案，免除了组合陀螺仪与加速器时之轴间差的问题，减少了大量的封装空间。mpu-6050整合了3轴陀螺仪、3轴加速器，并含可藉由第二个i2c端口连接其他厂牌之加速器、磁力传感器、或其他传感器的数位运动处理(dmp: digital motion processor)硬件加速引擎，由主要i2c端口以单一数据流的形式，向应用端输出完整的9轴融合演算技术。mpu-6050的角速度全格感测范围为250、500、1000与2000/sec (dps)，可准确追緃快速与慢速动作，并且，用户可程式控制的加速器全格感测范围为2g、4g8g与16g。产品传输可透过最高至400khz的i2c或最高达20mhz的spi。mpu6050大量使用在姿态检测、体感游戏与手机中。图2-7 hmc5883lhoneywell hmc5883l 是一种表面贴装的高集成模块，并带有数字接口的弱磁传感器芯片，应用于低成本罗盘和磁场检测领域。hmc5883l 包括最先进的高分辨率 hmc118x 系列磁阻传感器，并附带honeywell专利的集成电路包括放大器、自动消磁驱动器、偏差校准、能使罗盘精度控制在 12的 12 位模数转换器.简易的 i2c 系列总线接口。hmc5883l 是采用无铅表面封装技术，带有 16 引脚，尺寸为 3.0x3.0x0.9mm。hmc5883l 的所应用领域有手机、笔记本电脑、消费类电子、汽车导航系统和个人导航系统。 hmc5883l 用honeywell各向异性磁阻(amr)技术，该技术的优点是其他磁传感器技术所无法企及。这些各向异性传感器具有在轴向高灵敏度和线性高精度的特点.传感器带有的对于正交轴低敏感行的固相结构能用于测量地球磁场的方向和大小，其测量范围从毫高斯到 8 高斯(gauss)。图2-8 hc-06同时，为了让系统可以无线的与pc通信，发回姿态数据，我们选择hc-06蓝牙串口芯片作为通信芯片。hc-06有以下特点：蓝牙 2.0 带 edr， 2mbps-3mbps 调制度内置 2.4ghz 天线, 用户无需调试天线外置 8mbit flash低电压 3.3v 工作（3.1v4.2v）配对时3040ma波动，配对完毕通信 8ma可选 pio 控制标准 hci端口(uart or usb)作为一个完整的蓝牙模块，hc06给用户使用蓝牙提供了极大地便利。在使用hc-06芯片前，需要使用串口的tx和rx接入hc-06的tx和rx来配置hc-06的蓝牙配对信息。hc-06的默认波特率是9600，首先将电脑串口也设为9600波特率，然后发送”at”，如果hc-06返回”ok”，则代表hc-06工作正常。本系统使用的是115200，所以需要改变系统的波特率。改变波特率使用：at+baudn，n与波特率的对应关系是：1:1200;2:2400;3:4800;4:9600;5:19200;6:38400;7:57600;8:115200;9:230400;a:460800;b:921600;c:1382400所以本系统中需要执行”at+baud8”，如果设置正常，hc-06会返回”ok115200”，此后，hc-06的串口波特率已经改为115200，所以，pc上的波特率也需要做响应的修改。不建议用在超过115200的波特率，信号的干扰会使系统不稳定。设置超过115200后用电脑无法使用，要用单片机编程于高于115200才能使用此波特率。用at命令设好波特率后，下次上电使用不需再设，可以掉电保存波特率。发送”at+namexxx”可以设置搜索蓝牙时的蓝牙显示名。发送”at+pin1234”可以设置于其他蓝牙设备配对时的配对密码，出厂默认密码是1234。2.4.2 硬件电路图2-9 头部跟踪系统pcb预览设计头部跟踪系统时候，原理图使用cadence orcad绘制，pcb使用cadence allegro绘制。上图为pcb的3d预览和布局图。制版选择在深圳的pcb工厂，他们提供高质量快速的小批量pcb打样业务。pcb第一次打板时上有电路问题，第二次重新调整电路后加工得到现在使用的头部跟踪系统。电源部分：图2-10 电源电路电源电路部分主要有ams1117作为稳压芯片，保证系统的vcc维持在3.3v。ams1117边的10uf电容用于滤去电源低频噪声, 0.1uf电容用于滤去电源高频噪声。d4是一个发光二极管，提示系统是否上电。单片机电路：图2-11 单片机电路单片机电路中，系统使用8m外部晶振，内部通过pll倍频到72mhz；boot0和boot1都用10k电阻下拉到gnd，boot0和boot1这两个管脚在芯片复位时的电平状态决定了芯片复位后从哪个区域开始执行程序：boot1=x boot0=0 从用户闪存启动，这是正常的工作模式。boot1=0 boot0=1 从系统存储器启动，这种模式启动的程序功能由厂家设置。boot1=1 boot0=1 从内置sram启动，这种模式可以用于调试。蓝牙电路：图2-12 蓝牙电路hc-06模块将蓝牙封装成串口，只需要将单片机串口的tx和rx接过来，就可以使用蓝牙通信了。led d6的作用是提示蓝牙状态，蓝牙未配对时d6会不停闪烁，配对成功时d6会保持常亮状态。惯性测量模块电路：图2-13 惯性导航电路惯性测量部分主要是mpu6050和hmc5883l，它们的scl和sda上拉了10k的电阻接到stm32的pb6和pb7，通过这两个端口来传输数据。 图2-14 实物图 左图为背面 右图为正面2.4.3 软件实现头部跟踪软件部分使用keil 4开发，用j-link v8下载和调试程序。代码结构：srcalgorithm：imu算法srcalgorithmimu.c(.h)：惯性导航姿态检测实现srcdriver：驱动程序srcdriverdelay.c(h)：基于system tick的延迟srcdriverhmc5883l.c(h)：hmc5883l地磁传感器读取srcdriverioi2c.c(h)：gpio模拟i2c驱动srcdrivermpu6050.c(h)：mpu6050读取srcdrivertimer.c(h)：32位高精度计时器，用于算法时间差计算srcdriveruarts.c(h)：串口通信srcmain.c：主程序软件部分总体的框架是： while(1) 读取mpu6050 读取hmc5883l 姿态计算 imu_getyawpitchroll 蓝牙串口发送帧头标记位 蓝牙串口发送加计数值 ax ay az 蓝牙串口发送陀螺仪数值gx gy gz 蓝牙串口发送地磁计数值mx my mz 蓝牙串口发送解算的姿态yaw pitch roll 其中最核心的是姿态解算部分void imu_getyawpitchroll(int16_t * angles,int16_t *data)，*angles是指向解算结果的指针(angles0是yaw, angles1是pitch, angles2是roll)，*data是指向传感器数据的指针(data0, data1,data2是3轴加计数值, data3, data4,data5是3轴陀螺仪数值, data6, data7,data8是3轴地磁计数值)。在imu_getyawpitchroll中，首先对加速度计数据ax,ay,az和地磁计数值mx,my,mz进行归一化： norm = invsqrt(ax*ax + ay*ay + az*az); ax = ax * norm; ay = ay * norm; az = az * norm; norm = invsqrt(mx*mx + my*my + mz*mz); mx = mx * norm; my = my * norm; mz = mz * norm;然后，计算地磁的参考方向与地磁方向，具体解释请参考上一节： hx = 2*mx*(0.5 - q2q2 - q3q3) + 2*my*(q1q2 - q0q3) + 2*mz*(q1q3 + q0q2); hy = 2*mx*(q1q2 + q0q3) + 2*my*(0.5 - q1q1 - q3q3) + 2*mz*(q2q3 - q0q1); hz = 2*mx*(q1q3 - q0q2) + 2*my*(q2q3 + q0q1) + 2*mz*(0.5 - q1q1 - q2q2); bx = sqrt(hx*hx) + (hy*hy); bz = hz; wx = 2*bx*(0.5 - q2q2 - q3q3) + 2*bz*(q1q3 - q0q2); wy = 2*bx*(q1q2 - q0q3) + 2*bz*(q0q1 + q2q3); wz = 2*bx*(q0q2 + q1q3) + 2*bz*(0.5 - q1q1 - q2q2);之后，计算重力的参考方向： vx = 2*(q1q3 - q0q2); vy = 2*(q0q1 + q2q3); vz = q0q0 - q1q1 - q2q2 + q3q3;之后，计算场与传感器方向的差值 ex = (ay*vz - az*vy) + (my*wz - mz*wy); ey = (az*vx - ax*vz) + (mz*wx - mx*wz); ez = (ax*vy - ay*vx) + (mx*wy - my*wx);之后，使用pid更新差值，并把差值融合到gx,gy,gz中： exint = exint + ex*ki * halft; eyint = eyint + ey*ki * halft; ezint = ezint + ez*ki * halft; gx = gx + kp*ex + exint; gy = gy + kp*ey + eyint; gz = gz + kp*ez + ezint;之后，使用四元数微分方程(龙格库塔方法)来更新四元数，并重新标准化四元数： iq0 = (-q1*gx - q2*gy - q3*gz)*halft; iq1 = (q0*gx + q2*gz - q3*gy)*halft; iq2 = (q0*gy - q1*gz + q3*gx)*halft; iq3 = (q0*gz + q1*gy - q2*gx)*halft; q0 += iq0; q1 += iq1; q2 += iq2; q3 += iq3; norm = invsqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); q0 = q0 * norm; q1 = q1 * norm; q2 = q2 * norm; q3 = q3 * norm;最后，由四元数结果反算出欧拉角数值： angles0 = (int16_t)(atan2(2 * q1 * q2 - 2 * q0 * q3, 2 * q0*q0 + 2 * q1 * q1 - 1) * 1800/m_pi); angles1 = (int16_t)(atan(agx / sqrt(agy*agy + agz*agz) * 1800/m_pi); angles2 = (int16_t)(atan(agy / sqrt(agx*agx + agz*agz) * 1800/m_pi);以上为单片机端软件的主体内容。同时，为了验证姿态解算结果正确，在pc上，使用visual studio c#开发了一个可视化的查看器，通过蓝牙串口收取姿态解算数据，并使用木盒代替人头将姿态显示出来。3d姿态显示部分使用的是opengl。图2-15 ahrs viewer截图ahrs viewer (attitude heading reference system, 航向姿态参考系统) 中的3d view中会显示头部的实时姿态，三个进度条也会显示出yaw, pitch, roll三个方向上的角度大小。左上角的sensor result中，实时显示了各个传感器的数值。左下角的”fps”标签记录了当前秒一共收到了多少帧的数据。如果串口是在打开软件之后插入电脑的话，需要先点击”refresh ports”按钮，否则”serial port”列表中不会显示新加的串口。第三章 3d声音渲染系统3.1 3d声音渲染原理3d声音技术的本质是模拟声音从声源到听者位置所产生的振动。双耳间会有时间差（itd，interaural time difference）和强度差（ild，interaural level difference），同时声音在传入人耳的过程中，会被人的身体、头部、耳廓、中耳等吸收，从而改变声音的性质。人脑可以通过这些区别来感知声源的位置、方向。我们可以使用头相关传递函数（hrtf，head-related transfer function）来记录声音从声源到听者的这种变化，本质上，hrtf是从声源到听者的冲击响应。有了来自某方向的hrtf数据使用它与原始声音相卷积就可以得到来自该方向的声音。即：listener=hrtf*source图3-1 人头录音示例hrtf的获取使用的是人头录音方法。首先，在假人头中的左右耳道中放置麦克风，用来记录左右耳听到的声音。之后，围绕着假人头用扬声器播放冲击，与此同时，耳道里的双麦克风记录冲击的响应，这个响应便是这个方向上从声源到听者的hrtf。光有hrtf是不足够的，因为hrtf只是记录了声音的方位信息，对于一个具体的场景，声音不仅有方向，声音还会在场景中反射、折射，声源离听者的距离也会不固定。于是，我们提出使用光线追踪来模拟声音的传播，并模拟声线与听者的碰撞，碰撞时记录每个声线此时的传播时间、强度和方位，最终使用这些信息，生成双耳房间冲击响应brir。通过brir与原始声音的卷积，我们可以模拟在房间里，从声源到听者发出的声音。这样方式产生的声音可以给用户更加的沉浸感，3d声音听起来也越真实。3.2 光线追踪算法光线追踪（ray tracing），又称为光线跟踪，来自于几何光学的一项通用技术，它通过跟踪与光学表面发生交互作用的光线从而得到光线经过路径的模型6。光线追踪广泛用于计算机图形学中，使用这样的方法，可以比光栅化图像生成达到更加好的光学的效果。同时，光线追踪对于折射和反射都有非常好的效果，所以在对图像质量要求高的场合，例如电源渲染中，常常会采用这种方法。图3-2 光线追踪示例在实际应用中，可以将各种电磁波或者微小粒子看成光线，我们可以利用光线追迹来计算微小粒子在介质中传播的情况。首先发射一条光线到空间中，当光线碰撞到障碍物时，原光线传播方向改变，沿着反射角方向传播。声音是波，不过考虑到声波的波长远小于房间中墙壁的尺寸，我们认为使用光线追踪这样粒子性较强的算法是合适的。每个声线在房间中传播，遇到墙壁会反射，在空气中传播会有损耗，最终会撞上听者的声音是我们所需要关注的，使用这些声音的方位、时延、强度，我们可以获得brir。要实现光线追踪，首先需要定义几个类：primitive：场景的基元，本质上是一个三角形，类中有检测碰撞、求交、渲染三角形等功能。ray：声线，类中记录声线的起点、方向、强度、传播距离等信息。scene：场景类，是primitive的集合，包含读取obj模型、整个场景渲染等功能。vector：向量类plane：平面类，包含法相n和距离d光线追踪伪代码：对于每一个基元如果光线与该基元碰撞光线强度根据光线起点到碰撞位置的距离衰减将原光线的设为反射的方向光线强度再减去碰撞的衰减如果光线的强度为0，将光线设为不活跃状态重复上面过程直到活跃光线数目为0在primitive类中，主要的计算是：（1）根据三个顶点计算平面方程；（2）光线与三角形求交。根据三个顶点计算平面方程：设有三个点p1,p2,p3，p12=p1-p2, p23=p2-p3,则平面的法相量为normal=p12p23平面方程的d参数为：d=normalp3代码示意： vector3 p12,p23; p12=p1-p2; p23=p2-p3; vector3 normal = p12.cross(p23); normal.normalize(); float d = dot( normal, p3 ); m_plane.n=normal; m_plane.d=d;光线与三角形求交分为2步1.光线与三角形所在的平面求交，计算距离2.看交点是否在三角形内部对于1，设有光线p=p0+tv,平面pn+d=0，约掉p，则有p0+tvn+d=0其解为：t=-(p0n+d)/(vn)示例代码：float dist = dot(p1-ray.getorigin(),m_plane.n)/(dot(ray.getdirection(),m_plane.n);图3-3 三角形与光线碰撞示例之后检测交点是否在三角形内部，对于三角形的每一条边，v1=t1-pv2=t2-pn1=v1v2normalize(n1)d=-p0n1如果对每条边pn1+d10，则交点在三角形内部。示例代码： if (dist 0) if (dist = 0) return true; else return false; scene类主要存储了一个primitive列表std:vector primlist;。当声线ray做碰撞检测时都是访问scene里的intersect接口进行碰撞检测，scene里的intersect负责返回被碰撞的primitive的id与从声线起点到该primitive的距离。同时scene负责prim