# H_DplotVowels.R

# (c) 2011 by Qingyang Yan and her classmates in Ling 600.01.

# Script for plotting formants for vowel stimuli chosen by someone
# as closest to the vowels in English /h__d/ words when acting as 
# the "listener" in the Method of Adjustment experiment, and then 
# overlaying the vowel formants measured from tokens of the English
# /h__d/ words produced by the same someone when acting as the 
# "talker" in a production experiment. The plotting is done in an 
# ERB space, with choices of plot type as "s" (square) and of 
# xlim and ylim that will make the x- and y-axes cover the same 
# range of values. 

# Start by setting the working directory to the directory where 
# you have the txt file containing the vowel formant values that
# you want to plot (and where you want to store the output pictures.
# Change the path name as below.  
setwd('c:/Lx600.01')

# Read in the data as a data frame.  This assumes that you have 
# a file with three tab-separated columns, that looks something
# like this, which is the contents of the sampleVowels.txt file
# that is linked into the course web page.  Change the file name 
# as appropriate.
#
# vowel	F1_Hz	F2_Hz
# [i]	360	2400
# [e]	680	1900
# [a]	800	1300
# [o]	590	900
# [u]	390	900
formants=data.frame(read.table("CollinsFormants.txt",
	header=TRUE, sep="\t"))

# If the vowel column includes IPA font that won't show up in
# the plot, you may want to change the first column to be the
# orthographic form of the word, like this for Michael Collins's
# formants file. 
formants$vowel=c("had","hod","hawed","head","heard","hayed",
	"hid","heed","hoed","hood","HUD","whoed")
# Alternatively, if you've used only ASCII symbols, and want a
# cleaner plot, you could change the corresponding column in the
# MOA results data frame that you read in next.

# Here is the code that Jefferson Barlew made to read in the 
# text file of results from his experiment.  Note that this 
# file was one that already included the stimuli's formant 
# values, added by running the code in the addFormantValues.R
# script in the http://ling.osu.edu/~mbeckman/ChooseVowels 
# directory. 
results = read.table("jsb-perc-res.txt", header=TRUE, sep="\t")
# When Qingyang wrote the script in class, she already had read
# in results from Michael Collins's experiment, using the code 
# in the addFormantValues.R, so she didn't need to do that in
# this script.  

# Make two new columns for the formants data frame that are the ERB 
# values for the measured formant values in the produced tokens, 
# and name them appropriately.  The formula for calculating ERB here 
# is taken from the following locus classicus:
#   Moore, B., & Glasberg, B.R. (1987). "Formulae describing frequency
#   selectivity as a function of frequency and level, and their use 
#   in calculating excitation patterns." Hearing Research, 28, 
#   209-225.
# ERB stands for "equivalent rectangular bandwidth" and is a slightly 
# modified timbre scale that is close to Barks, differing mostly in 
# low frequencies.  
formants$F1_ERB=
	11.17*log(((formants$F1_Hz/1000)+0.312)/((formants$F1_Hz/1000)+14.675))+43
formants$F2_ERB=
	11.17*log(((formants$F2_Hz/1000)+0.312)/((formants$F2_Hz/1000)+14.675))+43

# Set up a window that is 6" wide by 6" high, using 12point type.
windows(w=6,h=6,pointsize=12)   # Like this when using Windows
quartz(w=6,h=6,pointsize=12)    # Like this when using MacOS
# Specify plotting parameters to make the font be a serif font and 
# make the plot be square, and have no outside margin around the plot.
par(family="serif",pty="s",oma=rep(0,4))

# Plot the formants from the <h__d> words that Michael produced, 
# in ERB space, with the x-axis going from 
# 23 to 14 and the y-axis going from 15 to 6.
plot(formants$F2_ERB,formants$F1_ERB,xlim=c(23,14),ylim=c(15,6),
	xlab="second formant measured at vowel mid point (ERB)",
	ylab="first formant measured at vowel mid point (ERB)",pch=19)
# Add text to label the points.
text(formants$F2_ERB,formants$F1_ERB,formants$vowel,adj=c(1.5,0))

# Make the results from Michael choosing vowel stimuli in the MOA 
# experiment be in ERB, too.
results$F1_ERB=
	11.17*log(((results$F1/1000)+0.312)/((results$F1/1000)+14.675))+43
results$F2_ERB=
	11.17*log(((results$F2/1000)+0.312)/((results$F2/1000)+14.675))+43

# Overlay the points for the stimuli that Michael chose onto the 
# same plot as the points for the formants he measured in his 
# own productions of the same words.
points(results$F2_ERB, results$F1_ERB, pch=19, col="red")
text(results$F2_ERB, results$F1_ERB, results$stimulus,adj=c(1.5,0), col="red")

# Now save the plot to a postscript file to embed in your tex file,
# or save it to a pdf to embed in your MS Word doc file, or copy it 
# to your "clipboard" to embed in your doc file, or ... 
# That is, do whatever you normally do to get the plot from the R
# graphics device into your report for the exercise.