clc; clear; close all;

%strings and variables set by user
word = 'obey'; %string for the word being analyzed
phoneme = 'o'; %string for the phoneme being analyzed
start_time_secs=0; %start and end times of phonemes so x-axis starts accurately
portion_size_variable = .025; %size of the portion taken from the phoneme to use in FFT
start_of_portion = .01; % offset of where to take portion from phoneme

wordfilename = strcat(word,'.wav');
phonemefilename = strcat(word,'PHONEME','.wav');

[s,fs, NBITS] = wavread(wordfilename);
phonemewav= wavread(phonemefilename);

%first figure: word and its specgrams
% SPECTROGRAM PLOTS

t=0:1/fs:(length(s)-1)/fs; %create time vector for word

T=0.025;    %let's set window size to 25 ms;
% Set digital window size to the closest power of 2 (e.g, 64, 128, 256, 512, ...). It makes computation more efficient
N=2^round(log2(T*fs));   
win=hamming(N);     %define a Hamming window
figure(1); 
subplot(311); plot(t,s); axis tight; grid on; ylabel('x(t)'); xlabel('time, s'); % Plot the word time waveform
title(strcat('Word: "',word,'" --  Phoneme under investigation: ',' /',phoneme,'/'));   %label the figure
subplot(312);specgram(s,N,fs,win,round(0.97*N)); xlabel('time, s'); ylabel('f, Hz');   %plot the narrowband spectrogram
title('Narrowband Spectrogram');
subplot(313), specgram(s,N,fs,hamming(N/4),round(0.9*N/4));xlabel('time, s'); ylabel('f, Hz');%plot the wideband spectrogram
title('Wideband Spectrogram');

%save figure and then write it to file
A = figure(1);
saveas(A, word, 'bmp')

%second figure, phoneme and its frequency plots
% DFT spectrum and LP smooth spectrum

figure(2); 
ph_time_vector = [start_time_secs:1/fs:start_time_secs+(length(phonemewav)-1)/fs];
subplot(311);plot(ph_time_vector,phonemewav); axis tight; grid on; ylabel('x(t)'); xlabel('time, s');
title(strcat('Phoneme /', phoneme, '/ from word "', word, '"'));

%take a small part of the phoneme and plot it before analyzing.
%20 ms covers both 3 periods in most voiced sounds
portion_size = portion_size_variable*fs;
portion = phonemewav(start_of_portion*fs:start_of_portion*fs+portion_size);
portion_time_vector = ((start_time_secs+start_of_portion):1/fs:(start_time_secs)+start_of_portion+(length(portion)-1)/fs);
subplot(312);plot(portion_time_vector,portion); axis tight; grid on; ylabel('x(t)'); xlabel('time, s')
title('25ms of phoneme');

%plot magnitude spectrum of small portion
bins=1024;     % number of FFT bins
fss=fft(portion,bins); fss=abs(fss(1:bins/2));     % Compute the baseband (d.c. to fs/s) magnitude spectrum
freqvector=0:fs/bins:fs/2-fs/bins; freqvector=freqvector';
subplot(313), plot(freqvector, db((fss/bins)));    % Plot the magnitude spectrum
axis tight; grid on; xlabel('f, Hz'); ylabel('|S(f)|, dB');
%A spectral envelope plot can be very useful
% Linear Prediction is a good way to compute the spectral envelope of voiced speech
p= 20; hold on;    % Set the order of prediction, p
[a,g]=lpc(portion,p);    % compute the Linear Prediction model
lspec=freqz(g,a, freqvector, fs);   % Compute the spectral envelope 
subplot(313), plot(freqvector, db(abs(lspec)), 'r'); % Superimpose the spectrial envelope on the Fourier magnitude spectrum
axis tight;
title('Magnitude Spectrum (Blue) and Smoothed Spectral Envelope (Red)');

B = figure(2);
saveas(B, strcat(word,'PHONEME'), 'bmp')


