《MATLAB仿真在通信与电子工程中的应用》课件第10章

第10章声学仿真试验10.1回声试验10.2多普勒效应10.3声学信号的滤波特性试验10.4交混回响10.5短时傅立叶变换10.1回声试验

10.1.1回声的基本原理

回声是一种物理现象。当直达声结束后，声音经过其他物体反射、折射后，又返回来收到的声波，称为回声。回声有时泛指一切反射回来的声音。当传到人耳的直达声和回声之间的时间差在1/20s以上时，可以很清楚地把它们区分开。在山谷或大厅中，常有回声现象发生。尤其是在影剧院里，回声往往会妨碍听音，所以建筑师必须考虑消除回声的影响。相反，直达声和回声之间的时间差在1/20s以下时，回声不但不妨碍听音，反而使声音更加响亮。应用MATLAB软件可以方便地进行仿真试验。回声现象取决于几个重要的物理参数：

(1)声源距离听者有多远。

(2)声音产生往复反射的两个反射体的距离有多远。

(3)声音的反射体在反射声音的过程中对声音的衰减有多大。(假设听者在声源与反射体之间的某个位置上。)10.1.2回声的仿真试验

图10-1所示是回声的仿真试验系统框图。图中的1800的延迟代表了声音在两个反射体之间走一个来回的延迟（1800/Fs秒）。它的值愈大，反映了两个反射体的距离愈大。图中0.8的增益代表了声音在两个反射体之间走一个来回的损耗，它愈小，回声衰减愈快。图中0.9的增益代表了传输损耗。图10-2所示是仿真试验时示波器上显示的波形图，上图是原声，下图是往复反射的回声波形图。如果计算机装有声卡和音箱，可以听到扬声器发出的回声声音。放大器增益和整数延迟可以看图标注。图10-1回声仿真试验系统框图图10-2回声仿真试验时域图回声的仿真试验的模型参数设置分析：回声的仿真试验的模型参考了MATLAB中Toolbox＼dspblks＼dspdemos＼dspafxr的DEMOS例子，在仿真系统中设置了Reshape（整型器），主要参数见表10-1；ToWaveDevice（扬声器）的主要参数见表10-2；FromWorkspace（声源）的主要参数见表10-3；它以帧的形式输出。整型器的设置目的是将信号转变为示波器能够接受的一维数据流。在10.1、10.4、10.5节中都需要一段声音文件来进行仿真试验。可以用以下两种方法来获得声音文件：

(1)调出MATLAB中的声音文件“loadmtlb.mat”。

(2)应用6.2.2节介绍的声音采集的方法，建立一个声音文件。譬如文件名是Sound11.mat。

通过下列程序调入仿真的声音数据并启动仿真。Echo11是图10-1所示的仿真系统的名字（文件名）。

程序10-1

loadsound11%读入声音数据文件sound11

sim（′echo11′）

10.2多普勒效应

10.2.1多普勒效应的基本原理

当波源或观察者相对于媒质运动时，或者说波源和观察者有相对运动时，观察者接受到的振动频率与波源振动频率不同的现象，称为多普勒效应。

