% Steven Usma
% EEN 540
% Project 3
% 4/24/04


function [] = cepsdemo()
clear; close all;

wav_file='usma.wav';
voiceds=.334;
voicedf=.357;
unvoiceds = .110;
unvoicedf = .130;

%read in wave file
[x, fs, nbits] = wavread(wav_file);
ts=1/fs;
t=(1:length(x))/fs;


%Plot the entire time waveform
figure(1);
plot(t,x);
axis tight;
grid on;
ylabel('Amp.');
title('spoken "Steven Usma"');
hold on;



%Narrowband spectrogram
T=0.020;
N=T*fs;
Lfft=512;
figure(2);
subplot(211);
overlap=round(.95*N);
specgram(x,Lfft,fs,hamming(N), overlap);
zz=specgram(x,Lfft,fs,hamming(N), overlap);
size(zz)
%spectrogram(x,hamming(N),overlap,Lfft,fs,'yaxis')
title(sprintf('Narrowband Spectrogram overlap=%i timepts=%i',overlap,size(zz,2)));


%Wideband spectrogram, 4 ms window
subplot(212);
overlap=round(.95*N/5);
%spectrogram(x,hamming(N/5),overlap,Lfft,fs,'yaxis')
specgram(x,Lfft,fs,hamming(N/5), overlap);
zz=specgram(x,Lfft,fs,hamming(N/5), overlap);
size(zz)
title(sprintf('Wideband Spectrogram overlap=%i timepts=%i',overlap,size(zz,2)));

dopause

%voiced
vstart=voiceds*fs;
vend=voicedf*fs;
voiced = x(vstart:vend)';
t_voiced=(vstart/fs:1/fs:vend/fs);

%unvoiced
ustart=unvoiceds*fs; uend=unvoicedf*fs;
unvoiced = x(ustart:uend)';
t_unvoiced=(ustart/fs:1/fs:uend/fs);

