语音信号的滤波与频谱分析

1、生物医学信号处理大作业题目：语音信号的滤波与频谱分析学生姓名 学院名称 精密仪器与光电子工程 专 业 学 号 一、实验目的语音信号的滤波与频谱分析录制自己的一段语音：“天津大学精密仪器与光电子工程学院生物医学工程X班XXX, College of precision instrument and opto-electronics engineering, biomedical engineering”，时间控制在15秒到30秒左右；利用wavread函数读入语言信号，记住采样频率。二、实验过程（1）求原始语音信号的特征频带：可以分别对一定时间间隔内，求功率谱（傅里叶变换结果取模的平方）并画出功

2、率谱。（2）根据语音信号频谱特点，设计FIR或IIR滤波器，分别画出滤波器幅频和相频特性曲线。说明滤波器特性参数。用设计的滤波器对信号滤波，画出滤波后时域波形。用sound函数回放语音信号。（3）求出特征频段语音信号随时间变化的曲线（每隔0.05秒求一次功率谱，连接成曲线）。 （4）选做：语谱图：横轴为时间，纵轴为频率，灰度值大小表示功率谱值的大小。(提示，可以采用spectrogram函数)（1）读入语音文件并画出其时域波形和频域波形，实现加窗fft并求出其功率谱。clcclear all;close all;x,Fs,bits=wavread('C:Users刘冰Desktop数字

3、信号处理liubing');x0=x(:,1);%将采集来的语音信号转换为一个数组sound(x0,Fs,bits);y=fft(x);figure;plot(x,b);title('原始语音信号时域波形');y1=fft(x0); y1=fftshift(y1); d = Fs/length(x); figure;plot(-Fs/2:d: Fs/2-d,abs(y1),b);title('原始语音信号的频域信号');% 画出原始语音信号的频谱图N=length(x);w1 = window(hann,N);w2 = window(blackman,N

4、);x1=x0.*w1;%对原始信号加汉宁窗处理x2=x0.*w2；%对原始信号加布兰克曼窗处理figure,plot(x1);title(加汉宁窗后的语音信号)%显示加窗后的时域语音信号s=floor(length(x0)/Fs);%计算原始语音信号的时间长度，这里得到的结果是18秒，因为floor是向下取整，所以信号的末尾一点会被去掉，但是因为最后一点没有声音信号，所以影响可以忽略。%加汉宁窗后功率谱，加布兰克曼窗后又可以得到一组图，只需要将下列循环中的x1改为x2，这里就不再显示%每两秒对语音信号求一次功率谱并显示for i=1:1:s/2 f=x1(i-1)*Fs*2+1:i*Fs*2

5、);%每两秒取出一段信号 l=length(f); q=fft(f,l); E=abs(q).*abs(q); %傅里叶变换结果取模的平方 figure,plot(E(1:3000),'b');title('第',num2str(i*2-1),'',num2str(i*2),'秒语音功率谱');%因为语音信号主要集中在低频段，所以这里只需要显示低频段即可，取（1:3000）end加窗后的时域波形频谱分析：人说话的频率基本集中在1200Hz以内的低频段，但我这里可以在高频段（2000Hz左右）可以观察到部分能量，这应该是在录音的时候

6、电脑本身的噪音以及录音设备的误差造成（2）根据功率谱线，大致可观察到语音信号在150400和500650之间有两个波峰，因此取这两个频宽作为我语音信号的特征频带，根据要求选择汉宁窗作为滤波器滤波。1）汉宁窗滤波器的幅频特性显示%汉宁窗 f1=150;%通带上下两个频率f2=400; %带通滤波器的通带范围w1=2*pi*f1/Fs; %归一化 w2=2*pi*f2/Fs; w=w1,w2;N=ceil(4*pi/(2*pi*200/Fs);%求出所需滤波器的阶数Windows=HANNING(N);b1=FIR1(N-1,w, Windows);%带通滤波器figure;freqz(b1,1,

7、512);title('汉宁窗带通（150400）滤波器的频率响应');%数字滤波器频率响应f1=500;f2=650; %带通滤波器的通带范围w1=2*pi*f1/Fs; w2=2*pi*f2/Fs; w=w1,w2;N=ceil(4*pi/(2*pi*200/Fs);Nw=N;Windows=HANNING(N);b2=FIR1(N-1,w, Windows);%带通滤波器figure;freqz(b2,1,512);title('汉宁窗带通（500650）滤波器的频率响应');%数字滤波器频率响应2）将原始语音信号分别通过两个频带，以下是滤波后的时域波形b

8、and1=fftfilt(b1,x0);sound(band1,Fs,bits);%滤波后听到的声音有了明显的变化subplot(2,1,1);plot(band1);title('汉宁窗（150400）滤波后的时域波形');band2=fftfilt(b2,x0);sound(band2,Fs,bits);subplot(2,1,2);plot(band2);title('汉宁窗（500650）滤波后的时域波形');（3）特征频段语音信号随时间变化的曲线，分别作出150400和500650特征频带的语音信号曲线a=zeros(s/0.05,1);%这一步是为了

9、求出所画信号曲线的长度，并定义空数组for i=1:s/0.05 f=x0(i-1)*Fs*0.05+1:i*Fs*0.05);%每0.05秒取出一段数据 l=length(f); q=fft(f,l); E=q.*conj(q)/l;%求出该段的功率 a(i)=mean(E(150400);%求平均endfigure,plot(0:0.05:s-0.05,a);title('150400Hz汉宁窗滤波后的功率谱曲线')b=zeros(s/0.05,1);for i=1:s/0.05 f=x0(i-1)*Fs*0.05+1:i*Fs*0.05); l=length(f); q=

10、fft(f,l); E=q.*conj(q)/l; b(i)=mean(E(500:650);endfigure,plot(0:0.05:s-0.05,b);title('500650Hz汉宁窗滤波后的功率谱曲线')(4)语谱图，也叫频谱分析视图，如果针对语音数据的话，叫语谱图。语谱图的横坐标是时间，纵坐标是频率，坐标点值为语音数据能量。由于是采用二维平面表达三维信息，所以能量值的大小是通过颜色来表示的，颜色深，表示该点的语音能量越强。所有函数：spectrogram功能：使用短时傅里叶变换得到信号的频谱图。语法：S,F,T,P=spectrogram(x,window,noverlap,nfft,fs)x-输入信号的向量。默认情况下，即没有后续输入参数，x将被分成8段分别做变换处理，如果x不能被平分成8段，则会做截断处理。默认情况下，其他参数的默认值为window-窗函数，默认为nfft长度的海明窗Hammingnoverlap-每一段的重叠样本数，默认值是在各段之间产生50%的重叠nfft-做FFT变换的长度，默认为256和大于每段长度的最小2次幂之