matlab实现语音信号的频域分析及应用_matlab 生成语音信号

1.语音信号本质上是非平稳信号。但我们可以假设语音信号在一个短时间内是平稳的，这样我们用稳态分析方法处理非平稳信号。应用在傅立叶分析就是短时傅立叶变换。

语音的频域分析：包括语音信号的频谱、功率谱、倒频谱、频谱包络等.

常用频域分析方法：带通滤波器组法、Fourier 变换法、同态分析、线性预测法等。

2.倒谱分析：语音信号同态处理方法是一种设法将非线性问题转化为线性问题来进行处理的方法。它能将两个信号通过乘法合成的信号，或通过卷积合成的信号分开。

这种由卷积结果求得参与卷积的各信号分量—解卷。

对语音信号进行同态分析后，将得到语音信号的倒谱参数，因此同态分析也称为倒谱分析。

3.同态信号处理的基本原理：

典型卷积同态系统由三部分组成：

特征系统D*[]、线性系统L及逆特征系统[]。

4.语音倒谱的应用

基音周期估计：浊音信号的倒谱中存在峰值，它的出现位置等于该语音段的基音周期，而清音的倒谱中不存在峰值。利用倒谱的这个特点，我们可以进行语音的清浊音判决，并且可以估计浊音的基音周期。首先计算语音的倒谱，然后在可能出现的基音周期附近寻找峰值。如果倒谱峰值超过了预先设置的门限，则输入语音判断为浊音，其峰值位置就是基音周期的估计值；反之，如果没有超出门限的峰值的话，则输入语音为清音。

5.共振峰估计

对倒谱进行滤波，取出低时间部分进行逆特征系统处理，可以得到一个平滑的对数谱函数，这个对数谱函数显示了输入语音段的共振峰结构，同时谱的峰值对应于共振峰频率。通过此对数谱进行峰值检测，就可以估计出前几个共振峰的频率和强度。对于

6.语谱图

语谱图反映了语音信号的动态频率特性，在语音分析中具有重要的实用价值。

语谱图的时间分辨率和频率分辨率是由窗函数的特性决定的。时间分辨率高，可以看出时间波形的每个周期及共振峰随时间的变化，但频率分辨率低，不足以分辨由于激励所形成的细微结构，称为宽带语谱图，而窄带语谱图正好与之相反。

宽带语谱图可以获得较高的时间分辨率，反映频谱的快速时变过程；窄带语谱图可以获得较高的频率分辨率，反映频谱的精细结构。两者相结合，可以提供两种语音特性相关的信息。语谱图上因其不同的灰度，形成不同的纹路，称之为“声纹”。

7.绘制函数如下

语谱图绘制函数specgram

调用格式: specgram(data,nfft,Fs,window,numoverlap)

Data是语音信号，nfft是fft的长度，一般取1024或者512，fs就是采样率。
window是指窗的长度，一般和nfft相同即可，numoverlap是帧重叠的长度，取1/4 * nff 就可以了。

复倒谱：cceps，实倒谱：rceps

调用格式：y= cceps(x)；y= rceps(x)

8.下面是代码实现：

短时谱：（将语音信号保存为txt文件，即保存矩阵即可，----save函数）

clc;clear;
fid=fopen('mathvoice.txt', 'rt');%以文本的形式打开文件
x=fscanf(fid, '%f');
fclose(fid);
 
s1=x(12000: 40:153320);%取数组 x 的前 320 个数字
N=320;
s2=s1/max(s1);
figure(1);subplot(4,1,1);plot(s2);
xlabel('样点数');ylabel('幅值');
axis([0, 320, -1, 1]);
title('浊音原信号');
x2=enframe(s2,100,128);%分帧
ee=(x2(1,:));
 
%加Hamming窗
f=ee'.*hamming(length(ee));          %对选取的语音信号加Hamming窗
f1=f/max(f);                        %对加窗后的语音信号的幅值归一化
subplot(412)                       %画第三个子图
plot(f1)                           %画波形
%axis([0,256,-1.5,1.5])                %限定横纵坐标范围
xlabel('样点数')                    %横坐标名称
ylabel('幅度')                      %纵坐标名称                 
title ('窗选语音')                  %文字标注
 
