2020年认证杯SPSSPRO杯数学建模
A题 听音辨位
原题再现：
把若干 (⩾ 1) 支同样型号的麦克风固定安装在一个刚性的枝形架子上 (架子下面带万向轮，在平地上可以被水平推动或旋转，但不会歪斜)，这样的设备称为一个麦克风树。不同的麦克风由于位置不同，录制到的声音往往也有细微的不同，所以通过对多支麦克风接收到的声音进行对比分析，可以得到更多的有关声源的信息。我们假设每个麦克风都是全向的，也就是单麦克风无法分辨声源的方向。
现在有一个地面、墙壁和天花板都是光滑大理石的大厅，大厅内空旷而安静。在大厅里只有一个走动的人，发出清晰的脚步声。我们准备在大厅里安放一个麦克风树，希望通过检测声音来进行一些测量，包括尽量准确地实时确定这个人的位置，也包括测量这个大厅的某些几何参数。请你建立合理的数学模型，设计一个成本尽量低、而且可以达到使用要求的麦克风树。要求给出每支麦克风的相对位置以及相对于地面的高度，至于枝形架子的具体力学结构则不需要考虑。
第二阶段问题： 在录制唱片时，有一种非常特殊的录音方式，称为人头录音。人头录音通常是使用一个人头的塑胶模型，将两个麦克风分别放在两个耳朵里，分别负责左右两个声道的录制，以此来仿效人耳的收听效果。人头录音的唱片在使用耳机收听时有极佳的定位效果，可以使收听者相当清楚地感受到声源在收听者的前后左右等不同位置，有如身临其境。当然它只能通过耳机来回放，使用音箱就失去了应有的效果。
已知大厅的平面形状是矩形，地面和天花板都是水平的。此时我们不知道大厅的轮廓尺寸。大厅中可能有少量几根方形柱子，柱子的侧壁和墙壁平行，但具体位置和尺寸未知。柱子表面也是光滑的大理石。
我们希望使用一个麦克风树来录制大厅里的音效，并将多个麦克风的录制信号进行后期处理，以尽量准确地模仿人头录音的效果。请给出成本尽量低的麦克风树的设计方案，以及对录制信号的后期处理算法。
整体求解过程概述(摘要)
麦克风树在听声辨位研究中逐渐得到了普及，许多兴趣人士也乐此不疲的在不同方向上展开研究探索，如对通过麦克风树来模仿人头录音效果问题的探究。本文针对麦克风树听声辨位模仿人头录音效果的问题，建立了方案寻优模型得出了以 5 支麦克风组成的十字形均排布的设计方案，同时针对大厅中存在方形柱子的情况给出了声源定位算法。最后，使用后期音频处理算法对录制声音信号进行三维音效渲染，生成了具有人头录音效果的双耳三维音频信号。
针对问题一，考虑到成本尽量低和尽可能准确地模仿人头录音的效果，建立了麦克风树的最优化模型。人头录音唱片在收听者使用耳机收听时，收听者会感受到声源在自身的前后左右不同位置，为满足这一定位效果，在成本尽可能低的约束下，对麦克风阵列在数量上和结构上进行分析研究，分别对 2 支，3 支，4 支和 5 支麦克风以及一维、二维排布的定位效果进行寻优。最终，给出的设计方案为：由 5 支麦克风组成二维水平十字型麦克风树最优化模型。
针对问题二，基于点声源传播特点和空间几何方法，给出了时实的声源点定位算法。基于声音在传播过程中为直达波的条件下，对问题一中的最优设计方案进行分析。由于每 2 支麦克风的接收途径可得两个测量坐标，录制信号不受方形柱子影响时，5 支麦克风组成的阵列最多会产生 20 种不同的声源位置。其中有 10 个坐标的位置是比较接近的，另外有 10 个坐标的位置是分散的，对 10 个坐标的位置是比较接近的求加权平均值，最终较为精确地定位出声源位置的具体坐标。录制信号受到方形柱子影响时，由于大厅柱子的影响，可能导致一至两个麦克风接收到的信号数据失真，将这些麦克风的失效数据全部舍弃，用剩余的几个麦克风来确定声源的位置坐标。
针对问题三，建立了录制信号的后期处理算法，并对音频信号进行渲染处理，进而得出了较准确的人头录音音效的音频信号。首先，基于头相关传输函数 HRTF 算法，结合录制得到的音频信号，针对声源信号的坐标位置进行三维立体音效渲染；然后，利用双耳三维音频合成算法合成具有空间方位感的三维音频信号；最后，将渲染后的三维音频信号通过双耳耳机播放，即可获得人头录音音效。
最后，使用 Matlab 编程进行仿真模拟实验，利用仿真模拟实验的数据对声源定位算法的精准度进行定量分析与误差检验，并对音频信号处理算法的误差进行了定性分析。综合分析得出：本模型具有较高的可靠性、准确性和精准度，且算法处理后的音频信号能够较准确的模拟人头录音效果。
问题分析：
问题的总体分析
本文的整体任务：给出成本尽量低、能够尽可能达到人头录音效果的麦克风设计方案，并给出声音信号后期处理算法。
由任务出发进行分析：考虑到人头录音效果有极佳的定位效果，可以使用头相关传输函数针对较为精准的声源空间位置进行三维音效渲染，进而使音效能够达到人头录音的效果。于是，问题转化成求解较为精准的声源定位问题，需给出声源定位算法，进而以声源定位的精度与成本作为约束条件求解最优化设计方案。故本文需解决的问题可以分为问题一、问题二、问题三去求解，即问题一为设计最优化麦克风数方案的求解，问题二为声源定位算法的求解，问题三后期声音信号处理算法的求解。
综上可知：本文通过对麦克风树听声辨位并达到人头录音效果的问题进行分析研究，首先需建立满足定位要求和人头录音效果的麦克风树优化模型，并通过所建立的麦克风树模型给出大厅中存在柱子影响时人的定位，最后，为满足更加贴近人头录音效果，即收听者在佩戴耳机时可感受到声源在自身前后左右等方位的情境，给出相应的录制信号的后期处理算法。
问题一的分析
本题麦克风树在满足成本尽可能低的约束条件下，在原有麦克风定位声源的基础上，要尽可能准确地模仿人头录音的效果，人头录音的唱片在收听者使用耳机收听时，可感受到声源在收听者的前后左右等位置。对于问题中麦克风树的设计方案可解读为在满足要求的前提下建立最优化模型。
问题二的分析
本题主要是对上述寻优模型中得出的优化方案的定位求解，对优化方案进行分析后，我们可以在问题二中建立模型，当录制信号不受柱子影响时，通过对 5 只麦克风录制声源数据的处理，可以求得更加贴近于真实位置的坐标点。当录制信号受柱子影响时，导致一至两个麦克风所接收到的数据失效，用剩余的麦克风接收到的数据来定位声源位置坐标
问题三的分析
问题三的目的：建立后期录制的音频信号处理算法，使处理后的音频信号具有较准确人头录音的音效。
思路分析：首先，考虑到头相关传输函数能够较好地渲染三维音效，且能够获得较好的人头录音效果；然后，给出头相关传输函数相关求解算法，并使用头相关传输函数针对声源空间位置进行三维音效渲染，再通过双耳三维音频合成算法合成渲染后的音频信号；最后，通过双耳耳机输出处理后的信号，即可得到具有人头录音音效的信号。
模型假设：
假设一：每支麦克风都是全向的；
假设二：录制的声音回放只通过耳机进行；
假设三：声音在空气中传播所衰减的能量可以忽略不计；
假设四：声音在大理石表面反射时会产生一定能量的衰减；
论文缩略图：
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
%本程序用于仿真模拟声源发出声音信号的实验
clear;clc;
H=1.8;%麦克风树的高度
R=1;%麦克风距离十字中心点的距离
%a=45;b=20;c=4;%墙长 45，宽 20，高 4
v=340;%取声速 340
n=[0,0,H;0,-R,H;R,0,H;0,R,H;-R,0,H];%麦克风的坐标位置
x1=-6:3:6;%取样本点来模拟
x2=-6:0.5:6;
y1=sqrt(36-x1.^2);%人行走的轨迹
y2=-sqrt(36-x1.^2);%人行走的轨迹
y3=sqrt(36-x2.^2);
y4=-sqrt(36-x2.^2);
X=[x2,x2];
Y=[y3,y4];
y=[y1,y2];
x=[x1,x1];
plot(X,Y,x,y,'o')
M=[x',y'];%声源点的在 xoy 面内的坐标位置位置
fprintf('选取行人轨迹中的点坐标（x,y）分别为：\n');
disp(M);
for i=1:10for j=1:5t1(i,j)=sqrt((x(i)-n(j,1))^2+(y(i)-n(j,2))^2+(n(j,3))^2)/v;end
end
fprintf('不同声源信号到达 n1、n2、n3、n4、n5 的时间\n');
disp(t1);
%以上为模拟仿真过程
#include<iostream>
#include<math.h>
#include<iomanip>
using namespace std;
int main(){
double v=340;
double L1,L2,L3,L4,L5;
double t1,t2,t3,t4,t5;
double r1,r2,r3,r4,r5;
double x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6;
double x7,y7,x8,y8,x9,y9,x10,y10,x,y;
double x11,y11,x22,y22,x33,y33,x44,y44,x55,y55,x66,y66;
double x77,y77,x88,y88,x99,y99,x1010,y1010;
double h,d;
double b1,b2,b3,b4,b5,b6,b7,b8,b9,b10;
double q1,q2,q3,q4,q5,q6,q7,q8,q9,q10;
cout<<"麦克风的高度 h"<<endl;
cout<<"h="<<"";
cin>>h;
cout<<"相临两个麦克风之间的距离"<<endl;
cout<<"d="<<"";
cin>>d;
for(int i=1;i<9;i++){
cout<<"请输入每个麦克风第一次接收到信号的时间 t1,t2,t3,t4,t5"<<endl;
cout<<"t1="<<"";
cin>>t1;
cout<<"t2="<<"";
cin>>t2;
cout<<"t3="<<"";
cin>>t3;
cout<<"t4="<<"";
cin>>t4;
cout<<"t5="<<"";
cin>>t5;
L1=v*t1;
r1=sqrt(L1*L1
-h*h);
L2=v*t2;
r2=sqrt(L2*L2
-h*h);
L3=v*t3;
r3=sqrt(L3*L3
-h*h);
L4=v*t4;
r4=sqrt(L4*L4
-h*h);
L5=v*t5;
r5=sqrt(L5*L5
-h*h);b1=(d*d+r2*r2
-r1*r1)/(2*d*r2);q1=sqrt(1
-b1*b1);x1=r2*q1;y1=r2*b1
-d;cout<<setprecision(4)<<"x1="<<x1<<endl;cout<<setprecision(4)<<"y1="<<y1<<endl;x11=
-1*x1;y11=y1;cout<<setprecision(4)<<"x11="<<x11<<endl;cout<<setprecision(4)<<"y11="<<y11<<endl;b2=(d*d+r1*r1
-r3*r3)/(2*d*r1);q2=sqrt(1
-b2*b2);x2=r1*b2;y2=r1*q2;cout<<setprecision(4)<<"x2="<<x2<<endl;cout<<setprecision(4)<<"y2="<<y2<<endl;x22=x2;y22=
-1*y2;cout<<setprecision(4)<<"x22="<<x22<<endl;cout<<setprecision(4)<<"y22="<<y22<<endl;b3=(d*d+r1*r1
-r4*r4)/(2*d*r1);q3=sqrt(1
-b3*b3);x3=r1*q3;y3=r1*b3;cout<<setprecision(4)<<"x3="<<x3<<endl;cout<<setprecision(4)<<"y3="<<y3<<endl;x33=
-1*x3;y33=y3;cout<<setprecision(4)<<"x33="<<x33<<endl;cout<<setprecision(4)<<"y33="<<y33<<endl;b4=(d*d+r5*r5
-r1*r1)/(2*d*r5);q4=sqrt(1
-b4*b4);x4=r5*b4
-d;y4=r5*q4;cout<<setprecision(4)<<"x4="<<x4<<endl;cout<<setprecision(4)<<"y4="<<y4<<endl;x44=x4;y44=
-y4;cout<<setprecision(4)<<"x44="<<x44<<endl;
cout<<setprecision(4)<<"y44="<<y44<<endl;b5=(2*d*d+r2*r2
-r3*r3)/(2*sqrt(2)*d*r2);q5=sqrt(1
-b5*b5);x5=r2*((sqrt(2)/2)*(b5
-
q5));y5=r2*((sqrt(2)/2)*(b5+q5))
-d;cout<<setprecision(4)<<"x5="<<x5<<endl;cout<<setprecision(4)<<"y5="<<y5<<endl;x55=(x5
-1*x5+2*y5+2)/2;
y55=(2*x5+y5
-
2
-y5)/2;
cout<<setprecision(4)<<"x55="<<x55<<endl;cout<<setprecision(4)<<"y55="<<y55<<endl; b6=(2*d*2*d+r2*r2
-r4*r4)/(4*d*r2);q6=sqrt(1
-b6*b6);x6=r2*q6;y6=r2*b6
-d;cout<<setprecision(4)<<"x6="<<x6<<endl;cout<<setprecision(4)<<"y6="<<y6<<endl;x66=
-1*x6;y66=y6;cout<<setprecision(4)<<"x66="<<x66<<endl;cout<<setprecision(4)<<"y66="<<y66<<endl;b7=(2*d*d+r5*r5
-r2*r2)/(2*sqrt(2)*d*r5);q7=sqrt(1
-b7*b7);x7=r5*((sqrt(2)/2)*(b7+q7))
-d;y7=r5*((sqrt(2)/2)*(q7
-b7));cout<<setprecision(4)<<"x7="<<x7<<endl;cout<<setprecision(4)<<"y7="<<y7<<endl;x77=(x7
-x7+(2*
-1*y7)
-2)/2;y77=(
-2*x7+y7
-
2
-y7)/2;cout<<setprecision(4)<<"x77="<<x77<<endl;cout<<setprecision(4)<<"y77="<<y77<<endl;b8=(2*d*d+r4*r4
-r3*r3)/(2*sqrt(2)*d*r4);q8=sqrt(1
-b8*b8);x8=r4*((sqrt(2)/2)*(b8+q8));y8=r4*((sqrt(2)/2)*(q8
-b8))+d;cout<<setprecision(4)<<"x8="<<x8<<endl;cout<<setprecision(4)<<"y8="<<y8<<endl;x88=(x8
-x8
-2*y8+2)/2;y88=(
-2*x8+y8+2
-y8)/2;cout<<setprecision(4)<<"x88="<<x88<<endl;cout<<setprecision(4)<<"y88="<<y88<<endl;b9=(4*d*d+r5*r5
-r3*r3)/(4*d*r5);q9=sqrt(1
-b9*b9);x9=r5*b9
-d;y9=r5*q9;cout<<setprecision(4)<<"x9="<<x9<<endl;cout<<setprecision(4)<<"y9="<<y9<<endl;x99=x9;y99=
-1*y9;cout<<setprecision(4)<<"x99="<<x99<<endl;cout<<setprecision(4)<<"y99="<<y99<<endl;b10=(2*d*d+r5*r5
-r4*r4)/(2*sqrt(2)*d*r5);q10=sqrt(1
-b10*b10);x10=r5*((sqrt(2)/2)*(b10
-q10))
-d;y10=r5*((sqrt(2)/2)*(b10+q10));cout<<setprecision(4)<<"x10="<<x10<<endl;cout<<setprecision(4)<<"y10="<<y10<<endl;x1010=(x10-x10+2*y10-2)/2;y1010=(2*x10+y10+2-y10)/2;cout<<setprecision(4)<<"x1010="<<x1010<<endl;cout<<setprecision(4)<<"y1010="<<y1010<<endl;}
return 0;
}
