% 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
Impulse Train
Speech Vowel and its Spectrum
speech = filter(1,A,sig); % impulse-train -> 'Ah' filter soundsc([sig,speech],fs); winspeech = w .* speech(1:length(w)); sspec = fft([winspeech,zeros(1,3*nplot)]); % interpolated spectrum dbsspecfull = 20*log(abs(sspec)); dbsspec = dbsspecfull(1:nspec); dbenv = 20*log(abs(freqz(1,A,nspec)')); dbsspecn = dbsspec + ones(1,nspec)*(max(dbenv) ... - max(dbsspec)); % normalize plot(f,[max(dbsspecn,-100);dbenv]); grid;
Spectral Envelope via Windowed Cepstrum
rcep = ifft(dbsspecfull); % real cepstrum imagerr = norm(imag(rcep))/norm(rcep) % check rcep = real(rcep); % imag part is just roundoff error 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-centered version wrcep = wzp .* rcep;
Rectangular-Windowed Cepstrum
% 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]);
Hanning-Windowed Cepstrum
Spectral Envelope via Linear Prediction
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 Toeplitx 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*log(abs(freqz(1,Alp,nspec)')); dbsspecn = dbsspec + ones(1,nspec)*(max(dbenvlp) ... - max(dbsspec)); % normalize plot(f,[max(dbsspecn,-100);dbenv;dbenvlp]); grid;