subspe.m

matlab sound processing

function Subspe=Subspe()
if nargin < 1 
[ file, path ] = uigetfile( '*.wav', '选择待处理的带噪语音' );
[ signal, fs, nbits ] = wavread(strcat(path,file));
end
W = 256;              % 帧长 fs =8kHz
step = 1/4*256;          % 帧移 
fft_size =256 ;              % 计算FFT的窗长
signal= signal-mean( signal);
wnd=hamming(W);%长为W的汉明窗
y1=segment(signal,W,step,wnd);    %分帧:返回分帧后的信号矩阵  窗长（hamming）*帧数   和下面的fft相对应
Y1=fft(y1,fft_size);                        %各帧幅度谱512(fft)*?
Y1=abs(Y1(1:128,:));
Y1=Y1.^2;
yy1=[sum(abs(Y1(2:32,:)));sum(abs(Y1(33:64,:)));sum(abs(Y1(65:end,:)))];
yy1=10*log10(yy1);


%--------------------------------------------------------------------------


NIS = 20;                     % 认为初始无语音的帧数  (256 + 64*(20-1) )/8000秒
temp1=[];
temp2=[];
temp3=[];
temp4=[];
state_all=[];
newdist=[];
prelatch=[];
predist=[];
mux1=[];
mux2=[];
D_all=[];
Dist_all=[];
flag=0;
NoiseMargin_all=[];
alpha1=0.9;
state=0;
signal= signal-mean( signal);
wnd=hamming(W);%长为W的汉明窗
y=segment(signal,W,step,wnd);    %分帧:返回分帧后的信号矩阵  窗长（hamming）*帧数   和下面的fft相对应
Y=fft(y,fft_size);                        %各帧幅度谱512(fft)*?
numberOfFrame=length(Y(1,:));
for tt=1:numberOfFrame-1
 sp=Y(:,tt);
 spangle=angle(sp) ;                         % 取出相位信息 get angle
 sp_am = abs(sp(1:fft_size/2+1));            % 取出幅度信息 get Specrogram
 sp=sp_am.^2 ;                           % 功率谱 we have the power spectrum here 
 if tt<NIS
    if  tt==1
            noise =sp;
    end
    noise = 0.55*noise +0.45*sp;    
    state=0;
    prestate=0;
    state_all=[state_all state];
    D_all=[D_all 0];
    if tt==NIS-1
        PreSpectralDist= 10*(log10(sp)-log10(noise));%SNR
        II=find(PreSpectralDist<0);
        PreSpectralDist(II)=0.01;
    end
    temp1=[temp1 [0.1 0.1 0.1]'];
 else
    SpectralDist= 10*(log10(sp)-log10(noise));%SNR
    II=find(SpectralDist<0);
    SpectralDist(II)=0.005;
    SpectralDist=alpha1*SpectralDist+(1-alpha1)*PreSpectralDist;
    PreSpectralDist=SpectralDist;
    temp_dist=SpectralDist(2:end);%舍去直流       
    temp_dist=reshape(temp_dist,32,4);%分成16个平均的子带
    band1(:,1)=temp_dist(:,1);
    band2(:,1)=temp_dist(:,2);
    band3(:,1)=[temp_dist(:,3); temp_dist(:,4)];    
    newdist=[];
    newmean=[];
      for i=1:3
      if i<2
       temp_band=band1(:,i);
       temp_band=sort(temp_band);    
       newdist=[newdist mean(temp_band(end-1:end))];
      elseif i<3
       temp_band=band2(:,i-1);
       temp_band=sort(temp_band);    
       newdist=[newdist mean(temp_band(end-1:end))];
    else
       temp_band=band3(:,i-2);
       temp_band=sort(temp_band);    
       newdist=[newdist mean(temp_band(end-2:end))];
    end    
  end
   newdist=newdist';
   if tt>NIS
   newdist=0.9*newdist+0.1*temp1(end);
   end
   temp1=[temp1 newdist];%做保存
end
end






figure;
% subplot(311),plot((yy1(1,:)-min(yy1(1,:)))/(max(yy1(1,:))-min(yy1(1,:))));hold ;
% plot((temp1(1,:)-min(temp1(1,:)))/(max(temp1(1,:))-min(temp1(1,:))),'r');
% subplot(312),plot((yy1(2,:)-min(yy1(2,:)))/(max(yy1(2,:))-min(yy1(2,:))));hold;
% plot((temp1(2,:)-min(temp1(2,:)))/(max(temp1(2,:))-min(temp1(2,:))),'r');
% subplot(313),plot((yy1(3,:)-min(yy1(3,:)))/(max(yy1(3,:))-min(yy1(3,:))));hold;
% plot((temp1(3,:)-min(temp1(3,:)))/(max(temp1(3,:))-min(temp1(3,:))),'r');

subplot(311),p1=plot((yy1(1,:)-min(yy1(1,:)))/(max(yy1(1,:))-min(yy1(1,:))),'b');grid on;
xlabel('Frame');
ylabel('归一化幅度值');
title('带噪语音信号第1个子带归一化后能量和(0--1000Hz)');
subplot(312),p2=plot((temp1(1,:)-min(temp1(1,:)))/(max(temp1(1,:))-min(temp1(1,:))),'r');grid on;
xlabel('Frame');
ylabel('归一化幅度值');
title('带噪语音信号第1个子带归一化后信噪比最大值(0--1000Hz)');
set(p1,'LineWidth',2);
set(p2,'LineWidth',2);

% subplot(311),plot(yy1(1,:));hold ;
% plot(temp1(1,:),'r');
% subplot(312),plot(yy1(2,:));hold ;
% plot(temp1(2,:),'r');
% subplot(313),plot(yy1(3,:));hold ;
% plot(temp1(3,:),'r');
subplot(313),p1=plot(yy1(1,:),'b');hold on;
p2=plot(temp1(1,:),'r--');grid on;
legend('归一化前能量和','归一化前信噪比最大值');
xlabel('Frame');
ylabel('幅度值');
title('带噪语音信号第1个子带未归一化两种方法的特征值(0--1000Hz)');
set(p1,'LineWidth',2);
set(p2,'LineWidth',2);




%====================分帧子函数============================================
function Seg=segment(signal,W,SP,Window)
   if nargin<4
      Window=hamming(W);
   end
  Window=Window(:); %make it a column vector列向量

  L=length(signal);
  SP;%帧移
  N=fix((L-W)/SP +1); %number of segments帧数

  Index=(repmat(1:W,N,1)+repmat((0:(N-1))'*SP,1,W))';%构造W*N的矩阵
  
  hw=repmat(Window,1,N);
  %Window是一个W*N矩阵。
  Seg=signal(Index).*hw;
  %取signal的与Index的元素分别对应的排列成W*N矩阵再与hw点*
%------------------------end 分帧----------------------------------