plp.m

这是一个用于语音信号处理的工具箱

%  subroutine which computes m-th order (max 15) PLP model (described by 
%
%  function [a,rc,gain] = plp(speech, nwind, m, sf)
%
%  m+1 autoregressive coefficients a and gain or by m reflection 
%  coefficients rc) from nwind-points (max 512) of speech signal with
%  sampling frequency sf (max 20 000 Hz)
% 
%  Input
%  speech : speech signal (of type float)
%  nwind  : window (512 points typically) extractracted from speech
%  m      : order of model. (m+1) autoregressive coefficients
%  sf     : sampling frequency

%  Output
%  a      : auto regressive coefficients
%  rc     : reflection coefficients
%  gain   : model gain parameters
% 

function [a,rc,gain,spectra,audspe,nfilt] = plp(speech, nwind, m, sf)

    % get fft size from the window length 
    nfft = floor((log(nwind) / 0.693148)) + 1;

    % get number of spectral points 
    npoint = 2^nfft/2 + 1; 	

    % compute spectral weights coefficients
    [nfilt cb ibegen]=plp_audw(npoint, sf);  

    % compute IDFT coefficient table
    wcos = plp_cosf(m, nfilt);	   

    %     window input speech  
    speech = hamming(nwind) .* speech;

    %  compute speech power spectrum 
     spectra=plp_fft(speech,nwind,nfft);

    % compute auditory spectrum
    audspe=zeros(1,nfilt);
    for jfilt = 2:nfilt-1,
	for kk = ibegen(jfilt,1):ibegen(jfilt,2),
		icb=ibegen(jfilt,3)-ibegen(jfilt,1)+kk;
		audspe(jfilt)=audspe(jfilt)+spectra(kk)*cb(icb);
	end
    end

    % cubic root intensity-loudness compression */
    for ii = 2:nfilt-1,
	audspe(ii)=audspe(ii)^0.33;
    end

    % the first and last point of auditory spectrum */
    audspe(1) = audspe(2);
    audspe(nfilt) = audspe(nfilt-1);

    % inverse discrete fourier transform */
    nspt=2*(nfilt-1);
    for kk = 1:m+1,
	r(kk)=audspe(1);
	for ll=2:nfilt,
		r(kk)=r(kk)+audspe(ll)*wcos(ll,kk);
	end
	r(kk)=r(kk)/nspt;
    end


    % solution for autoregressive model */
    a(1)   = 1.0;
    alp(1) = r(1);
    rc(1)  = -r(2) / r(1);
    a(2)   = rc(1);
    alp(2) = r(1) + r(2) * rc(1);

    for mct = 2:m,
	s = 0.0;
	mct2 = mct + 2;
	alpmin = alp(mct);

	for ip = 1:mct,
		idx=mct2-ip;
		s=s+r(idx)*a(ip);
	end 

	rcmct = -s / alpmin;
	mh = mct/2 + 1;

	for ip = 2:mh,
	    ib = mct2 - ip;
	    aip = a(ip);
	    aib = a(ib);
	    a(ip) = aip + rcmct * aib;
	    a(ib) = aib + rcmct * aip;
	end	
	a(mct+1) = rcmct;
	alp(mct+1) = alpmin-alpmin*rcmct*rcmct;
	rc(mct) = rcmct;
    end 
    gain = alp(m+1);
    a=a(2:length(a));