MATLAB语音信号采集与处理(DOC)

1、MATLAB课程设计报告课题：语音信号采集与处理MATLAB 仿真 4.1、频谱分析：4.2、调制与解调：5.3、信号变化：8.快放： 8.慢放： 8.倒放： 8.回声： 8.男女变声：9.4、信号加噪 1.0.5、用窗函数法设计FIR 滤波器 1.1FIR低通滤波器：1.2FIR高通滤波器：1.3FIR带通滤波：1.4本次课程设计的课题为基于MATLAB的语音信号采集与处理，学会运用MATLAB的信号处理功能，采集语音信号，并对语音信号进行滤波及变换处理，观察其时域和频域特性，加深对信号处理理论的理解，并为今后熟练使用MATLAB进行系统的分析仿真和设计奠定基础。此次实习课程主要是为了进一步

2、熟悉对matlab 软件的使用，以及学会利用matlab 对声音信号这种实际问题进行处理，将理论应用于实际，加深对它的理解。二、实践原理：利用MATLAB对语音信号进行分析和处理，采集语音信号后，利用MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。语音信号的“短时谱”对于非平稳信号 , 它是非周期的 , 频谱随时间连续变化, 因此由傅里叶变换得到的频谱无法获知其在各个时刻的频谱特性。 如果利用加窗的方法从语音流中取出其中一个短断, 再进行傅里叶变换, 就可以得到该语音的短时谱。三、课题要求：1 利用wind

3、ows 自带的录音机或者其它录音软件，录制几段语音信号（要有几种不同的声音，要有男声、女声）。2 对录制的语音信号进行频谱分析，确定该段语音的主要频率范围，由此频率范围判断该段语音信号的特点（低沉or 尖锐）。3 利用采样定理，对该段语音信号进行采样，观察不同采样频率（过采样、欠采样、临界采样）对信号的影响。4 对采集到的语音信号进行调制与解调，观测调制与解调前后信号的变化。5 实现语音信号的快放、慢放、倒放、回声、男女变声。6 对语音信号加噪，然后进行滤波，分析不同的滤波方式对信号的影响。7 利用 MATLAB GUI 制作语音信号采集与分析演示系统。MATLAB 仿真1 、频谱分析：用 W