%Cepstrum and LPC for voiced
vcepstrum = fft(voiced'.*hamming(length(voiced)),512);
vcepstrum=log(abs(vcepstrum));
vcepstrum=real(ifft(vcepstrum));
Len=length(vcepstrum);

%Set Quefrency for Slow/Fast Liftering
Quefreq=42;

%Plot Voiced waveform and Real Cepstrum
figure(3);
subplot(211);
plot(t_voiced,voiced);
axis tight;
title('Voiced Waveform "i"');
ylabel('Amp.');
xlabel(['sec']);

subplot(212);
t_cep=((1:length(vcepstrum))-257)/fs * 1000;
plot(t_cep, real(fftshift(vcepstrum)));
axis tight;
title('Real Cepstrum of Voiced Frame');
xlabel('Quefreq. (ms)');
ylabel('Amp.');

%Slow liftering
vsystem = [vcepstrum(1);
2*vcepstrum(2:Quefreq);
vcepstrum(Quefreq+1);
zeros(Len/2-Quefreq-1,1)];
figure;
vsenvelope =real(fft(vsystem,Len));
vsenvelope=exp(vsenvelope);
Freq = ((1:length(vsenvelope))/length(vsenvelope)*fs)';
figure(4);
subplot(211);
plot(Freq(1:Len/2),(vsenvelope(1:Len/2)));
axis tight; title('Spectral Envelope of Voiced Frame ');
xlabel('Freq. (Hz)');
ylabel('Amp. ');

%Fast liftering
vsource = [zeros(Quefreq-1,1);
vcepstrum(Quefreq);
2*vcepstrum(Quefreq+1:(Len)/2-1);
vcepstrum((Len)/2)];
voiced_spec_details=real(fft(vsource,Len));
voiced_spec_details = exp(voiced_spec_details);
subplot(212);
plot(vsource);
plot(Freq(1:Len/2), (voiced_spec_details(1:Len/2)));
axis tight;
title('Spectral Details of Voiced Frame ');
xlabel('Freq. (Hz)');
ylabel('Amp.');

%Linear Prediction Coefficients
P=14;
[a,g]=lpc(voiced,P);
voiced_predicted=filter(a,1,voiced);
figure(5);
subplot(211);
plot(t_voiced,voiced);
title('Voiced waveform with Linear Prediction Error');
axis tight;
grid on;
xlabel('sec');
ylabel('Amp.')
hold on;
plot(t_voiced,voiced_predicted,'r');
axis tight;

%Magnitude spectrum of the Error
M=512;
error=abs(fft(voiced_predicted,M));
error=error(:,1:M/2);
L=((1:length(error)))/M*fs;

%Frequency Response of filter with coefficients from lpc
lspec=freqz(1,a,L,fs);

%Plot spectral envelope on same axis as error magnitude spectrum
subplot(212);
plot(L, (abs(lspec/max(lspec))));
title('Linear Prediction Spectrum with Error Signal Spectrum');
xlabel(['Freq. (Hz)']);
ylabel('Amp. ');
hold on;
plot(L,(error),'r');
axis tight;
grid on;


%Cepstrum and LPC for unvoiced
ucepstrum = fft(unvoiced'.*hamming(length(unvoiced)),512);
ucepstrum=log(abs(ucepstrum));
ucepstrum=real(ifft(ucepstrum));
Len=length(ucepstrum);


Quefreq=31;       %Quefrequency for Slow and Fast Liftering

%Plot Unvoiced waveform and Real Cepstrum
figure(6);
subplot(211);
plot(t_unvoiced,unvoiced);
axis tight;
title('Unvoiced Waveform "s" ');
ylabel('Amp');
xlabel(['sec']);
subplot(212);
t_cep_un=((1:length(ucepstrum))-257)/fs * 1000;
plot(t_cep_un, real(fftshift(ucepstrum)));
axis tight; title('Real cepstrum of Unvoiced Frame');
xlabel('Quefreq. (ms)');
ylabel('Amp.');

%Window area less than quefreq for slow liftering
usystem = [ucepstrum(1);
2*ucepstrum(2:Quefreq);
ucepstrum(Quefreq+1);
zeros(Len/2-Quefreq-1,1)];
unvsenvelope =real(fft(usystem,Len));
unvsenvelope=exp(unvsenvelope);
Freq = ((1:length(unvsenvelope))/length(unvsenvelope)*fs)';
figure(7);
subplot(211);
plot(Freq(1:Len/2),(unvsenvelope(1:Len/2)));
axis tight; title('Spectral Envelope of Unvoiced Frame');
xlabel('Freq. (Hz)');
ylabel('Amp. ');

%Window area greater than quefreq for fast liftering
unvsource = [zeros(Quefreq-1,1);
ucepstrum(Quefreq);
2*ucepstrum(Quefreq+1:(Len)/2-1);
ucepstrum((Len)/2)];
unvoiced_spec_details=real(fft(unvsource,Len));
unvoiced_spec_details = exp(unvoiced_spec_details);
subplot(212);
plot(Freq(1:Len/2), (unvoiced_spec_details(1:Len/2)));
axis tight;
title('Spectral Detail of Unvoiced Frame');
xlabel('Freq. (Hz)');
ylabel('Amp. ');

%Linear prediction coefficients
P=14;
[a,g]=lpc(unvoiced,P);
unvoiced_predicted=filter(a,1,unvoiced);
figure(8);
subplot(211);
plot(t_unvoiced,unvoiced);
title('Unvoiced waveform with Linear Prediction Error');
axis tight;
hold on;
xlabel('sec');
ylabel('Amp.');
plot(t_unvoiced,unvoiced_predicted,'r');
axis tight;
grid on;

%Magnitude spectrum of the Error
M=512;
error2=abs(fft(unvoiced_predicted,M));
error2=error2(:,1:M/2);
L=((1:length(error2)))/M*fs;

%Frequency Response of filter with coefficients from lpc
lspec=freqz(1,a,L,fs);

%Plot spectral envelope on same axis as error magnitude spectrum
subplot(212);
plot(L, (abs(lspec/max(lspec))));
title('Linear Prediction Spectrum with Error Signal Spectrum');
xlabel('Freq. (Hz)');
ylabel('Amp.');
hold on;
plot(L,(error2),'r');
axis tight;
grid on;

end

function [] = dopause()
disp('press any key to continue')
pause
end