% 矩形窗傅立叶变换
r=fft(s2,1024);                 %对信号ee进行1024点傅立叶变换
r1=abs(r);                    %对r取绝对值 r1表示频谱的幅度值
r1=r1/max(r1);                %幅值归一化
yuanlai=20*log10(r1);          %对归一化幅值取对数
signal(1:256)=yuanlai(1:256);    %取256个点，目的是画图的时候，维数一致
pinlv=(0:1:255)*8000/512;       %点和频率的对应关系
subplot(413)                  %画第五个子图
plot(pinlv,signal);              %画幅值特性图
xlabel('f/Hz')                  %横坐标名称
ylabel('对数幅度/dB')          %纵坐标名称
title ('加矩形窗时语音谱')     %文字标注
axis([0,4000,-80,15])           %限定横纵坐标范围
 
 
%加Hamming窗傅立叶变换
r=fft(f,1024);                      %对信号ee进行1024点傅立叶变换
r1=abs(r);                        %对r取绝对值 r1表示频谱的幅度值
r1=r1/max(r1);                    %幅值归一化
yuanlai=20*log10(r1);              %对归一化幅值取对数
signal(1:256)=yuanlai(1:256);        %取256个点，目的是画图的时候，维数一致
pinlv=(0:1:255)*8000/512;          %点和频率的对应关系
subplot(414)                      %画第七个子图
plot(pinlv,signal);                  %画幅值特性图
xlabel('f/Hz')                      %横坐标名称
ylabel('对数幅度/dB')              %纵坐标名称

倒谱图：

[s,fs,nbit]=wavread('123.wav');           %读入一段语音
b=s';                                   %将s转置
x=b(5000:5399);                         %取400点语音
N=length(x);                            %读入语音的长度
S=fft(x);                               %对x进行傅立叶变换
Sa=log(abs(S));                         %log为以e为底的对数
sa=ifft(Sa);                             %对Sa进行傅立叶逆变换
ylen=length(sa); 
for i=1:ylen/2
    sa1(i)=sa(ylen/2+1-i);
end
for i=(ylen/2+1):ylen
    sa1(i)=sa(i+1-ylen/2);
end
%绘图
figure(1);
subplot(2,1,1);
plot(x);
axis([0,400,-0.5,0.5])
title('截取的语音段');
xlabel('样点数');
ylabel('幅度');
subplot(2,1,2);
time2=[-199:1:-1,0:1:200];
plot(time2,sa1);
axis([-200,200,-0.5,0.5])
title('截取语音的倒谱');
xlabel('样点数');
ylabel('幅度');

共振峰检测：

waveFile='qinghua.wav '; 
[y, fs, nbits] = wavread(waveFile); 
time=(1:length(y))/fs; 
frameSize=floor(40*fs/1000);              % 帧长
startIndex=round(15000);                 % 起始序号
% startIndex=round(20000);                 % 起始序号
% endIndex=startIndex+frameSize-1;          % 结束序号 
endIndex=startIndex+frameSize-1;          % 结束序号 
frame = y(startIndex:endIndex);            % 取出该帧 
frameSize=length(frame);
frame2=frame.*hamming(length(frame));    % 加 hamming window 
rwy= rceps(frame2);                     % 求倒谱 
%ylen=length(rwy); 
ylen=max(size(rwy)) ;
cepstrum=rwy(1:floor(ylen/2)); 
 
%基音检测 
LF=floor(fs/500);
HF=floor(fs/70);
cn=cepstrum(LF:HF);
[mx_cep ind]=max(cn); 
 
%共振峰检测核心代码： 
% 找到最大的突起的位置 
NN=ind+LF; 
han= hamming (NN); 
cep=cepstrum(1:NN); 
ceps=cep.*han;                           % hamming window 
formant1=20*log(abs(fft(ceps))); 
formant(1:2)=formant1(1:2); 
for t=3:NN 
%--do some median filtering 
    z=formant1(t-2:t); 
    md=median(z); 
    formant2(t)=md; 
end 
for t=1:NN-1 
    if t<=2 
       formant(t)=formant1(t); 
    else
       formant(t)=formant2(t-1)*0.25+formant2(t)*0.5+formant2(t+1)*0.25;
    end 
