MATLAB语音信号采集与处理.doc

MATLAB课程设计报告课题：语音信号采集与处理目录一、实践目的3二、实践原理：3三、课题要求：3四、MATLAB仿真41、频谱分析：42、调制与解调：53、信号变化：8快放：8慢放：8倒放：8回声：9男女变声：94、信号加噪115、用窗函数法设计FIR滤波器12FIR低通滤波器：13FIR高通滤波器：14FIR带通滤波：15一、实践目的本次课程设计的课题为基于MATLAB的语音信号采集与处理，学会运用MATLAB的信号处理功能，采集语音信号，并对语音信号进行滤波及变换处理，观察其时域和频域特性，加深对信号处理理论的理解，并为今后熟练使用MATLAB进行系统的分析仿真和设计奠定基础。此次实习课程主要是为了进一步熟悉对matlab软件的使用，以及学会利用matlab对声音信号这种实际问题进行处理，将理论应用于实际，加深对它的理解。二、实践原理：利用MATLAB对语音信号进行分析和处理，采集语音信号后，利用MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。语音信号的“ 短时谱”对于非平稳信号, 它是非周期的, 频谱随时间连续变化, 因此由傅里叶变换得到的频谱无法获知其在各个时刻的频谱特性。如果利用加窗的方法从语音流中取出其中一个短断, 再进行傅里叶变换, 就可以得到该语音的短时谱。三、课题要求：利用windows 自带的录音机或者其它录音软件，录制几段语音信号（要有几种不同的声音，要有男声、女声）。对录制的语音信号进行频谱分析，确定该段语音的主要频率范围，由此频率范围判断该段语音信号的特点（低沉or 尖锐）。利用采样定理，对该段语音信号进行采样，观察不同采样频率（过采样、欠采样、临界采样）对信号的影响。对采集到的语音信号进行调制与解调，观测调制与解调前后信号的变化。实现语音信号的快放、慢放、倒放、回声、男女变声。对语音信号加噪，然后进行滤波，分析不同的滤波方式对信号的影响。利用MATLAB GUI 制作语音信号采集与分析演示系统。四、MATLAB仿真1、频谱分析：用WINDOWS下的录音机，用单声道录制一段音乐或声音，时间在5S内。然后MATLAB软件平台下，利用函数WAVREAD对语音信号进行采样，记住采样频率和采样点数。对语音信号进行快速傅立叶变换，在一个窗口同时画出信号的时域波形图和频谱图，分析语音信号的频谱特点程序：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);plot(t/2,x) %做原始语音信号的时域图形title(原始信号时域波形图);subplot(2,1,2);plot(f,abs(y2);title(原始信号频谱图)仿真波形：门铃：和弦：男女声： 2、调制与解调：首先画出语音信号的时域波形，然后对语音信号进行频谱分析。在Matlab中可以利用函数fft对信号行快速傅里叶变换,得到信号的频谱特性，从而加深对频谱特性的理解。程序：clear; dt=1/44100; fs=44100; f1,fs,nbits=wavread(D:1huan.wav); figure(1); subplot(1,1,1); N=length(f1); t=0:1/fs:(N-1)/fs;plot(t,f1); title(信息信号的时域波形);fy1=fft(f1);w1=0:fs/(N-1):fs;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(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=5000; 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);title(滤波后的解调信号时域波形);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=wavread(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); %对三次回声信号做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调用函数：function Y=voice(x,f) %更改采样率使基频改变 f1降低;f1升高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-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 1降调,1升调4、信号加噪程序：fs=22050; %语音信号采样频率为22050x1=wavread(3yang.wav); %读取语音信号数据赋值给x1f=fs*(0:511)/1024; %将0到511，步长为1的序列的值与fs相乘并除以1024的值赋给ft=0:1/fs:(length(x1)-1)/fs; %将0到x1的长度减1后的值除以fs的值，且步长为1/fs的值，的序列的值，赋予tAu=0.05; %噪声幅值d=Au*cos(2*pi*2000*t); %干扰信号构建命令函数，构建了一个余弦函数x2=x1(:,1)+d; %取原始语音信号的单声部信号然后和噪声信号相加wavwrite(x2,22050,jiazaoyang);%生成wav文件sound(x2,22050); %播放语音信号y1=fft(x1,1024); %对信号做1024点的FFT变换y2=fft(x2,1024);figure(1); %创建图形窗subplot(2,1,1);plot(t,x1); %做原始语音信号的时域波形title(加噪前的信号);xlabel(timen); %x轴的名字是timeylabel(fuzhin); %y轴的名字是fuzhisubplot(2,1,2); %创建两行一列绘图区间的第1个绘图区间plot(t,x2)title(加噪后的信号);xlabel(timen);ylabel(fuzhin);figure(2)subplot(2,1,1);plot(f,abs(y1(1:512);title(原始语音信号频谱);xlabel(Hz);ylabel(fuzhi);subplot(2,1,2);plot(f,abs(y2(1:512);title(加噪后的信号频谱);xlabel(Hz);ylabel(fuzhi);仿真波形： 5、用窗函数法设计FIR滤波器根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为w，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出： 在一般情况下，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子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。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.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);%对信号进行低通滤波,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高通滤波器：程序：x1,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高通滤波器(时域);subplot(212);plot(-Fs/2:derta_Fs: