%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]);

%since T=1/F, the bigger the freq, the smaller thme time
%if yopu are testing 440 Hz (A4 tone), set to 400-500Hz
high_freq=500
low_freq=100;

ms1_high_freq=fs/high_freq; % maximum speech Fx at 1000Hz, 
ms20_low_freq=fs/low_freq;  % minimum speech Fx at 50Hz
%sound(x,fs)
%
figure(1)
clf
%-----------------------------1 subplot(4,1,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)));
%
%-----------------------------2 subplot(4,1,2);---------------------------

% 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));

%-----------------------------3 subplot(4,1,3);---------------------------
%to m,ake sure you haev high frequcy , you must set range betwwen 
%high/low freqcy close to teh target freq. Ie, 400-500Hz for 
%typical tunign freq 440 (A4 in teh pinao scale)
%If your plot is between 1ms (=1000Hz) and 20ms (=50Hz), it is not enough
% set to 400-550Hz is much better

% ms1_high_freq=round(ms1_high_freq); %added by khw
% ms20_low_freq=round(ms20_low_freq);%added by khw
q=(ms1_high_freq:ms20_low_freq)/fs;
%plot lower side[freq:1000(ms1_high_freq:1milli-sec)->500(ms20_low_freq)Hz]
subplot(4,1,3);
size(C(ms1_high_freq:ms20_low_freq))
plot(q,abs(C(ms1_high_freq:ms20_low_freq)));
legend('Cepstrum,limited range,:Amplitude vs.Quefrency (sec)   khw');
xlabel('Quefrency (s)');
ylabel('Amplitude ');

%-----------------------------4 subplot(4,1,1);---------------------------
subplot(4,1,4);
C_length=length(C)
plot([1:C_length]/fs,abs(C));
legend('Cepstrum, all range,:Amplitude vs.Quefrency (sec)   khw');
xlabel('Quefrency (s)');
ylabel('Amplitude ');
%----------------------------5 ---calculate pitch ------------------
%usage : [max_val, index]=max(item)
[c,fx]=max(abs(C(ms1_high_freq:ms20_low_freq))); %fx is the sample index (not time)  
fprintf('>>>>>>>>>>>pitch found Fx=%gHz\n',fs/(ms1_high_freq+fx-1)); 