end 
 
subplot(3,1,1); 
plot(cepstrum); 
title('倒谱'); 
xlabel('样点数');
ylabel('幅度')
axis([0,220,-0.5,0.5])
 
spectral=20*log10(abs(fft(frame2))); 
subplot(3,1,2); 
xj=(1:length(spectral)/2)*fs/length(spectral); 
 plot(xj,spectral(1:length(spectral)/2));  
title('频谱'); 
xlabel('频率/Hz');
ylabel('幅度/dB')
axis([0,5500,-100,50])
 
subplot(3,1,3); 
xi=(1:NN/2)*fs/NN; 
plot(xi,formant(1:floor(NN/2))); 
title('平滑对数幅度谱'); 
xlabel('频率/Hz');
ylabel('幅度/dB')
axis([0,5500,-80,0])

基音检测：

waveFile='beijing.wav '; 
[y, fs, nbits] = wavread(waveFile); 
time1=1:length(y); 
time=(1:length(y))/fs; 
frameSize=floor(50*fs/1000);               % 帧长
startIndex=round(5000);                   % 起始序号
endIndex=startIndex+frameSize-1;          % 结束序号 
frame = y(startIndex:endIndex);             % 取出该帧 
 
frameSize=length(frame);
frame2=frame.*hamming(length(frame));     % 加 hamming window 
rwy= rceps(frame2);                      % 求倒谱 
ylen=length(rwy); 
cepstrum=rwy(1:ylen/2); 
 
for i=1:ylen/2
    cepstrum1(i)=rwy(ylen/2+1-i);
end
for i=(ylen/2+1):ylen
    cepstrum1(i)=rwy(i+1-ylen/2);
end
 
%基音检测 
LF=floor(fs/500);                      %基音周期的范围是70Hz~500Hz
HF=floor(fs/70);
cn=cepstrum(LF:HF);
[mx_cep ind]=max(cn);
if mx_cep>0.08&ind>LF 
a= fs/(LF+ind);
else 
a=0; 
end 
pitch=a 
 
% 画图
figure(1); 
subplot(3,1,1); 
plot(time1, y); 
title('语音波形'); 
axis tight 
ylim=get(gca, 'ylim'); 
line([time1(startIndex),time1(startIndex)],ylim,'color','r');
line([time1(endIndex), time1(endIndex)],ylim,'color','r');
xlabel('样点数');
ylabel('幅度');
 
subplot(3,1,2); 
plot(frame); 
axis([0,400,-0.5,0.5])
title('一帧语音'); 
xlabel('样点数');
ylabel('幅度')
 
subplot(3,1,3); 
time2=[-199:1:-1,0:1:200];
plot(time2,cepstrum1); 
axis([-200,200,-0.5,0.5])
title('一帧语音的倒谱'); 
xlabel('样点数');
ylabel('幅度');

清浊音频谱图：

% 浊音的波形和短时频谱图(窗长256)
y=wavread('beijing.wav');
e=fra(256,128,y);              %对y分帧，帧长256，帧移128
ee=e(45,:);                   %选取第10帧
subplot(421)                  %画第一个子图
ee1=ee/max(ee);               %幅值归一化
plot(ee1)                     %画波形
xlabel('样点数')               %横坐标名称
ylabel('幅度')                 %纵坐标名称
title ('浊音')             %文字标注
axis([0,256,-1.5,1.5])           %限定横纵坐标范围
 
% 矩形窗傅立叶变换
r=fft(ee,1024);                 %对信号ee进行1024点傅立叶变换
r1=abs(r);                    %对r取绝对值 r1表示频谱的幅度值
r1=r1/max(r1);                %幅值归一化
yuanlai=20*log10(r1);          %对归一化幅值取对数
signal(1:256)=yuanlai(1:256);    %取256个点，目的是画图的时候，维数一致
pinlv=(0:1:255)*8000/512;       %点和频率的对应关系
subplot(425)                  %画第五个子图
plot(pinlv,signal);              %画幅值特性图
xlabel('f/Hz')                  %横坐标名称
ylabel('对数幅度/dB')          %纵坐标名称
title ('加矩形窗时语音谱')     %文字标注
axis([0,4000,-80,15])           %限定横纵坐标范围
 
