回声信号地产生与消除

数字信号处理课程设计

回声信号的产生与消除

姓名张针海

学号 10300123

专业电子信息工程

指导教师樊玲

年级 10级电信2班

日期 2013 .5 . 25

【摘要】本课程是利用Windows下的录音机，录制一段自己不小于10s的语音，然后在Matlab 软件平台下，利用函数wavread对语音信号进行采样，并记录采样频率和采样点数。在抽样信号的基础上，通过采样后的的信号与原信号实现一次及多次延迟、叠加产生回波信号，再使用Matlab绘出有回声及无回声语音信号的时域波形和频谱图。再分别用频率抽样法设计的FIR滤波器和冲激相应不变法设

计设计的IIR滤波器消除回声，并记录滤波器的频域响应，再绘制滤波后信号的时域波形和频谱，并对前后信号进行对比，分析信号的变化。

[关键词] 录音 matlab 采样滤波抽样

[Abstract] this course is to use a tape recorder to record voice under Windows, a section of their own not less than 10s, then in Matlab software platform, sampling of the speech signal using the function wavread, and record the sampling frequency and sampling points. Based on the sampling signal, through its implementation of single and multiple superposition delay, echo, and use Matlab to draw the echo and echo free speech signal time-domain waveform and spectrum. FIR filter respectively by frequency sampling design method and impulse corresponding invariant IIR filter design to eliminate echo, and record the response of the filter in frequency domain, and then draw the time-domain waveform and spectrum of the filtered signal, and compared before and after the signal, analysis of signal changes

目录

1 设计目的及要求 (3)

1.1设计回音目的及要求 (3)

1.2设计滤波器目的及要求 (3)

1.2.1 FIR滤波器 (3)

1. 2. 2 巴特沃兹滤波器 (3)

1. 2. 3 距离估计要求 (3)

2 设计原理 (4)

3设计内容 (4)

3．1语音采集........ (4)

3. 2信号分析 (4)

3．3制作回音 (5)

3.4设计滤波器及滤波 (8)

3. 4. 1 设计FIR滤波器及滤波 (8)

3.4.1.1单回声的滤波 (8)

3. 4.2设计巴特沃兹滤波器及滤波 (11)

3.4.2.1设计巴特沃斯数字低通滤波器 (11)

3．5估算距离 (12)

3.5.1通过理论计算法 (12)

3.5.2程序返回测量法 (13)

4总结 (14)

5、参考文献 (15)

1 设计目的及要求

1.1设计回音目的及要求

现代通信中回波是影响通信质量的噪声，本课程设计是在matble库元件中搜索一段不小于10s的录音，再利用函数wavread对语音信号进行采样，并自身实现一次及多次延迟、叠加产生回波信号，再使用Matlab绘出有回声及无回声语音信号的时域波形和频谱图。在此过程中必须灵活运用matlab中 wavread函数对信号进行采样，同时加深对声频信号中噪声的认识。

1.2设计滤波器目的及要求

1.2.1 FIR滤波器

FIR滤波器是有限长冲激响应滤波器，它是通信，语音与图象处理模式识别及频谱分析等应用中一种基本的处理部件。他可以满足滤波器幅度和相位特性的严格要求避免模拟滤波器无法克服的电压漂移，温度漂移和噪声等问题。数字滤波器的作用是滤除信号中某一部分频率分量。信号经过滤波处理。就是相当于信号频谱与滤波响应相乘的结果。时域上看就是信号与滤波器的冲激响应卷积的结果。有限长冲激响应滤波器可以保证任何幅频特性的同时具有严格的相频特性，同时其单位冲激响是有限的。没有输入到输出的反馈，是稳定的系统。

1. 2. 2 巴特沃兹滤波器

