%http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html
% get a section of vowel


%[x,fs]=wavread('six.wav',[24120 25930]);
%[xx,fs]=wavread('sor1.wav');
%[xx,fs]=wavread('trumpet.wav');
%[xx,fs]=wavread('violin3x.wav');
%[xx,fs]=wavread('tz2.wav');
[xx,fs]=wavread('A4_violin.wav');
%[xx,fs]=wavread('A4_oboe.wav');
%[xx,fs]=wavread('A5_flute.wav');


%range=1810
duration_ms=30 %pick durmation_ms=30 ms
range=(duration_ms/1000)/(1/fs)  
length_xx=length(xx)
start=round(length_xx/2)-round(range/2) %pick middle
x=xx([start:start+range]);

ms1=fs/1000;                 % maximum speech Fx at 1000Hz
ms20=fs/50;                  % minimum speech Fx at 50Hz
%sound(x,fs)
%
%-----------------------------1 ------------------------------------
% plot waveform
t=(0:length(x)-1)/fs;        % times of sampling instants
subplot(4,1,1);
plot(t,x);
legend('Waveform: Amplitude vs. Time (second)');
xlabel('Time (s)');
ylabel('Amplitude');
%
% do fourier transform of windowed signal
Y=fft(x.*hamming(length(x)));
%
% plot spectrum of bottom 5000Hz, pitch should be below this freq.
hz5000=5000*length(Y)/fs;
f=(0:hz5000)*fs/length(Y);
subplot(4,1,2);
plot(f,20*log10(abs(Y(1:length(f)))+eps));
legend('Spectrum: Magnitude vs. Frequency');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
%
% cepstrum is DFT of log spectrum
C=fft(log(abs(Y)+eps));

%plot between 1ms (=1000Hz) and 20ms (=50Hz)
ms1=round(ms1); %added by khw
ms20=round(ms20);%added by khw
q=(ms1:ms20)/fs;
subplot(4,1,3);
size(C(ms1:ms20))
plot(q,abs(C(ms1:ms20)));
legend('Cepstrum,limited range,:Amplitude vs.Quefrency (s)   khw');
xlabel('Quefrency (s)');
ylabel('Amplitude ');

%-----------------------------4 ------------------------------------
subplot(4,1,4);
C_length=length(C)
plot([1:C_length]/fs,abs(C));
legend('Cepstrum, all range,:Amplitude vs.Quefrency (s)   khw');
xlabel('Quefrency (s)');
ylabel('Amplitude ');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[c,fx]=max(abs(C(ms1:ms20)));
CC=C(ms1:ms20);
size(C)
abs(CC)
size_C=size(C)
c
fx
fs
ms1
ms20
fprintf('Fx=%gHz\n',fs/(ms1+fx-1));