%加Hamming窗
f=ee'.*hamming(length(ee));          %对选取的语音信号加Hamming窗
f1=f/max(f);                        %对加窗后的语音信号的幅值归一化
subplot(423)                       %画第三个子图
plot(f1)                           %画波形
axis([0,256,-1.5,1.5])                %限定横纵坐标范围
xlabel('样点数')                    %横坐标名称
ylabel('幅度')                      %纵坐标名称                 
title ('窗选语音')                  %文字标注
 
%加Hamming窗傅立叶变换
r=fft(f,1024);                      %对信号ee进行1024点傅立叶变换
r1=abs(r);                        %对r取绝对值 r1表示频谱的幅度值
r1=r1/max(r1);                    %幅值归一化
yuanlai=20*log10(r1);              %对归一化幅值取对数
signal(1:256)=yuanlai(1:256);        %取256个点，目的是画图的时候，维数一致
pinlv=(0:1:255)*8000/512;          %点和频率的对应关系
subplot(427)                      %画第七个子图
plot(pinlv,signal);                  %画幅值特性图
xlabel('f/Hz')                      %横坐标名称
ylabel('对数幅度/dB')              %纵坐标名称
title ('加Hamming窗时语音谱')    %文字标注
axis([0,4000,-80,15])               %限定横纵坐标范围
 
%清音的波形和短时频谱图(窗长256)
y=wavread('beijing.wav');
e=fra(256,128,y);                   %对y分帧，帧长256，帧移128
ee=e(5,:);                         %选取第2帧
subplot(422)                       %画第二个子图
ee1=ee/max(ee);                    %幅值归一化
plot(ee1)                          %画波形
xlabel('样点数')                    %横坐标名称
ylabel('幅度')                      %纵坐标名称
title ('清音')                      %文字标注
axis([0,256,-1.5,1.5])                %限定横纵坐标范围
 
% 矩形窗傅立叶变换
r=fft(ee,1024);                      %对信号ee进行1024点傅立叶变换
r1=abs(r);                          %对r取绝对值 r1表示频谱的幅度值
r1=r1/max(r1);                      %幅值归一化
yuanlai=20*log10(r1);                 %对归一化幅值取对数
signal(1:256)=yuanlai(1:256);           %取256个点，目的是画图的时候，维数一致
pinlv=(0:1:255)*8000/512;             %点和频率的对应关系
subplot(426)                        %画第六个子图
plot(pinlv,signal);                    %画幅值特性图
xlabel('f/Hz')                        %横坐标名称
ylabel('对数幅度/dB')                 %纵坐标名称
title('加矩形窗时语音谱')           %文字标注
axis([0,4000,-80,1])                   %限定横纵坐标范围
 
%加Hamming窗
f=ee'.*hamming(length(ee));            %对选取的语音信号加Hamming窗
f1=f/max(f);                         %对加窗后的语音信号的幅值归一化
subplot(424)                    %画第四个子图
plot(f1)                        %画波形
axis([0,256,-1.5,1.5])             %限定横纵坐标范围
xlabel('样点数')                 %横坐标名称
ylabel('幅度')                   %纵坐标名称
title ('窗选语音')               %文字标注
 
%加Hamming傅立叶变换
r=fft(f,1024);                   %对信号ee进行1024点傅立叶变换
r1=abs(r);                      %对r取绝对值 r1表示频谱的幅度值
r1=r1/max(r1);                  %幅值归一化
yuanlai=20*log10(r1);            %对归一化幅值取对数
signal(1:256)=yuanlai(1:256);      %取256个点，目的是画图的时候，维数一致
pinlv=(0:1:255)*8000/512;        %点和频率的对应关系
subplot(428)                    %画第八个子图
plot(pinlv,signal);                %画幅值特性图
xlabel('f/Hz')                    %横坐标名称
ylabel('对数幅度/dB')             %纵坐标名称
title ('加Hamming窗时语音谱')  %文字标注
axis([0,4000,-80,1])               %限定横纵坐标范围fid=fopen('voice2.txt','rt');      

代码仅供参考，部分代码来源于网络，但找不到出处了。侵删