4、INDOWS下的录音机，用单声道录制一段音乐或声音，时间在5S 内。然后MATLAB软件平台下，利用函数WAVREA对语音信号进行采样，D记住采样频率和采样点数。对语音信在一个窗口同时画出信号的时域波形图和频谱图，分析语音信号的频谱特点程序：fs =22050;Nbits =16;x,fs,Nbits =wavread('D:matlab22hexian.wav') ;%读声音文件n=length(x);t=0:1/fs:(length(x)-1)/fs;y1=fft(x,n) ;y2=fftshift(y1);f=0:fs/n:fs*(n-1)/n;subplot(2,1,1

5、);plot(t/2,x)求出语音信号的长度 傅里叶变换 对频谱图进行平移得出频点做原始语音信号的时域图形title（' 原始信号时域波形图'）;subplot(2,1,2); plot(f,abs(y2); title(' 原始信号频谱图 仿真波形：')2 和弦：3 男女声：2、调制与解调：fy1=fft(f1);w1=0:fs/(N-1):fs;首先画出语音信号的时域波形，然后对语音信号进行频谱分析。在Matlab 中可以利程序：clear;dt=1/44100;fs=44100;f1,fs,nbits=wavread('D:1huan.wav

6、9;);figure(1);subplot(1,1,1);N=length(f1);t=0:1/fs:(N-1)/fs;plot(t,f1);title(' 信息信号的时域波形');用函数 fft对信号行快速傅里叶变换, 得到信号的频谱特性，从而加深对频谱特性的理解。figure(2);subplot(1,1,1);plot(w1,abs(fy1);title(' 信息信号的频谱');f2=cos(22000*pi*t);figure(3);subplot(1,1,1);fy2 = fft(f2);N2=length(f2);w2=fs/N*0:N-1;plot

7、(w2,abs(abs(fy2);title(' 载波信号的频谱');f1=f1(:,1);f3=f1'.*f2;figure(4);subplot(1,1,1);fy3 = fft(f3);plot(w1,abs(abs(fy3);title(' 已调信号的频谱');sound(f3,fs,nbits);f4=f3.*f2;figure(5);subplot(1,1,1);fy4=fft(f4);plot(w1,abs(abs(fy4);title(' 解调信号的频谱');sound(f4,fs,nbits);fp1=0;fs1=500

8、0;As1=100;wp1=2*pi*fp1/fs;ws1=2*pi*fs1/fs;BF1=ws1-wp1;wc1=(wp1+ws1)/2;M1=ceil(As1-7.95)/(2.286*BF1)+1;N1=M1+1;beta1=0.1102*(As1-8.7);Window=(kaiser(N1,beta1);b1=fir1(M1,wc1/pi,Window);figure(6);subplot(1,1,1);freqz(b1,1,512);title('FIR 低通滤波器的频率响应');f4_low = filter(b1,1, f4);plot(t,f4_low);ti

9、tle(' 滤波后的解调信号时域波形');sound(f4_low,fs,nbits);f5=fft(f4_low);figure(7);subplot(1,1,1);plot(w1,abs(f5);title(' 滤波后的解调信号频谱');3、信号变化：快放：x,fs,nbits=wavread('D:1huan.wav')w=2；M=w*fs;wavplay(x,M);慢放：x,fs,nbits=wavread('D:1huan.wav')w=0.8M=w*fs;wavplay(x,M);倒放：x,fs,Nbits=wavre

10、ad('D:1huan.wav');y0=flipud(x);sound(y0,fs);程序：x,fs,bits=wavread('D:3yang.wav',1 40000);% 读取语音信号 n1=0:2000;b=x(:,1);%产生单声道信号N=3;yy2=filter(1,1,zeros(1,80000/(N+1),0.7,b',zeros(1,40000); figure(3)% 三次回声滤波器时域波形subplot(2,1,1);plot(yy2);title(' 三次回声滤波器时域波形');YY2=fft(yy2);%对三次

11、回声信号做FFT 变换subplot(2,1,2);plot(n1(1:1000),YY2(1:1000);%三次回声滤波器频谱图title(' 三次回声滤波器频谱图');figure(4)subplot(2,1,1);plot(abs(YY2);%经傅里叶变换之后的信号的幅值title(' 幅值');subplot(2,1,2);plot(angle(YY2);%经傅里叶变换之后的信号的相位title(' 相位');sound(2*yy2,fs,bits);% 经三次回声滤波器后的语音信号，乘以2是为了加强信号仿真波形：男女变声：Voice 调

12、用函数：function Y=voice(x,f)%更改采样率使基频改变f>1 降低 ;f<1 升高f=round(f*1000);d=resample(x,f,1000);%时长整合使语音文件恢复原来时长W=400;Wov=W/2;Kmax=W*2;Wsim=Wov;xdecim=8;kdecim=2; X=d'F=f/1000;Ss =W-Wov;xpts = size(X,2);ypts = round(xpts / F);Y = zeros(1, ypts);xfwin = (1:Wov)/(Wov+1);ovix = (1-Wov):0; newix = 1:(W

13、-Wov);simix = (1:xdecim:Wsim) - Wsim;padX = zeros(1, Wsim), X, zeros(1,Kmax+W-Wov);Y(1:Wsim) = X(1:Wsim); lastxpos = 0; km = 0;for ypos = Wsim:Ss:(ypts-W)xpos = round(F * ypos);kmpred = km + (xpos - lastxpos);lastxpos = xpos;if (kmpred <= Kmax)km = kmpred;elseysim = Y(ypos + simix);rxy = zeros(1,

14、 Kmax+1);rxx = zeros(1, Kmax+1);Kmin = 0;for k = Kmin:kdecim:Kmaxxsim = padX(Wsim + xpos + k + simix);rxx(k+1) = norm(xsim);rxy(k+1) = (ysim * xsim');endRxy = (rxx = 0).*rxy./(rxx+(rxx=0);km = min(find(Rxy = max(Rxy)-1);endxabs = xpos+km;Y(ypos+ovix) = (1-xfwin).*Y(ypos+ovix) + (xfwin.*padX(Wsim

15、+xabs+ovix);Y(ypos+newix) = padX(Wsim+xabs+newix);endend变声程序：%读取声音文件%读入的y 矩阵有两列,取第1 列%调整 voice()第 2 个参数转换音调,>1 降调 ,<1 升调y,fs,nbits=wavread('3yang.wav');x=y(:,1);sound(voice(x,1.4),fs,nbits);4、信号加噪程序：fs=22050;%x1=wavread('3yang.wav'); % f=fs*(0:511)/1024;%以 1024 的值赋给ft=0:1/fs:(l

16、ength(x1)-1)/fs;%为 1/fs 的值，的序列的值，赋予tAu=0.05;%d=Au*cos(2*pi*2000*t)'%x2=x1(:,1)+d;%语音信号采样频率为将 0 到 511，步长为22050x11 的序列的值与fs 相乘并除将 0 到 x1 的长度减1 后的值除以fs 的值，且步长噪声幅值干扰信号构建命令函数，构建了一个余弦函数 取原始语音信号的单声部信号然后和噪声信号相加生成wav文件wavwrite(x2,22050,'jiazaoyang');%sound(x2,22050);y1=fft(x1,1024);y2=fft(x2,1024

17、);figure(1);%subplot(2,1,1);plot(t,x1);%title(' 加噪前的信号');xlabel('time n'); ylabel('fuzhi n');subplot(2,1,2);%plot(t,x2) title(' 加噪后的信号');xlabel('time n'); ylabel('fuzhi n'); figure(2) subplot(2,1,1); plot(f,abs(y1(1:512);title(' 原始语音信号频谱');xlab

18、el('Hz'); ylabel('fuzhi'); subplot(2,1,2);plot(f,abs(y2(1:512); title(' 加噪后的信号频谱');%播放语音信号%对信号做1024 点的FFT变换创建图形窗做原始语音信号的时域波形%x轴的名字是time%y轴的名字是fuzhi创建两行一列绘图区间的第1 个绘图区间xlabel('Hz');ylabel('fuzhi');5、用窗函数法设计FIR 滤波器根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1） ，窗函数类型可根

19、据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为 w，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N 和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd( n) 。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n) ， 如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下

20、，hd(n) 是不能用封闭公式表示的，需要采用数值方法表示；从 w=0到 w=2采样 N 点，采用离散傅里叶反变换(IDFT)即可求出。用窗函数wd(n) 将 hd(n) 截断，并进行加权处理，得到如果要求线性相位特性，则 h(n) 还必须满足：根据上式中的正、负号和长度N 的奇偶性又将线性相位FIR 滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n) 一类，而不能选h(n)=-h(N-1-n) 一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用 freqz 子程序， 如果不满足要求，可根据具体情况，

21、调整窗函数类型或长度，直到满足要求为止。FIR 低通滤波器：程序：x1,Fs,bits=wavread('jiazapyang');derta_Fs = Fs/length(x1);% 设置频谱的间隔，分辨率， 这里保证了x轴的点数必须和y 轴点数一致fs=Fs;fp1=1000;fs1=1200;As1=100;wp1=2*pi*fp1/fs;%ws1=2*pi*fs1/fs;%BF1=ws1-wp1;wc1=(wp1+ws1)/2;M1=ceil(As1-7.95)/(2.286*BF1)+1;% 按凯泽窗计算滤波器阶数N1=M1+1;beta1=0.1102*(As1-8

22、.7);Window=(kaiser(N1,beta1);%求凯泽窗窗函数b1=fir1(M1,wc1/pi,Window); %wc1/pi 为归一化, 窗函数法设计函数figure(2);freqz(b1,1,512);%H,w=freqz(B,A,N),(1) 中 B和A分别为离散系统的系统函数分子、分母多项式的系数向量，返回量H 则包含了离散系统频响在0pi 范围内N个频率等分点的值 (其中 N为正整数)， w则包含了范围内N个频率等分点。调用默认的N时， 其值是512。title('FIR 低通滤波器的频率响应');x1_low = filter(b1,1, x1);

23、% 对信号进行低通滤波,Y = filter(B,A,X) ， 输入 X为滤波前序列，Y为滤波结果序列，B/A 提供滤波器系数，B为分子，A为分母sound(x1_low,Fs,bits);figure(3);subplot(2,1,1);plot(x1_low);title(' 信号经过FIR 低通滤波器( 时域 )');subplot(2,1,2);plot(-Fs/2:derta_Fs:Fs/2-derta_Fs,abs(fftshift(fft(x1_low);title(' 信号经过FIR 低通滤波器(频域)');仿真波形：FIR 高通滤波器：程序：x

24、1,Fs,bits=wavread('jiazaoyang.wav');derta_Fs = Fs/length(x1);% 设置频谱的间隔，分辨率， 这里保证了x轴的点数必须和y 轴点数一致fs=Fs;As2=100;fp2=4000;fs2=2000;wp2=2*pi*fp2/fs;ws2=2*pi*fs2/fs;BF2=wp2-ws2;wc2=(wp2+ws2)/2;M2=ceil(As2-7.95)/(2.286*BF2)+1;% 按凯泽窗计算滤波器阶数N2=M2+1;beta2=0.1102*(As2-8.7);Window=(kaiser(N2,beta2);%求凯泽窗窗函数b2=fir1(M2,wc2/pi,'high',Window);figure(4);freqz(b2,1,512);% 数字滤波器频率响应title('FIR 高通滤波器的频率响应');x1_high = filter(b2,1,x1);% 对信号进行高通滤波sound(x1_high,Fs,bits);figure(5);subplot(211);plot(x1_high);title(' 信号经过FIR 高通滤波器( 时域 )