第10章-声学仿真试验

1、第10章 声学仿真试验,10.1 回声试验 10.2 多普勒效应 10.3 声学信号的滤波特性试验 10.4 交混回响 10.5 短时傅立叶变换,10.1 回声试验,10.1.1 回声的基本原理 回声是一种物理现象。当直达声结束后，声音经其它物体反射、折射，又返回来收到的声波，称为回声。回声有时泛指一切反射回来的声音。当传到人耳的直达声和回声之间的时间差在1/20s以上时，可以很清楚地把它们区分开。在山谷或大厅中，常有回声现象发生。尤其是在影剧院里，回声往往会妨碍听音，所以建筑师必须考虑消除回声的影响。相反，直达声和回声之间的时间差在1/20s以下时，回声不但不妨碍听音，反而使声音更加响亮。,

2、应用MATLAB软件可以方便地进行仿真试验。回声现象取决于几个重要的物理参数： (1)声源距离听者有多远。 (2)声音产生往复反射的两个反射体距离有多远。 (3)声音的反射体在反射声音的过程中对声音的衰减有多大。(假设听者在声源与反射体之间的某个位置上。),10.1.2 回声的仿真试验 图10-1所示是回声的仿真试验系统。图中1800的延迟代表了声音在两个反射体之间走一个来回的延迟（1800/Fs秒），它的值愈大，反映了两个反射体的距离愈大。图中0.8的增益代表了声音在两个反射体之间走一个来回的损耗，它愈小，回声衰减愈快。图中0.9的增益代表了传输损耗。图10-2所示是仿真试验时示波器上显示的

3、波形图，上图是原声，下图是往复反射的回声波形图。如果计算机装有声卡和音箱，可以听到扬声器发出的回声声音。放大器增益和整数延迟，可以看图标注。,图10-1 回声仿真试验系统框图,图10-2 回声仿真试验时域图,回声的仿真试验的模型参数设置分析：回声的仿真试验的模型参考了MATLAB中Toolboxdspblksdspdemosdspafxr的DEMOS例子，在仿真系统中设置了Reshape（整形器），主要参数见表10-1，ToWaveDevice（扬声器）的主要参数见表10-2，FromWorkspace（声源）的主要参数见表10-3，它以帧的形式输出。整形器的设置目的是将信号转变为示波器能够接

4、受的1维数据流。,表10-1 Reshape（整形器）的主要参数,表10-2 ToWaveDevice（扬声器）的主要参数,表10-3 FromWorkspace（来自工作空间）的主要参数,在10.1、10.4、10.5节中都需要一段声音文件来进行仿真试验。可以有两种方法获得声音文件： (1)调出MATLAB中的声音文件“loadmtlb.mat”。 (2)应用6.2.2节介绍的声音采集的方法，建立一个声音文件。譬如文件名是Soud11.mat。 通过下列程序调入仿真的声音数据并启动仿真。Echo11是图10-1所示的仿真系统的名字（文件名）。 程序10-1 loadsoud11 %读入声音数

5、据文件soud11 sim（echo11）,10.2 多普勒效应,10.2.1 多普勒效应的基本原理 当波源或观察者相对于媒质运动时，或者说波源和观察者有相对运动时，观察者接受到的振动频率与波源振动频率不同的现象，称为多普勒效应。 当听者与声源的相对速度大到与声速可以比拟时，就可以明显地感觉到声音频率的变化。,譬如坐火车时，听到相向开来的火车汽笛声，当火车由远而近开来，汽笛声的音调变高；由近而远离去，汽笛声的音调变低。这是日常生活中的一个多普勒效应的例子。在天文、通信等领域还有众多的例子。 以下的公式描述了该现象的各个物理量之间的定量关系： 其中,f0是声源发出的声音的频率，v是听者与声源相对

