LPC Envelope Example: Speech vowel ``ah''

% Let's make an "ah" [a] vowel:
% Ref: Dennis H. Klatt, "Software for a 
% cascade/parallel formant synthesizer,"
%      JASA, vol. 67, pp. 13-33, 1980.
F = [700, 1220, 2600]; % Formant frequencies in Hz
B = [130,   70,  160]; % Formant bandwidths in Hz

fs = 8192;  % Sampling rate in Hz 
	    % ("telephone quality" for speed)
R = exp(-pi*B/fs);     % Pole radii
theta = 2*pi*F/fs;     % Pole angles
poles = R .* exp(j*theta) % Poles
[B,A] = zp2tf(0,[poles,conj(poles)],1); % control/

f0 = 200; % Pitch in Hz
w0T = 2*pi*f0/fs;

nharm = floor((fs/2)/f0); % number of harmonics
sig = zeros(1,nsamps);
n = 0:(nsamps-1);
% Synthesize bandlimited impulse train
for i=1:nharm,
    sig = sig + cos(i*w0T*n);
end;
sig = sig/max(sig);
soundsc(sig,fs); % Let's hear it

% Now compute the speech vowel
speech = filter(1,A,sig);
soundsc([sig,speech],fs);
winspeech = w .* speech(1:length(w));

% and its spectrum
sspec = fft([winspeech,zeros(1,3*nplot)]);
dbsspecfull = 20*log10(abs(sspec));
dbsspec = dbsspecfull(1:nspec);
dbenv = 20*log10(abs(freqz(1,A,nspec)'));
dbsspecn = dbsspec + ones(1,nspec)*(max(dbenv) ...
                   - max(dbsspec)); % normalize
plot(f,[max(dbsspecn,-100);dbenv]); grid;

% Now measure spectral envelope various ways 
% First, the cepstral smoothing method:

rcep = ifft(dbsspecfull);  % real cepstrum
imagerr = norm(imag(rcep)) % check
   7.9798e-12
rcep = real(rcep);

period = round(fs/f0)
    41

aliasing = norm(rcep(nspec-10:nspec+10))/norm(rcep)
    0.0229

nw = 2*period-4; % almost 1 period left and right
if floor(nw/2) == nw/2, nw=nw-1; end; % make it odd
w = boxcar(nw)'; % rectangular window
wzp = [w(((nw+1)/2):nw),zeros(1,nfft-nw), ...
       w(1:(nw-1)/2)];  % zero-phase version
wrcep = wzp .* rcep;

% Display cepstral envelope
rcepenv = fft(wrcep);
imagerr = norm(imag(rcepenv)) % should be zero
rcepenvp = real(rcepenv(1:nspec));
rcepenvp = rcepenvp + ones(1,nspec)*(mean(dbenv)...
                    - mean(rcepenvp)); % normalize
plot(f,[max(dbsspecn,-100); dbenv; rcepenvp]);

% Repeat using Hanning window:

% Finally, let's do an LPC window.
% It had better be good because the LPC model 
% is exact for this example.

M = 6; % three formants

% compute Mth-order autocorrelation function:
rx = zeros(1,M+1)';
for i=1:M+1,
  rx(i) = rx(i) + speech(1:nsamps-i+1) ...
                * speech(1+i-1:nsamps)';
end

% prepare the M by M Toeplitz covariance matrix:
covmatrix = zeros(M,M);
for i=1:M,
  covmatrix(i,i:M) = rx(1:M-i+1)';
  covmatrix(i:M,i) = rx(1:M-i+1);
end

% solve "normal equations" for prediction coeffs

Acoeffs = - covmatrix \ rx(2:M+1)

Alp = [1,Acoeffs']; % LP polynomial A(z)

dbenvlp = 20*log10(abs(freqz(1,Alp,nspec)'));
dbsspecn = dbsspec + ones(1,nspec)*(max(dbenvlp) ...
                   - max(dbsspec)); % normalize
plot(f,[max(dbsspecn,-100);dbenv;dbenvlp]); grid;

---

Order a Hardcopy of Spectral Audio Signal Processing

Previous: Spectral Envelope
Next: LPC Functions in Matlab and Octave

---

written by Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.