欲消除回声，则需设计滤波器，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化，理想低通滤波器的一个非因果系统，如果允许低通滤波器的通带和阻带之间有一定的过渡带，且通带和阻带允许有一定的衰减，我们就可以用物理可实现的系统去逼近理想低通滤波器的频率特性，从而获得较好的滤波效果；巴特沃兹滤波器就是工程中常用的频率响应逼近理想低通滤波器的物理可实现系统，还原出原有音频信号。

1. 2. 3 距离估计要求

从信号y中估计反射物的距离，即回声延迟的时间，可理解为估算原有信号中的N值。也就是，估测y(n)中的原始声音信号x(n)与其延时衰减分量kx(n-N)的相关联的程度。

2 设计原理

语音采集：通过搜索在pc机的中matble库函数中chimes.wav，绘制其时域波形，对此音频信号用FFT作谱分析。

制作回音：利用MATLAB软件中的wavread函数进行采样分析得到的源噪声，通过叠加到原始语音

信号中，模仿语音信号被污染，并对其进行频谱分析；

消除回音：设计FIR和IIR数字滤波器的滤波器，对带有回声的声音信号进行滤波，恢复出原始信号。绘制所设计滤波器的幅频和相频特性，及滤波后的信号的

时域波形和频谱图，分析滤波后信号的时域和频域特征，回放语音信号。

估算距离：计算信号y中估计反射物的距离，可理解为估计延时序列中的N值。即估计y(n)中的原始声音信号x(n)与其延时衰减分量kx(n-N)的相关联的程度；

相关是指两个确定信号或两个随机信号之间的相互关系，对于随机信号，信

号一般是不确定的，但是通过对它的规律进行统计，它们的相关函数往往是

确定的，故在随机信号处理中，可以用相关函数来描述一个平稳随机信号的

统计特性，因而可通过相关分析法估算反射物的距离。

3设计内容

3．1语音采集

读取本地matble库中的自带的文件chimes.Wav。

3. 2信号分析

在matbal中绘制其时域波形，利用wavread()函数将其提取出来绘制其时域波形对此音频信号用FFT作谱分析，用plot()函数绘制其图形。

%声音信号的提取

fs=22050;

[x,fs,Nbits]=wavread('D:\chimes.wav'); %把语音信号进行加载入Matlab仿真软件平台中wavplay(x,fs);% 回放语音信号。或者sound(x,fs)

figure(1);

N=length(x);%求语音信号的长度

subplot(3,1,1);

plot(x(1:N));

title('原始信号波形');

y=fft(x,N);%傅立叶变换

subplot(3,1,2);

plot((0:N-1)/N*fs,abs(y));

title('原始信号幅值');

subplot(3,1,3);

plot(angle(y));

title('原始信号相位');

2，产生的原信号的波形，以及其幅度、相位谱如图1所示：

02000400060008000100001200014000

-0.50

0.5原始信号波形

0.5

1

1.5

2

2.5x 10

4

0100

200原始信号幅值

02000400060008000100001200014000

-5

5原始信号相位

图1 原信号波形图

3．3制作回音

1 参数的设置：

因为人耳能分辨出的声音延迟至少是0.1s ，因此，最小延迟量不能小于0.1s 。在此先先延迟时间为0.2s,即最小延迟量N=0.2*fs=0.2*16000=3200。

在已有声音信号x 的基础上产生带回声的声音信号，可以表达为在于娜信号的基础上叠加其延时的分量。假设只有一个回声的情况下，可简化其模型为

y(n)=x(n)+kx(n-m) k 为反射系数；m 为延迟时间。 这里设m=3200 ； k=0.41 2 声音延迟：

,利用矩阵置零产生x 的延迟，以及得到y 信号：

为了保证图像的完整性，对读取的信号先延长4000个采样点，将原始信号延长3200个采样点，然后再在后面补上800个点，得到如下代码： 1—— %一次回波的产生

[x,fs]=wavread('D:\chimes.wav');%把语音信号进行加载入Matlab 仿真软件平台中。

N=length(x);%语音信号的长度。 x1=x(1:N); x2=x(1:N);