6、运动的速度，为速度矢量与声源和听者的连线夹角，vs为声音在空气中传播的速度，f是听者听到的声音频率。 ,(10-1）,10.2.2 多普勒效应的仿真试验 Simulink仿真模型是参考MATLAB中的Toolboxdspblksdspdemosdspstfft建立的。 下面是一段用MATLAB的M文件编写的程序。它产生一个文件名为u1的表现多普勒效应的声音文件，描述火车向一个距铁路30m(垂直距离）、距火车510m（水平距离）的听者开来时他听到的声音。将u1送入图10-3所示的仿真系统，再一次听到该声音，并且看到用两种方法（频谱仪和矩阵图仪）表现的声音的时频图。表10-4表10-11分别是图1

7、0-3所示仿真系统中各个模块的主要参数。图10-4所示是M文件程序运行结束时绘出的声音信号的原声波形（上图）和多普勒效应的波形（下图）。图10-5所示是某个时刻图10-3中的频谱仪（左图）和矩阵仪（右图）上显示的时频图。,图10-3 显示声音信号(多普勒响应)的仿真系统模型,表10-4 FromWorkspace（来自工作空间）的主要参数,表10-5 Rebuffer（缓存器）的主要参数,表10-6 ShortTimeFFT（短时傅立叶变换）的主要参数,表10-7 FrequencyFrameScope（频率帧示波器）的主要参数,表10-8 MatrixViewer（矩阵显示图仪）的主要参数,

8、表10-9 Select（选择器）的主要参数,表10-10 dB（分贝）的主要参数,表10-11 ToWaveDevice（扬声器）的主要参数,图10-4 声源发出的信号（上图）和 听者接收到的信号（下图，多普勒效应）,图10-5 频谱仪显示（左图）和矩阵表示（右图）,相应程序如下： 程序10-2 %多普勒效应程序 x0=550; v=150;%车速 y=30; c=330;%声音在空气中的传播速度 w=2200;%声音的频率 t=0:.0005:30;,r=sqrt（x0-v*t）.2+（y）.2）;%计算声源与听者距离 t1=t-r/c;%经距离迟延后听者的等效时间 u=1.5*sin（w

9、*t）;%声源发出的信号 u1=1.5*sin（w*t1）;%听者接收到的信号 sound（u，10000）;pause（5）;sound（u1，10000）;%将原信号和接收到的信号恢复为声音 figure（1） subplot（2，1，1） plot（u）;axis（11500，11700，-2.2，2.2） subplot（2，1，2） plot（u1，r）;axis（5500，7800，-2.2，2.2）,10.3 声学信号的滤波特性试验,10.3.1 声音滤波的应用 数字信号处理课程中详细地讨论了各种数字滤波器的设计方法，本例是数字滤波器应用在声学领域的一个例子。原始的声音由三个单音

10、频组成，采用低通滤波器，可以将高音频的信号滤除；采用带通滤波器，可以只让某一个频率的音频信号通过。编程的要点是在正确地选择采样频率的基础上，正确地选取通带和阻带的截止频率，以保证滤掉和保留下预想的频率。,10.3.2 声音滤波的仿真试验 下面是一段声音滤波的仿真试验的程序。300Hz、500Hz、1000Hz的三个正弦波信号叠加构成本例的原信号，通过低通滤波器后滤掉了1000Hz的信号。图10-6中的上面两幅图，是表示低通滤波器传输特性的幅频图和相频图。通过带通滤波器后只剩下1000Hz的信号。图10-6中的下面两幅图，是表示带通滤波器传输特性的幅频图和相频图。,图10-6 滤波器的幅频特性和

11、相频特性(上，低通；下，带通）,图10-7所示是信号的频谱特性。从上到下依次是三个音频信号；通过低通滤波器后，1000Hz的信号滤掉了，仅剩300Hz和500Hz的信号；通过中心频率为1000Hz带通滤波器后，仅剩1000Hz的信号。程序运行结束时，可以听到依次发出的三种信号的声音。,图10-7 信号的频谱特性 （上，原信号；中，通过低通；下，通过带通）,程序10-3 t=1/10000:1/10000:3;%设置三个频率的正弦信号 %300Hz，500Hz，1000Hz s1=sin（2*pi*300*t）; s2=sin（2*pi*500*t）; s3=1.2*sin（2*pi*1000*

12、t）; s=s1+s2+s3;%三个正弦波信号叠加 figure（1）,subplot（3，1，1）; ss=fft（s，4096）;SS=abs（ss（1:1:2049）;%求合成信号的频谱 k1=0:2048;w1=（1/.4096）*k1;%取0.Fs/2的部分 plot（w1，SS）;grid%画频谱图 axis（0，2000，-100，2200） title（原信号频谱图） %通过低通滤波器 ws1=1000;%设计一个通带为600Hz、阻带为1000Hz的低通滤波器,wp1=600; wc=5000; wp=wp1/wc;ws=ws1/wc; n，Wn=cheb2ord（wp，ws

13、，1，30）%切比雪夫II型滤波器设计 b，a=cheby2（n，30，Wn）; subplot（3，1，2） sb=3*filter（b，a，s）;%合成信号通过低通滤波器 ssb=fft（sb，4096）;SSb=abs（ssb（1:1:2049）;%求频谱 k1=0:2048;w1=（1/.4096）*k1;%画频谱图 plot（w1，SSb）;grid axis（0，2000，-100，6500） title（经过低通滤波器后的信号频谱图）,%通过带通滤波器 ws1=6001400;%设计一个阻带为6001400Hz、 wp1=9801020;%通带为9801020Hz的带通滤波器 w

14、c=5000; wp=wp1/wc;ws=ws1/wc; n，Wn=cheb1ord（wp，ws，4，45，s）%切比雪夫I型滤波器 bd，ad=cheby1（n，4，Wn，bandpass）; subplot（3，1，3） sd=3*filter（bd，ad，s）;%合成信号通过低通滤波器,ssd=fft（sd，4096）;SSd=abs（ssd（1:1:2049）;%求频谱 k1=0:2048;w1=（1/.4096）*k1; plot（w1，SSd）;grid%画频谱 axis（0，2000，-100，4000） title（经过带通滤波器后的信号频谱图） figure（2）%画低通滤波

15、器的幅频、相频图 freqz（b，a，4096，10000）;axis（0，4000，-70，1）,figure（3）%画带通滤波器的幅频、相频图 freqz（bd，ad，4096，10000）;axis（0，4000，-70，1） %依次用扬声器播放合成音、通过低通后、通过带通后的声音 sound（s，10000）;pause（5）;sound（sb，10000）; pause（5）;sound（sd，10000）,10.4 交混回响,10.4.1 交混回响的应用 在一间有若干个扬声器的礼堂里，处在某一个特定位置听到的声音是什么效果，可以用仿真的方法来实现。为了简化问题，假定墙面上有吸音的材

16、料，即忽略了墙壁的反射。编程仿真的思想是：电信号在电线里传输的时间可以忽略，但是当相距若干米的扬声器发出的声音到达听者的位置时，不同扬声器的声音经过了不同的路程，路程的差别产生相位差，所有带有不同相位的声音的叠加就产生了交混回响的效果。,10.4.2 交混回响的仿真试验 图10-8所示是交混回响仿真试验中的扬声器布置图，A、B、C是三个接到同一声源的扬声器，它距离听者的水平和垂直距离标注在图上。下面是用M文件编制的程序，它的运行绘出了原声和交混回响的声音时域波形。图10-9所示是声音信号的时域图，并且依次播放出原声和交混回响的声音效果。,图10-8 交混回响仿真试验中的扬声器布置,图10-9声

17、音信号的时域图 （上，原信号；下，听到的交混回响的信号）,程序10-4 loadsoud11%读入声音数据文件soud11 y=soud11;fs=8000; N1=8600; v=330;%声音速度 b1=（202+452）.5;%计算听者到扬声器B的距离 yb1=zeros（1，(fs*b1/v），y;%听者听到扬声器B的声音 yb=yb1（1:N1）; a1=（202+152）.5;%计算听者到扬声器A的距离 ya1=zeros（1，(fs*a1/v），y;%听者听到扬声器A的声音,ya=ya1（1:N1）; c1=（352+652）.5;%计算听者到扬声器C的距离 yc1=zeros（

18、1，(fs*c1/v），y;%听者听到扬声器C的声音 yc=yc1（1:N1）; ye=1.6*ya+2*yb+1.2*yc;%听者听到三个扬声器的声音叠加 subplot（211）;plot（y）;grid subplot（212）;plot（ye）;grid sound（y，fs）;pause（2）;sound（1.2*ye，fs)%播放原声后播放合成的声音,10.5 短时傅立叶变换,10.5.1 短时傅立叶变换的基本原理 傅立叶变换将信号系统的时间域和频率域的特性联系起来，这一工具在处理平稳信号时，因为信号的统计特性（平均值、方差以及频率特性）都不是时间的函数，因此，在求其频谱特性时，应

19、用了从负无穷到正无穷的时间积分。,(10-2）,当信号是一个非平稳信号时，它的频率特性是时间的函数，上面的方法就不能正确地反映出信号的特性，因为在不同的时段内有不同的频率特性。应用短时傅立叶变换： (10-3） 可以正确地反映出非平稳信号的频率特性，此时的频谱特性已经是时间和频率的二维函数。 公式中的W可以看成是一个窗函数，短时傅立叶变换的基本思路是把非平稳信号分割为若干小的时段，在每个小时段内把信号看成是平稳的。各个时段的频率特性的排列，就构成了信号的时频特性。时窗愈窄，它的时间分辨率愈高。,10.5.2 短时傅立叶变换的仿真试验 Simulink仿真模型是参考MATLAB中的Toolbox

20、dspblksdspdemosdspstfft构建的，并且设置了ToWaveDevice（扬声器）和WaterfallScope（瀑布图仪），如图10-10所示。,图10-10 显示声音信号(用于短时傅立叶变换)的仿真系统模型,下面程序说明：首先调出声音文件Soud11；然后运行mdl文件“dspstfftme3”，把“dspstfftme3”的声音信号的数据用ToWorkspace模块收集，变量名为sA，将数据序列sA作分段傅立叶变换，目的是将声音信号的时间频率特性用三维图形表示出来。 程序10-5 loadsoud11 sim（dspstfftme3） S，F，T=specgram（sA，128，1/8000，16）;%求解短时傅立叶变换%将数据序列sA作N点傅立叶变换，Fs，Nbit的意