生物医学信号处理大作业题目：语音信号的滤波与频谱分析学生姓名学院名称精密仪器与光电子工程专业学号一、实验目的语音信号的滤波与频谱分析录制自己的一段语音：“天津大学精密仪器与光电子工程学院生物医学工程X班XXX, College of precision instrument and opto-electronics engineering, biomedical engineering”，时间控制在15秒到30秒左右；利用wavread 函数读入语言信号，记住采样频率。

二、实验过程（1）求原始语音信号的特征频带：可以分别对一定时间间隔内，求功率谱（傅里叶变换结果取模的平方）并画出功率谱。

（2）根据语音信号频谱特点，设计FIR或IIR滤波器，分别画出滤波器幅频和相频特性曲线。

说明滤波器特性参数。

用设计的滤波器对信号滤波，画出滤波后时域波形。

用sound 函数回放语音信号。

（3）求出特征频段语音信号随时间变化的曲线（每隔0.05秒求一次功率谱，连接成曲线）。

（4）选做：语谱图：横轴为时间，纵轴为频率，灰度值大小表示功率谱值的大小。

(提示，可以采用spectrogram函数)（1）读入语音文件并画出其时域波形和频域波形，实现加窗fft 并求出其功率谱。

clcclear all; close all;[x,Fs,bits]=wavread('C:\Users\刘冰\Desktop\数字信号处理\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('原始语音信号的频域信号'); % 画出原始语音信号的频谱图123456789x 105-1-0.8-0.6-0.4-0.200.20.40.60.81原始语音信号时域波形N=length(x);w1 = window(@hann,N); w2 = window(@blackman,N); x1=x0.*w1; %对原始信号加汉宁窗处理 x2=x0.*w2； %对原始信号加布兰克曼窗处理 figure,plot(x1);title(加汉宁窗后的语音信号) %显示加窗后的时域语音信号 s=floor(length(x0)/Fs);%计算原始语音信号的时间长度，这里得到的结果是18秒，因为floor 是向下取整，所以信号的末尾一点会被去掉，但是因为最后一点没有声音信号，所以影响可以忽略。

%加汉宁窗后功率谱，加布兰克曼窗后又可以得到一组图，只需要将下列循环中的x1改为x2，这里就不再显示%每两秒对语音信号求一次功率谱并显示for i=1:1:s/2f=x1((i -1)*Fs*2+1:i*Fs*2); %每两秒取出一段信号 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-2.5-2-1.5-1-0.500.51 1.52 2.5x 104050010001500200025003000原始语音信号的频域波形加窗后的时域波形0123456789x 105-0.8-0.6-0.4-0.200.20.40.60.80500100015002000250030000500100015002000250030003500400045005000第1~2秒语音功率谱频谱分析：人说话的频率基本集中在1200Hz以内的低频段，但我这里可以在高频段（2000Hz左右）可以观察到部分能量，这应该是在录音的时候电脑本身的噪音以及录音设备的误差造成（2）根据功率谱线，大致可观察到语音信号在150~400和500~650之间有两个波峰，因此取这两个频宽作为我语音信号的特征频带，根据要求选择汉宁窗作为滤波器滤波。

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,512);title('汉宁窗带通（150~400）滤波器的频率响应');%数字滤波器频率响应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('汉宁窗带通（500~~650）滤波器的频率响应');%数字滤波器频率响应00.10.20.30.40.50.60.70.80.91-2000-100001000Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )-300-200-1000100Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )汉宁窗带通（150~400）滤波器的频率响应0.10.20.30.40.50.60.70.80.91-1500-1000-5000500Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s)00.10.20.30.40.50.60.70.80.91-300-200-1000Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )汉宁窗带通（500~~650）滤波器的频率响应2）将原始语音信号分别通过两个频带，以下是滤波后的时域波形band1=fftfilt(b1,x0);sound(band1,Fs,bits); %滤波后听到的声音有了明显的变化subplot(2,1,1);plot(band1);title('汉宁窗（150~400）滤波后的时域波形'); band2=fftfilt(b2,x0); sound(band2,Fs,bits);subplot(2,1,2);plot(band2);title('汉宁窗（500~650）滤波后的时域波形');（3）特征频段语音信号随时间变化的曲线，分别作出150~400和500~650特征频带的语音信号曲线a=zeros(s/0.05,1);%这一步是为了求出所画信号曲线的长度，并定义空数组 for i=1:s/0.05f=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(150~400)); %求平均 endfigure,plot(0:0.05:s -0.05,a);title('150~400Hz 汉宁窗滤波后的功率谱曲线')123456789x 105-2-1012汉宁窗（150~400）滤波后的时域波形123456789x 105-0.4-0.200.20.4汉宁窗（500~650）滤波后的时域波形b=zeros(s/0.05,1); for i=1:s/0.05f=x0((i -1)*Fs*0.05+1:i*Fs*0.05); l=length(f); q=fft(f,l);E=q.*conj(q)/l;b(i)=mean(E(500:650)); endfigure,plot(0:0.05:s -0.05,b);title('500~650Hz 汉宁窗滤波后的功率谱曲线')0.511.522.5-3150~400Hz 汉宁窗滤波后的功率谱曲线-5500~650Hz 汉宁窗滤波后的功率谱曲线(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次幂之间的最大值。

fs---采样频率，默认值归一化频率程序：figure,spectrogram(x0,512,500,512,Fs,'yaxis');colormap('g ray');colorbar;title('语谱图');这里由于语谱图的数据量太大，很容易造成死机，因此这里采用截屏得到。