x1=[x1,zeros(1,4000)];%zeros(1,4000)产生1行3000列全零矩阵加到x1后面。 x2=[zeros(1,3200),0.4*x2,zeros(1,800)];%N+3200+x=N+4000，得x=800 y=x1+x2;%加入回音的信号。 figure(2);

subplot(3,1,1);

plot(y(1:3200.+N)); title('含回声信号波形'); y1=fft(y);

subplot(3,1,2); plot(abs(y1));

title('含回声信号幅值'); subplot(3,1,3); plot(angle(y1)); title('含回声信号相位'); sound(y,fs);

加入回声以后的波形，如图2所示：

020004000600080001000012000140001600018000

-0.50

0.5含回声信号波形

020004000600080001000012000140001600018000

0100

200含回声信号幅值

020004000600080001000012000140001600018000

-5

5含回声信号相位

图2 含一次回声波形图

2——%多次回声的产生

[x,fs]=wavread('D:\chimes.wav');%把语音信号进行加载入Matlab 仿真软件平台中。

N=length(x);%语音信号的长度。 x1=x(1:N); x2=x(1:N); x3=x(1:N);

x1=[x1,zeros(1,6400)];%zeros(1,4000)产生1行3000列全零矩阵加到x1后面。 x2=[zeros(1,3200),0.4*x2,zeros(1,3200)];%N+3200+x=N+4000，得x=800 x3=[zeros(1,3200),0.7*x3,zeros(1,3200)]; y=x1+x2+x3;%加入回音的信号。 figure(2);

subplot(3,1,1);

plot(y(1:3200.+N)); title('含回声信号波形'); h=fft(y,N);%傅立叶变换 subplot(3,1,2);

plot((0:N-1)/N*fs,abs(h)); title('含回声信号幅值'); subplot(3,1,3); plot(angle(y));

title('含回声信号相位'); sound(y,fs);

程序设计如图3所示：

020004000600080001000012000140001600018000

-0.50

0.5含回声信号波形

0.5

1

1.5

2

2.5x 10

4

0200

400含回声信号幅值

00.51 1.52 2.5x 10

4

-5

5含回声信号相位

图3 含多次回声波形图

3.4设计滤波器及滤波

3. 4.1 设计FIR滤波器及滤波

3.4.1.1单回声的滤波

%单个回声的滤波

[x,fs]=wavread('D:\chimes.wav');%把语音信号进行加载入Matlab仿真软件平台中。

N=length(x);%语音信号的长度。

x1=x(1:N);

x2=x(1:N);

x3=x(1:N);

subplot(6,1,1);

plot(x);

title('原信号波形');

h=fft(x,N);%傅立叶变换

subplot(6,1,2);

plot((0:N-1)/N*fs,abs(h));

title('原始信号幅值');

x1=[x1,zeros(1,4000)];%zeros(1,4000)产生1行4000列全零矩阵加到x1后面。

x2=[zeros(1,3200),0.8*x2,zeros(1,800)];%N+3200+x=N+4000，得x=800

y=x1+x2;%加入回音的信号。

subplot(6,1,3);

plot(y(1:3200.+N));

title('含回声信号波形');

y1=fft(y);

subplot(6,1,4);

plot(abs(y1));

title('含回声信号幅值');

b=1;

a=zeros(1,N);

a(1)=1;

a(3201)=0.4;

z=filter(b,a,y);

subplot(6,1,5);

plot(z);

title('滤波后信号波形');

z1=fft(z);

subplot(6,1,6); %3行1列排列第1个图

plot(abs(z1)); %画出abs（z1）

title('滤波后信号的幅值'); %绘图区标题

sound(z1,fs);

程序设计后运行结果如图4所示：

2000

4000

6000800010000

12000

14000

-0.50

0.5原信号波形

00.51 1.52 2.5x 10

4

100

200原始信号幅值

2000

4000

6000

80001000012000140001600018000

-0.50

0.5含回声信号波形

2000

4000

