江大MATLAB语音信号的采集与处理课设

课程名称： Matlab课程设计 . 设计名称： 语音信号的采集和处理 . 小组成员：朱玲、 相国俊、张杰、周毅 、胥杰 .

班 级： J电信1301 .

指导教师： 毛彦欣 . 起止日期： 2015 1.11-1.17 .

题目一：基于MATLAB的语音信号的采集和处理

一、实践的目的和要求

本次课程设计的课题为《基于MATLAB的语音信号采集与处理》，学会运用MATLAB的信号处理功能，采集语音信号，并对语音信号进行滤波及变换处理，观察其时域和频域特性，加深对信号处理理论的理解，并为今后熟练使用MATLAB进行系统的分析仿真和设计奠定基础。

此次实习课程主要是为了进一步熟悉对matlab软件的使用，以及学会利用matlab对声音信号这种实际问题进行处理，将理论应用于实际，加深对它的理解。

二．本课题的设计要求及设计方案概

1：使用wavrecord录入自己的语音信号，使用save函数进行保存后使用wavplay 函数进行播放。

2：使用plot再画出该语音信号的时域波形，对原始波形进行用fft函数傅里叶

变换后，使用plot画出其频谱。

3：用窗函数法和双线性法分别设计低通，高通，带通滤波器对原始信号进行滤波。

4： 画出滤波后的信号时域、频域图 三．设计过程

1．使用wavrecord录入自己的语音信号.简单音乐信号和和弦音乐信号。并且分析各自语音信号的时域的不同点并且用fft函数对信号进行快速傅立叶变换并进行频谱分析，实现程序如下：