当听者与声源的相对速度大到与声速可以比拟时，就可以明显地感觉到声音频率的变化。譬如坐火车时，听到相向开来的火车汽笛声，当火车由远而近开来时，汽笛声的音调变高；由近而远离去时，汽笛声的音调变低。这是日常生活中的一个多普勒效应的例子。在天文、通信等领域还有众多的例子。以下的公式描述了该现象的各个物理量之间的定量关系:(10-1）其中,f0是声源发出的声音的频率，v是听者与声源的相对运动速度，θ为速度矢量与声源和听者的连线夹角，vs为声音在空气中传播的速度，f是听者听到的声音频率。10.2.2多普勒效应的仿真试验

下面是一段用MATLAB的M文件依据公式（10-1）原理编写的程序。它会产生一个文件名为dopp.wav的表现多普勒效应的声音文件，描述火车向一个距铁路30m（垂直距离）、距火车510m（水平距离）的听者开来时他听到的声音。将u1送入图10-3所示的Simulink仿真系统，再一次听到该声音，并且看到用频谱仪表现的声音频率随时间变化的情况。表10-4和表10-5给出了仿真系统中两模块的主要参数设置。

图10-4所示是M文件程序运行结束时绘出的听者接收到信号的频率变化曲线以及用信号处理工具箱函数spectrogram绘制的时间频率图。图10-5所示是某个时刻图10-3中的频谱仪上显示的功率谱。图10-3显示声音信号（多普勒效应）的仿真系统模型图10-4听者接收到的信号的频率变化曲线图10-5频谱仪显示（左图）和时间频率图表示（右图）程序10-2

%多普勒效应程序

x0=150;%m

v=50;%车速m/s

y0=30;%m

c=330;%声音在空气中的传播速度m/s

f0=2000;%声音的频率Hz

fs=8000;%采样率

t=0:1/fs:6;

x_t=－x0+v.*t;

y=y0;

r=sqrt(x_t.^2+y.^2);%计算声源与听者的距离

costheta=x_t./r;

f=f0./(1+v.*costheta./c);

fmin=min(f)

fmax=max(f)

voltage=(f-fmin)./(fmax-fmin)*2－1;%归一化调频电压在-1~+1之间

signal=0.5*vco(voltage,［fminfmax］,fs);

u=sin(2*pi*f0*t);%声源发出的信号

sound(u,fs);pause(5);%播放声源发出的信号

sound(signal,fs);%播放听者接收到的信号

wavwrite(signal,fs,′dopp.wav′);

figure(1);

plot(t,f);xlabel(′Time/s′);

ylabel(′Freq/Hz′);

figure(2);

gram(signal,kaiser(256,5),220,512,fs,′yaxis′);%作出时间频率图

axis(［06,1000,spectro3000］)

10.3声学信号的滤波特性试验

10.3.1声音滤波的应用

数字信号处理课程中详细讨论了各种数字滤波器的设计方法，本例是数字滤波器应用在声学领域的一个例子。原始的声音由三个单音频组成，采用低通滤波器，可以将高音频的信号滤除；而采用带通滤波器，可以只让某一个频率的音频信号通过。编程的要点是在正确地选择采样频率的基础上，正确地选取通带和阻带的截止频率，以保证滤掉和保留下预想的频率。10.3.2声音滤波的仿真试验

下面是一段声音滤波的仿真试验的程序。300Hz、500Hz、1000Hz的三个正弦波信号叠加构成本例的原信号，通过低通滤波器后滤掉了1000Hz的信号。图10-6中上面的两幅图分别是表示低通滤波器传输特性的幅频图和相频图。通过带通滤波器后只剩下1000Hz的信号。图10-6中下面的两幅图分别是表示带通滤波器传输特性的幅频图和相频图。图10-6滤波器的幅频特性和相频特性（左图：低通；右图:带通）图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*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，1，30）%切比雪夫Ⅱ型滤波器设计

［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］;

%小于等于600Hz、大于等于1400Hz的区间为阻带

wp1=［9801020］;%980~1020Hz区间为通带

wc=5000;

wp=wp1/wc;ws=ws1/wc;［n，Wn］=cheb1ord（wp，ws，4，45，′s′）%切比雪夫Ⅰ型滤波器

［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）%画低通滤波器的幅频、相频图

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交混回响的应用

在一间有若干个扬声器的礼堂里，处在某一个特定位置听到的声音是什么效果，可以用仿真的方法来实现。为了简化问题，假定墙面上有吸音的材料，即忽略了墙壁的反射。编程仿真的思想是：电信号在电线里的传输时间可以忽略，但是当相距若干米的扬声器发出的声音到达听者的位置时，不同扬声器的声音经过了不同的路程，路程的差别产生相

位差，所有带有不同相位的声音的叠加就产生了交混回响的效果。10.4.2交混回响的仿真试验

图10-8所示是交混回响仿真试验中的扬声器布置图，A、B、C是三个接到同一声源的扬声器，它距离听者的水平和

垂直距离标注在图上。下面是用M文件编制的程序，它的运行绘制出了原声和交混回响的声音时域波形。图10-9所示是声音信号的时域图，并且依次播放出原声和交混回响的声音效果。图10-8交混回响仿真试验中的扬声器布置图10-9声音信号的时域图（上图：原信号；下图：听到的交混回响的信号）程序10-4

loadsound11%读入声音数据文件sound11

y=sound11;fs=8000;

N1=8600;

v=330;%声音的速度

b1=（20^2+45^2）^.5;%计算听者到扬声器B的距离

yb1=［zeros（1，（fs*b1/v）），y′］′;

%听者听到扬声器B的声音

yb=yb1（1:N1）;

a1=（20^2+15^2）^.5;%计算听者到扬声器A的距离

ya1=［zeros（1，（fs*a1/v）），y′］′;

%听者听到扬声器A的声音

ya=ya1（1:N1）;

c1=（35^2+65^2）^.5;%计算听者到扬声器C的距离

yc1=［zeros（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短时傅立叶变换的基本原理

傅立叶变换将信号系统的时间域和频率域的特性联系起来，这一工具在处理平稳信号时，由于信号的统计特性（平均值、方差以及频率特性）都不是时间的函数,因此，在求其频谱特性时，应用了从负无穷到正无穷的时间积分。(10-2）当信号是一个非平稳信号时，它的频率特性是时间的函数，上面的方法就不能正确地反映出信号的特性，因为在不同的时段内信号有不同的频率特性。

应用短时傅立叶变换：(10-3）可以正确地反映出非平稳信号的频率特性，此时的频谱特性已经是时间和频率的二维函数。公式中的W可以看成是一个窗函数，短时傅立叶变换的基本思路是把非平稳信号分割为若干小的时段，在每个小时段内把信号看成是平稳的。各个时段的频率特性的排列，就构成了信号的时频特性。时窗愈窄，它的时间分辨率愈高。10.5.2短时傅立叶变换的仿真试验

Simulink仿真模型是参考MATLAB中的Toolbox＼dspblks＼dspdemos＼dspstfft构建的，并且设置了ToWaveDevice（扬声器）和WaterfallScope（瀑布图仪），如图10-10所示。

下面的程序说明：首先调出声音文件Sound11；然后运行mdl文件“dspstfftme3”，把“dspstfftme3”的声音信号的数据用ToWorkspace模块收集，变量名为sA，