6000

80001000012000140001600018000

0200

400含回声信号幅值

2000

4000

6000800010000120001400016000

18000

-0.50

0.5滤波后信号波形

2000

4000

6000

8000

10000120001400016000

18000

0100

200滤波后信号幅值

图4 滤除单个回波后的波形图

3.4.1.2多次回声的滤波

%多个回声的滤除

[x,fs]=wavread('D:\chimes.wav');%把语音信号进行加载入Matlab 仿真软件平台中。 N=length(x);%语音信号的长度。 x1=x(1:N); x2=x(1:N); x3=x(1:N); subplot(6,1,1); plot(x);

title('原信号波形');

h=fft(x,N);%傅立叶变换

subplot(6,1,2);

plot((0:N-1)/N*fs,abs(h));

title('原始信号幅值');

x1=[x1,zeros(1,2400)];%zeros(1,4000)产生1行3000列全零矩阵加到x1后面。x2=[zeros(1,3200),0.4*x2,zeros(1,3200)];%N+3200+x=N+6400，得x=3200

x3=[zeros(1,5200),0.7*x3,zeros(1,5200)];

y=x1+x2+x3;%加入回音的信号。

subplot(6,1,3);

plot(y(1:3200.+N));

title('含回声信号波形');

y1=fft(y);

subplot(6,1,4);

plot(abs(y1));

%plot((0:N-1)/N*fs,abs(y1));

title('含回声信号幅值');

%subplot(6,1,5);

%plot(angle(y1));

%title('含回声信号相位');

b=1;

a(1)=1;

a(3201)=0.4;

a(5201)=0.7;

z=filter(b,a,y);

subplot(6,1,5);

plot(z);

title('滤波后信号波形');

z1=fft(z);

subplot(6,1,6); %3行1列排列第1个图

%plot((0:N-1)/N*fs,abs(z1));

plot(abs(z1)); %画出abs（z2

title('滤波后信号的幅值'); %绘图区标题

sound(z1,fs);

程序设计后运行结果如图5所示：

:

02000400060008000100001200014000 -0.5

0.5

原信号波形

00.51 1.52 2.5

x 104 0

100

200

原始信号幅值

020004000600080001000012000140001600018000 -0.5

0.5

含回声信号波形

00.51 1.52 2.5

x 104 0

200

400

含回声信号幅值

00.51 1.52 2.5

x 104

-0.5

0.5

滤波后信号波形

00.51 1.52 2.5

x 104 0

100

200

滤波后信号的幅值

图5 滤除多个回波后的波形

3.4. 2设计巴特沃兹滤波器及滤波

3.4.2.1设计巴特沃斯数字低通滤波器

fp=2100;

fs=8000;

Fs=20000;

Rp=0.5;

Rs=30;

T=1/Fs; %设计指标

W1p=fp/Fs*2;

W1s=fs/Fs*2;%求归一化频率

[N,Wn]=buttord(W1p,W1s,Rp,Rs,'s'); %确定butterworth的最小介数N和频率参数Wn

[z,p,k]=buttap(N); %设计模拟低通原型的零极点增益参数 [bp,ap]=zp2tf(z,p,k); %将零极点增益转换成分子分母参数 [bs,as]=lp2lp(bp,ap,Wn*pi*Fs);%将低通原型转换为模拟低通

[bz,az]=impinvar(bs,as,Fs); %用脉冲响应不变法进行模数变换 sys=tf(bz,az,T);%给出传输函数H(Z)

[H,W]=freqz(bz,az,512,Fs); %生成频率响应参数 subplot(2,1,1);

plot(W,20*log10(abs(H))); %绘制幅频响应 grid on; %加坐标网格 xlabel('频率/Hz'); ylabel('振幅/dB'); subplot(2,1,2); plot(W,abs(H)); grid on;

xlabel('频率/Hz');

ylabel('振幅/H')；

运行后的结果所得波形如图6所示：

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

-40

-30-20-100