[x,fs,Nbits] =wavread('h15.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.语音信号的时域处理

（1）对信号进行时域尺度变换，抽取于插值,实现语音信号的快放,慢放，倒放；在matlab中，可以用sound函数对声音进行回放。实现

程序如下：

[x,fs,Nbits] =wavread('h15.wav') ; t=0:1/fs:(length(x)-1)/fs; a=1

sound(x,a*fs) ; %对加载的语音信号进行回放 pause(3) a=2

sound(x,a*fs) ; %对加载的语音信号进行快放 pause(3) a=0.5

sound(x,a*fs) ; %对加载的语音信号进行慢放 倒放程序实现如下：

[y,fs,nbits]=wavread('h15.wav'); t=0:1/fs:(length(y)-1)/fs; a=flipud(y); sound(a,fs);

(3)实现语音信号的调制与解调，实现程序如: clear; dt=1/44100; fs=44100;

[f1,fs,nbits]=wavread('h15.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);

频谱波形如下：

4理解傅里叶变换的性质

（1）对比信号时域的尺度变换，抽取与插值变换，观察其频域中频谱的变化，回放语音信号，体会时域语音信号的变化。 （2）信号的调制与解调

语音信号与高频正弦载波相调制，比较其频谱变化，回放信号，比较时域中语音信号的变化，将调制后的信号进行调解，回放信号，比较时域中语音的变化

解调

从高频已调信号中恢复出调制信号的过程称为解调(demodulation )，又称为检波(detection )。对于振幅调制信号，解调(demodulation )就是从它的幅度变化上提取调制信号的过程。解调(demodulation )是调制的逆过程。

可利用乘积型同步检波器实现振幅的解调，让已调信号与本地恢复载波信号相乘并通过低通滤波可获得解调信号。 程序

5.设计数字滤波器； 给定滤波器的性能指标如下： （1）低通滤波器的性能指标：

fb=1000Hz,fc=1200Hz,As=100dB,Ap=1dB.

(2)高通滤波器的性能指标：fb=4800Hz,fc=5000Hz,As=100dB,Ap=1dB. (3)带通通滤波器的性能指标：

fb1=1200Hz,fb2=3000Hz,fc1=1000Hz,fc2=3200Hz,As=100dB,Ap=1dB. 采用窗函数法和线性变换法设计上面要求的3种滤波器，并画出滤波器的频率响应；

采用窗函数法和双线性变换法设计上面要求的3种滤波器，并画出滤波器的频率响应。

窗函数带通滤波器程序： wn=kaiser(49);

fc1=1000;fc2=3200;fs=8000; wc1=2*fc1/fs;wc2=2*fc2/fs;

b=fir1(48,[wc1 wc2],wn); freqz(b,1);

title('窗函数带通滤波器频率响应图 ') 窗函数带通滤波器频率响应图：

窗函数高通滤波器程序： wn=kaiser(49);

fc=4800;fs=22050;wc=2*fc/fs; b=fir1(48,wc,'high',wn);freqz(b,1); figure(1); freqz(b);

title('窗函数高通滤波器频率响应图 ') 窗函数高通滤波器频率响应图：

窗函数低通滤波器程序： wn=kaiser(49);

fc=1200;fs=8000;wc=2*fc/fs; b=fir1(48,wc,wn);freqz(b,1) freqz(b);

title('窗函数低通滤波器频率响应图'); 窗函数低通滤波器频率响应图：

双线性变换法带通滤波器程序：

wp1=1200/8000*2*pi;wp2=3000/8000*2*pi; ws1=1000/8000*2*pi;ws2=3200/8000*2*pi; Rp=1;Rs=100;

Wp1=tan(wp1/2);Wp2=tan(wp2/2); Ws1=tan(ws1/2);Ws2=tan(ws2/2);

BW=Wp2-Wp1;W0=Wp1*Wp2;W00=sqrt(W0); WP=1,WS=WP*(W0^2-Ws1^2)/(Ws1*BW); [N,Wn]=cheb1ord(WP,WS,Rp,Rs,'s'); [B,A]=cheby1(N,Rp,Wn,'s'); [BT,AT]=lP2BP(B,A,W00,BW); [num,den]=bilinear(BT,AT,0.5); [h,omega]=freqz(num,den,32); subplot(2,2,1);stem(omega/pi,abs(h)); xlabel('\\omega/\\pi');ylabel('|H(z)|');

subplot(2,2,2);stem(omega/pi,20*log10(abs(h))); xlabel('\\omega/\\pi');ylabel('增益,db');

双线性变换法带通滤波器频率响应图：

双线性变换法高通滤波器程序： FS=8000

F1=2800;FH=3000; wp=(FH/FS)*2*pi; ws=(F1/FS)*2*pi; OmegaP=2*FS*tan(WP/2); OmegaS=2*FS*tan(WS/2);

[N,Wn]=buttord(OmegaP,OmegaS,1,100,'s'); [b,a]=butter(N,Wn,'s'); [bz,az]=bilinear(b,a,FS); freqz(bz,az,4096,FS,'whole'); title('双线性变换法高通')

双线性变换法高通滤波器频率波形图：

双线性变换法低通滤波器程序： FS=8000

F1=1000;Fh=1200; WP=(F1/FS)*2*pi; WS=(Fh/FS)*2*pi; OmegaP=2*FS*tan(WP/2); OmegaS=2*FS*tan(WS/2);

[n,Wn]=buttord(OmegaP,OmegaS,1,100,'s'); [b,a]=butter(n,Wn,'s'); [bz,az]=bilinear(b,a,FS); freqz(bz,az,4096,FS,'whole'); title('双线性变换法低通滤波器') 双线性变换法低通滤波器频率波形图：

6. 用滤波器对信号进行滤波

然后用自己设计的滤波器对采集到的信号进行滤波，画出滤波后信号的时域波形及频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号，分析滤波前后的语音变化；

窗函数低通滤波器滤波实现程序： [x1,Fs,bits]=wavread('h15.wav');

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,'low',Window); %wc1/pi为归一化,窗函数法设计函数 figure(2); freqz(b1,1,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低通滤波器（频域）');

窗函数高通滤波器滤波的程序： clc; clear; close all;

[x1,Fs,bits]=wavread('h15.wav');

derta_Fs = Fs/length(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致 fs=Fs; As2=100; fp2=3000; fs2=2700; 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(2,1,1); plot(x1_high);

title('信号经过FIR高通滤波器(时域)'); subplot(2,1,2);

plot([-Fs/2:derta_Fs: Fs/2-derta_Fs],abs(fftshift(fft(x1_high)))); title('信号经过FIR高通滤波器（频域）');

窗函数带通滤波器滤波的程序： clc; clear; close all;

[x1,Fs,bits]=wavread('h15.wav');

derta_Fs = Fs/length(x1);%设置频谱的间隔，分辨率 ，这里保证了x轴的点数必须和y轴点数一致 fs=Fs; As3=100;

fp3=[1200,3000];fs3=[1000,3200];

wp3=2*pi*fp3/fs; ws3=2*pi*fs3/fs; BF3=wp3(1)-ws3(1); wc3=wp3+BF3/2;

M3=ceil((As3-7.95)/(2.286*BF3))+1;%按凯泽窗计算滤波器阶数 N3=M3+1;

beta3=0.1102*(As3-8.7);

Window=(kaiser(N3,beta3)); %求凯泽窗窗函数 b3=fir1(M3,wc3/pi,'bandpass',Window);%带通滤波器 figure(6);

freqz(b3,1,512);%数字滤波器频率响应 title('FIR带通滤波器的频率响应');

x1_daitong = filter(b3,1,x1);%对信号进行带通滤波 sound(x1_daitong,Fs,bits); figure(7); subplot(2,1,1); plot(x1_daitong);

title('信号经过FIR带通滤波器(时域)'); subplot(2,1,2);

plot([-Fs/2:derta_Fs: Fs/2-derta_Fs],abs(fftshift(fft(x1_daitong)))); title('信号经过FIR带通滤波器（频域）');

双线性法低通滤波器滤波后的程序：

clear all; Ft=22050; Fp=1000; Fs=1200;

[y,fs,nbits]=wavread('h15.wav'); %sound(y,fs,nbits); y=y-mean(y); as=100; ap=1;

wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fs1=1;

fp=2*tan(wp/2); fs=2*tan(ws/2);

[n11,wn11]=ellipord(wp,ws,ap,as,'s'); [b11,a11]=ellip(n11,ap,as,wn11,'low','s'); [num11,den11]=bilinear(b11,a11,fs1); [h,w]=freqz(num11,den11,512,Ft); subplot(2,1,1);

plot( w,20*log10(abs(h)) ); %画对数幅度谱 hold on;

xlabel( '归一化频率w/pi' ); ylabel( '幅度(dB)' ); title( 'IIR-幅度响应');

subplot(2,1,2);

plot( w,angle(h) ); hold on;

xlabel( '归一化频率w/pi' ); ylabel( '相位' ); title( 'IIR-相位响应');

f1=filter(num11,den11,y); figure(2)

subplot(2,1,1)

plot(y) %画出滤波前的时域图 title('IIR低通滤波器滤波前的时域波形'); subplot(2,1,2)

plot(f1); %画出滤波后的时域图 title('IIR低通滤波器滤波后的时域波形');

sound(f1); %播放滤波后的信号 F0=fftshift(abs(fft(f1)));

w = linspace(-Ft/2,Ft/2,length(F0)); figure(3)

y2=fftshift(abs(fft(y)));

subplot(2,1,1);

plot(w,y2); %画出滤波前的频谱图 title('IIR低通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值');

w = linspace(-Ft/2,Ft/2,length(y2)); subplot(2,1,2)

plot(w,abs(F0)); %画出滤波后的频谱图 title('IIR低通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');

双线性法高通滤波器滤波的程序：

clear all; Ft=22050; Fp=5000; Fs=4800;

[y,fs,nbits]=wavread('h15.wav'); %sound(y,fs,nbits); y=y-mean(y); as=100; ap=1;

wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fs1=1;

fp=2*tan(wp/2); fs=2*tan(ws/2);

[n11,wn11]=ellipord(wp,ws,ap,as,'s'); [b11,a11]=ellip(n11,ap,as,wn11,'high','s'); [num11,den11]=bilinear(b11,a11,fs1); [h,w]=freqz(num11,den11,512,Ft); subplot(2,1,1);

plot( w,20*log10(abs(h)) ); %画对数幅度谱 hold on;

xlabel( '归一化频率w/pi' ); ylabel( '幅度(dB)' ); title( 'IIR-幅度响应');

subplot(2,1,2);

plot( w,angle(h) ); hold on;

xlabel( '归一化频率w/pi' ); ylabel( '相位' ); title( 'IIR-相位响应');

f1=filter(num11,den11,y); figure(2)

subplot(2,1,1)

plot(y) %画出滤波前的时域图 title('IIR高通滤波器滤波前的时域波形'); subplot(2,1,2)

plot(f1); %画出滤波后的时域图 title('IIR高通滤波器滤波后的时域波形');

sound(f1); %播放滤波后的信号 F0=fftshift(abs(fft(f1)));

w = linspace(-Ft/2,Ft/2,length(F0)); figure(3)

y2=fftshift(abs(fft(y))); subplot(2,1,1);

plot(w,y2); %画出滤波前的频谱图 title('IIR高通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值');

w = linspace(-Ft/2,Ft/2,length(y2)); subplot(2,1,2)

plot(w,abs(F0)); %画出滤波后的频谱图 title('IIR高通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');

双线性法设计带通滤波器滤波的程序：

clear all; Ft=22050; Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200;

[y,fs,nbits]=wavread('h15.wav'); %sound(y,fs,nbits); y=y-mean(y); Fp=[Fp1,Fp2]; Fs=[Fs1,Fs2]; as=100; ap=1; T=2;

wp1=2*pi*Fp1/Ft; wp2=2*pi*Fp2/Ft; ws1=2*pi*Fs1/Ft; ws2=2*pi*Fs2/Ft;

Wp1=(2/T)*tan(wp1/2); Wp2=(2/T)*tan(wp2/2);

Ws1=(2/T)*tan(ws1/2); Ws2=(2/T)*tan(ws1/2)

W0=Wp1*Wp2; w0=sqrt(W0);

BW=Wp2-Wp1; %带通滤波器的通带宽度 lp=1; %归一化处理

ls=Ws1*BW/(W0-Ws1^2); [N,Wn]=ellipord(lp,ls,ap,as,'s'); [B,A]=ellip(N,1,40,Wn,'s'); [BT,AT]=lp2bp(B,A,w0,BW); [num,den]=bilinear(BT,AT,0.5); [z,p,k]=tf2zp(num,den);

[h,w]=freqz(num,den,512,Ft); figure(1)

subplot(2,1,1);

plot( w,20*log10(abs(h)) ); %画对数幅度谱 hold on;

xlabel( '频率w/pi' ); ylabel( '幅度(dB)' ); title( 'IIR-幅度响应'); subplot(2,1,2);

plot( w,angle(h) ); hold on;

xlabel( '频率w/pi' ); ylabel( '相位' ); title( 'IIR-相位响应');

f1=filter(num,den,y); figure(2)

subplot(2,1,1)

plot(y) %画出滤波前的时域图 title('IIR带通滤波器滤波前的时域波形'); subplot(2,1,2)

plot(f1); %画出滤波后的时域图 title('IIR带通滤波器滤波后的时域波形');

sound(f1); %播放滤波后的信号 F0=fftshift(abs(fft(f1)));

w = linspace(-Ft/2,Ft/2,length(F0)); figure(3)

y2=fftshift(abs(fft(y))); subplot(2,1,1);

plot(w,y2); %画出滤波前的频谱图 title('IIR带通滤波器滤波前的频谱') xlabel('频率/Hz'); ylabel('幅值');

w = linspace(-Ft/2,Ft/2,length(y2)); subplot(2,1,2)

plot(w,abs(F0)); %画出滤波后的频谱图 title('IIR带通滤波器滤波后的频谱') xlabel('频率/Hz'); ylabel('幅值');

7. 设计系统界面

为了使编制的程序操作方便，设计处理系统的用户界面，在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数、显示滤波器的频率响应，选择信号等。

题目六：语音信号变声系统处理

电视台经常针对某些事件的知情者进行采访，为了保护知情者，

经常改变说话人的声音，请利用所学的知识，将其实现。

要求处理后的语音信号基本不受影响正常收听与理解； 对处理参数能够通过matlab界面进行调节，以对比不同处理结果；

语音信号变声系统程序：

[y,fs,nbits]=wavread('D:\\yuyin\\hc'); %读取声音文件 x=y(:,1); %读入的y矩阵有两列,取第1列

sound(voice(x,1.4),fs,nbits); %调整voice()第2个参数转换音调,>1降调,<1升调

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-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; else

ysim = Y(ypos + simix); rxy = zeros(1, Kmax+1); rxx = zeros(1, Kmax+1); Kmin = 0;

for k = Kmin:kdecim:Kmax

xsim = padX(Wsim + xpos + k + simix); rxx(k+1) = norm(xsim); rxy(k+1) = (ysim * xsim');

end

Rxy = (rxx ~= 0).*rxy./(rxx+(rxx==0)); km = min(find(Rxy == max(Rxy))-1); end

xabs = xpos+km; Y(ypos+ovix)

=

((1-xfwin).*Y(ypos+ovix))

+

(xfwin.*padX(Wsim+xabs+ovix));

Y(ypos+newix) = padX(Wsim+xabs+newix); end End

心得体会：

通过这次MATLAB的学习，我对MATLAB有了一个基础的认识，matlab

是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图示和精确计算于一身，在应用数学、物理、化工、机电工程、医药、金融和其他需要进行复杂数值计算的领域得到了广泛应用。MATLAB是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程的特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂程序（M文件）后再一起运行。 在这短短的一周内从开始的一头雾水，到自己看书学习，到同学讨论，再进行整个题目的理论分析和计算，参考课程上的代码，写出自己的代码。

我们也明白了学无止尽的道理，在我们所查的很多参考书中，很多知识是我们从没有接触过的，我们对它的了解还仅限于皮毛，对它的很多功能以及函数还不是很了解，所以在这个学习的过程中我们穿越在知识的海洋中，一点一点吸取着它的知识。在MATLAB编程中需要很多的参考书，要尽量多的熟悉matlab自带的函数及其作用，因为matlab的自带函数特别多，基本上能够满足一般的数据和矩阵的计算，所以基本上不用你自己编函数。这一点对程序非常有帮助，可以使程序简单，运行效率高，可以节省很多时间。本次课设中用了很多MATLAB自带的函数，使程序变得很简单。 把基本的知识看过之后，就需要找一个实际的程序来动手编一下，不要等所有的知识都学好之后再去编程，你要在编程的过程中学习，程序需要什么知识再去补充，编程是一点一点积累的，所以你要需做一些随手笔记什么的。 在编写程序代码时，需要什么函数，需要什么模块就应该去着重看那个知

识点，不要一步登天，一步一步学，如果太急于把所有东西都学到，也是不好的，更是实现不了的。所以那时一天一天积累的，慢慢地学通这个软件。 总之，通过这次学习，我了解了一下这个软件总体的功能，以及通过自己编写一些代码也学到了一些用法和知识。更了解到了，我们还有好多东西去学，学无止尽。

本文来源：https://www.bwwdw.com/article/w5b6.html