%http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html
% get a section of vowel
% source : http://www.angelfire.com/in2/yala/t4scales.htm
% Calculation for Equal-Tempered tuning [A4 = 440Hz]
% Hertz      Octave=1 Octave=2    Octave=3    Octave=4  Octave=5 Octave=6
% 
% 0   A     55.000     110.000     220.000     440.000   880.000  1,760.000
% 1 A#/Bb   58.270     116.541     233.082     466.164   932.328  1,864.655
% 2   B     61.735     123.471     246.942     493.883   987.767  1,975.533
% 3   C     65.406     130.813     261.626     523.251 1,046.502  2,093.005
% 4 C#/Db   69.296     138.591     277.183     554.365 1,108.731  2,217.461
% 5   D     73.416     146.832     293.665     587.330 1,174.659  2,349.318
% 6 D#/Eb   77.782     155.563     311.127     622.254 1,244.508  2,489.016
% 7   E     82.407     164.814     329.628     659.255 1,318.510  2,637.020
% 8   F     87.307     174.614     349.228     698.456 1,396.913  2,793.826
% 9 F#/Gb   92.499     184.997     369.994     739.989 1,479.978  2,959.955
% 10  G     97.999     195.998     391.995     783.991 1,567.982  3,135.963
% 11G#/Ab  103.826     207.652     415.305     830.609 1,661.219  3,322.438
% 12  A    110.000     220.000     440.000     880.000 1,760.000  3,520.000
clear
%[x,fs]=wavread('six.wav',[24120 25930]);
%[xx,fs]=wavread('sor1.wav');
%[xx,fs]=wavread('trumpet.wav');%measured by this program  E, 331 Hz
%[xx,fs]=wavread('violin3x.wav');%measured A2
%[xx,fs]=wavread('tz2.wav'); %measured B, around 493Hz)
%[xx,fs]=wavread('A4_violin.wav');%measured A4:436Hz
[xx,fs]=wavread('A4_oboe.wav');%measured A4,441Hz
%[xx,fs]=wavread('A5_flute.wav');%measured A5, 846Hz


%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=900
low_freq=200;

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