10振幅/d B

01000200030004000

5000600070008000900010000

0.5

1

1.5

频率/Hz

振幅/H

图6 巴特沃兹滤波器图

运行结果如下： N =4

bz = 0.0000 0.0999 0.1914 0.0252

az= 1.0000 -1.4336 1.0984 -0.4115 0.0627 可以得出：只需编程，结果非常直观

3．5估算距离

3.5.1通过理论计算法

（1） %计算单个回声方法

m=3200，而采样频率为fs=25050Hz ，可得延时t=m/fs=0.145s 设反射物距离为s ，声速为v ，则有2*s=v*t 得：s=v*t/2；

声速大概在340m/s左右，把t=0.13s,v=340m/s代入求得s=24.65m，即反射物的距离为24.65米。

3.5.2程序返回测量法

（2）采用相关分析法从带有多个回声的声音信号中估计反射物的距离。

%相关函数法测障碍物距离

运行后的结果所得波形如图7所示：

h=xcorr(y);

figure(4); %建立图形窗口

subplot(4,1,1) %4行1列第1个图

plot(abs(h)); %画出abs（h）

title('最值'); %绘图区标题

h1=h(17800:18000); %h1矩阵

[r1,t1]=max(h1'); %最大值

t1=t1+17800-1;

subplot(4,1,2); %4行1列第2个图

plot(h1) %画出h1

title('点1'); %绘图区标题

h2=h(21000:22000); %h2矩阵

[r2,t2]=max(h2); %最大值

t2=t2+21000-t1;

subplot(4,1,3); %4行1列第3个图

plot(h2); %画出h3

title('点2'); %绘图区标题

h3=h(23000:23200); %h3矩阵

[r3,t3]=max(h3); %最大值

t3=t3+23000-t1;

subplot(4,1,4); %4行1列第4个图

plot(h3); %画出h3title('点3'); %绘图区标题

t=[t2,t3];

00.51 1.52 2.53 3.5

4

x 10

4

50最值

50

100

150

200

250

300

350

400

450

-200

20点1

50

100

150

200

250

300

350

400

450

-200

20点2

50

100

150

200

250

300

350

400

450

-200

20

根据程序返回的结果，我们发现找到了三个延时的值（3200，5200)2个采样点，这和

我们起初产生回波时的2个回波的延时是一致的。假如我们得到延时的采样点个数为n ，那么我们就可以用t=n/fs 算出延时的时间，再用L=v*t 就可以算出距离，即L=v*n/fs ，最后根据实际中反射物，声源，接收器三者的相对位置就可以确定反射物与声源之间的距离了。

此次设计所选用的音频信号fs=22050Hz,故得声源距二个反射物的距离分别为： L1=340*3200/22050/2=24.7m L2=340*5200/22050/2=40.0m

4总结

这次的课程设计的大体经过：采集一个语音信号，加入回声，再设计滤波器并恢复原

信号，估算出反射物的距离，在误差允许的范围内，基本上可以实现此功能，通过这次课程设计的学习，我对于MATLAB 的强大功能有了初步的了解，同时也熟悉了如何用MATLAB 进行编程来解决一些声音信号的相关问题。这次课程设计的难点那就是后面设计滤波器滤波以及估算反射物的距离，这让我下了不少功夫，同时也非常的感谢我的指导老师，因此我也在该过程中学的了不少知识，同时我相信进一步加强对MATLAB 的学习与研究对我今后的学习将会起到很大的帮助。

5、参考文献

（1）李敏编，《MATLAB与信号与系统实验指导书》，2013；

（2）梁虹等，《信号与系统分析及MATLAB实现》，电子工业出版社，2002

（3）陈生谭等，《信号与系统》第三版，西安电子科技大学出版社，2003

（4）刘树棠译《数字信号处理——使用MATLAB》西安:西安交通大学出版社，2002 （5）程佩青．数字信号处理教程［M］．北京：清华大学出版社，2002．