SII
Calculate ANSI S3.5-1997 Speech Intelligibility Index

Package index
Search the SII package
Vignettes
Functions
14
Source code
8
Man pages
4
Home
/
R-Forge
/
SII
/
inst/doc/SII.R

inst/doc/SII.R
In SII: Calculate ANSI S3.5-1997 Speech Intelligibility Index 

### R code from vignette source 'SII.Rnw'

###################################################
### code chunk number 1: load.package
###################################################
library(SII)


###################################################
### code chunk number 2: load.data
###################################################
## Table 1: Critical band SII procedure constants
data("critical")
head(critical)

## Table 2:Equally contributing (17 band) critical band SII 
## procedure constants
data("equal")
  
## Table 3: One-third octave band SII procedure constants
data("onethird")

## Table 4: Octave band SII procedure constants
data("octave")
      
## Overall SPL constants
data("overall.spl")
overall.spl


###################################################
### code chunk number 3: load.alternative.data
###################################################
## Table B.1: Critical band importance functions for various speech tests.
data(sic.critical)
head(sic.critical)

## Table B.2: One-third octave band importance functions for various speech tests.
data(sic.onethird)

## Table B.3: Octave band importance functions for various speech tests.
data(sic.octave)


###################################################
### code chunk number 4: sicplot
###################################################
data(sic.critical)
ngroup <- ncol(sic.critical)
matplot(x=sic.critical[,1], y=sic.critical[,-1],
        type="o",
        xlab="Frequency, Hz",
        ylab="Weight",
        log="x",
        lty=1:ngroup,
        col=rainbow(ngroup)
)
legend(
       "topright",
       legend=names(sic.critical)[-1],
       pch=as.character(1:ngroup),
       lty=1:ngroup,
       col=rainbow(ngroup)
       )



###################################################
### code chunk number 5: header
###################################################
args(sii)


###################################################
### code chunk number 6: ExampleC.1
###################################################
sii.C1 <- sii(
              speech   = c(50.0, 40.0, 40.0, 30.0, 20.0,  0.0),
              noise    = c(70.0, 65.0, 45.0, 25.0,  1.0,-15.0),
              threshold= c( 0.0,  0.0,  0.0,  0.0,  0.0,  0.0),
              method="octave"
              #, importance="SII"
              #, importance=octave$Ii
              , importance="CST"
              )
round(sii.C1$table[,-c(5:7,13)],2)
sii.C1$sii


###################################################
### code chunk number 7: ExampleC.1
###################################################
sii.C2 <- sii(
              speech   = rep(54.0, 18),
              noise    = c(40.0, 30.0, 20.0, rep(0, 18-3) ),
              threshold= rep(0.0,  18),
              method="one-third"
              )
sii.C2$table[1:3,1:8]
sii.C2$sii


###################################################
### code chunk number 8: ExampleC.1
###################################################
sii.C1 <- sii(
              speech   = c(50.0, 40.0, 40.0, 30.0, 20.0,  0.0),
              noise    = c(70.0, 65.0, 45.0, 25.0,  1.0,-15.0),
              threshold= c( 0.0,  0.0,  0.0,  0.0,  0.0,  0.0),
              method="octave"
              #, importance="SII"
              #, importance=octave$Ii
              , importance="CST"
              )
round(sii.C1$table[,-c(5:7,13)],2)
sii.C1$sii


###################################################
### code chunk number 9: LoadData (eval = FALSE)
###################################################
## library(gdata)
## patInfo <- read.xls("../AI subject list.xls")


###################################################
### code chunk number 10: ExamineData (eval = FALSE)
###################################################
## ## Information about variables
## str(patInfo)


###################################################
### code chunk number 11: ExamineData (eval = FALSE)
###################################################
## ## First 6 rows
## head(patInfo)


###################################################
### code chunk number 12: vars (eval = FALSE)
###################################################
## ## measured frequencies
## freq <- c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000)
## 
## ## columns containing frequencies for the right/left ear
## rt.cols <- paste("PTR",freq, sep="")
## lt.cols <- paste("PTL",freq, sep="")
## 
## rt.vals <- patInfo[,rt.cols]
## lt.vals <- patInfo[,lt.cols]


###################################################
### code chunk number 13: missing (eval = FALSE)
###################################################
## rt.vals[rt.vals==-888] <- NA
## lt.vals[rt.vals==-888] <- NA


###################################################
### code chunk number 14: fun (eval = FALSE)
###################################################
## fun <- function(X)
##   {
##     ret <- try( 
##                sii(X, 
##                    speech="raised",
##                    threshold=c(15,15,20,25,35,35,45,50), 
##                    freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
##                    importance="SII",
##                    interpolate=TRUE
##                    )$sii
##                )
##   if("try-error" %in% class(ret))
##     return(NA)
##   else
##     return(ret)
##   }
## 
## ## Test it
## fun( rt.vals[1,] )


###################################################
### code chunk number 15: CalculateSII (eval = FALSE)
###################################################
## sii.right <- apply(rt.vals, 1, fun )
## 
## sii.left  <- apply(lt.vals, 1, fun )


###################################################
### code chunk number 16: AddSIIToTable (eval = FALSE)
###################################################
## patInfo$"SII.right" <- sii.right
## patInfo$"SII.left" <- sii.left
## 
## tail(patInfo)


###################################################
### code chunk number 17: SavePatInfo (eval = FALSE)
###################################################
## write.table(patInfo, 
##             file="../AI subject list-SII.xls",
##             row.names=FALSE, 
##             sep="\t"
##             )


###################################################
### code chunk number 18: SIIFileFun (eval = FALSE)
###################################################
## sii.dina <- function(infile, outfile, verbose=TRUE) 
##   {
##     if(verbose)
##       cat("\nLoading data file '", infile, "'...\n", sep="")
##     ## Load the data
##     library(gdata)
##     patInfo <- read.xls(infile)
##     
##     ## measured frequencies
##     freq <- c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000)
## 
##     if(verbose)
##       cat("\nExtracting hearing thresholds...\n")
## 
##     ## columns containing frequencies for the right/left ear
##     rt.cols <- paste("PTR",freq, sep="")
##     lt.cols <- paste("PTL",freq, sep="")
##     
##     rt.vals <- patInfo[,rt.cols]
##     lt.vals <- patInfo[,lt.cols]
##     
##     ## Handle missing code '-888'
##     rt.vals[rt.vals==-888] <- NA
##     lt.vals[rt.vals==-888] <- NA
##     
##     ## define function to compute SII with our defaults
##     fun <- function(X)
##       {
##         ret <- try( 
##                    sii(X, 
##                        speech="raised",
##                        threshold=c(15,15,20,25,35,35,45,50), 
##                        freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
##                        importance="SII",
##                        interpolate=TRUE
##                        )$sii #$
##                    )
##       if("try-error" %in% class(ret))
##         return(NA)
##       else
##         return(ret)
##       }
##     
##     ## Calculate SII
##     if(verbose)
##       cat("\nCalculating right ear SII...\n")
##     sii.right <- apply(rt.vals, 1, fun )
##     if(verbose)
##       cat("\nCalculating left ear SII...\n")
##     sii.left  <- apply(lt.vals, 1, fun )
##     
##     ## Add back onto the table
##     patInfo$"SII.right" <- sii.right
##     patInfo$"SII.left" <- sii.left
## 
##     if(verbose)
##       cat("\nWriting new data table as '", outfile, "'...\n", sep="")
## 
##     ## Save file
##     write.table(patInfo, 
##                 file=outfile,
##                 row.names=FALSE, 
##                 sep="\t"
##                 )
## 
##     if(verbose)
##       cat("\nDone.\n\n")
##   }


###################################################
### code chunk number 19: RunFun (eval = FALSE)
###################################################
## sii.dina(infile="../AI subject list.xls",
##          outfile="../AI subject list-SII.xls") 


###################################################
### code chunk number 20: CompareAItoSII (eval = FALSE)
###################################################
## library(xtable)
## xt <- xtable(patInfo[,c("AI","SII.right","SII.left")], 
##              caption="Comparison of original AI and new right and left SII values \\label{table}", 
##              digits=2)
## print(xt)


###################################################
### code chunk number 21: PlotAIandSII (eval = FALSE)
###################################################
## ## put histograms on the diagonal
## panel.hist <- function(x, ...)
##   {
##     usr <- par("usr"); on.exit(par(usr))
##     par(usr = c(usr[1:2], 0, 1.5) )
##     h <- hist(x, plot = FALSE)
##     breaks <- h$breaks; nB <- length(breaks)
##     y <- h$counts; y <- y/max(y)
##     rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
##   }
## pairs.2 <- function(x)
##   pairs(x, 
##       panel=panel.smooth,
##       cex = 1.5, 
##       pch = 24, 
##       bg="light blue",
##       diag.panel=panel.hist, 
##       cex.labels = 2, 
##       font.labels=2
##       )
## 
## 
## pairs.2(patInfo[,c("AI","SII.right","SII.left")])


###################################################
### code chunk number 22: SII.Rnw:438-479
###################################################
latexhelp <- function(topics, package=NULL)
  {
    if(missing(topics) && missing(package))
      stop("Help on what?")
    else if (missing(topics))
      {
        callObj <- call("library", help=package)
        topicStrings <- eval(callObj)$info[[2]]
        topicStrings <- lapply( topicStrings, strsplit, split=" ")
        topicStrings <- sapply( topicStrings, function(x) x[[1]][1])
        topicStrings <- topicStrings[topicStrings > " "]
        topics <- topicStrings
        }
    
    for(topic in topics)
      {
        if(!is.character(topic))
          topic <- deparse(substitute(topic))
       
        callObj <- call("help", topic=topic, package=package)
        
        ## Get the file path
        help.file<- as.character(eval(callObj))

        ## Get R to translate it to latex and display it
        tools::Rd2latex(utils:::.getHelpFile(help.file))
        
        ## new page
        cat("\n")
        cat("\\clearpage")
        cat("\n")
      }
  }


tmp <- function (x, ...) 
{
}

latexhelp(package="SII")



###################################################
### code chunk number 23: splineFit
###################################################

THDI=c(25,25,30,35,45,45,55,60)
freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000)

sii.freqs <- SII:::sii.constants[,1]

xlist <- sort(c(SII:::sii.constants[, 1], freq))
ylist <- rep(NA, length=length(xlist))
names(ylist) <- xlist
ylist[ as.character(freq) ] <- THDI


###################################################
### code chunk number 24: interSpline
###################################################
library(splines)
ispl <- interpSpline( THDI ~ freq )
ispl <- predict(ispl, sii.freqs)$y
      
ispl.l <- interpSpline( THDI ~ log(freq) )
ispl.l <- predict(ispl.l, log(sii.freqs) )$y

approx.l <- function(x,y,xout,...)
  { 
    retval <- approx(log(x), y, log(xout), ...)
    retval$x <- xout
    retval
  }



###################################################
### code chunk number 25: Spline_comparison_function
###################################################
doplot <- function(logx=FALSE, separate=FALSE)
  {
    if(separate)
      layout( cbind( c(1,2,3,7), c(4,5,6,7) ))
    tmp <- function() {
      plot(x=freq, y=THDI, col="black", cex=2, log=if(logx) "x" else "", 
           xlab="Frequency (Herz)", ylab="Threshhold of Detection (dB)",
           xlim=range(sii.freqs), ylim=c(20,60) )
      lines(SII:::sii.constants[, 1], SII:::sii.constants[, "Ti'.THDN"], col="red", lwd=3)
      abline(v=sii.freqs,lty=3)
    }
    
    tmp()
    lines( x=SII:::sii.constants[, 1], y=ispl, col="blue", lwd=2)
    
    if(separate) tmp()
    lines( x=SII:::sii.constants[, 1], y=ispl.l, col="cyan", lwd=2)
    
    if(separate) tmp()
    lines( spline  (x=freq, y=THDI), col="green", lwd=2)
    
    if(separate) tmp()
    lines( spline  (x=freq, y=THDI, method="natural", xout=sii.freqs), col="orange", lwd=2)
    
    if(separate) tmp()
    lines( approx  (x=freq, y=THDI, method="linear",  xout=sii.freqs, rule=2), col="magenta", lwd=2)
    
    if(separate) tmp()
    lines( approx.l(x=freq, y=THDI, method="linear",  xout=sii.freqs, rule=2), col="brown",   lwd=3)

    if(separate) plot.new()
    legend(if(separate) "center" else { if(logx)  "topleft" else "bottomright"  },
           legend=c(
             "Measured data",
             "Matlab INTERP1(X, Y, XI, 'spline', 'extrap')",
             "R's predict( interSpline(X,Y), XI )",
             "R's predict( interSpline(log(X),Y), log(XI) )",
             "R's spline(x, y)",
             "R's spline(x, y, method='natural')",
             "R's approx(x, y, method='linear')",
             "R's approx(x, y, method='linear') (on log scale)",
             "Frequencies used for SII calculation"
             ),
           col=c("black", "red", "blue","cyan","green","orange","magenta","brown", "black"),
           pch=c(      1,    NA,    NA,    NA,     NA,      NA,       NA,      NA,     NA),
           lty=c(     NA,     1,     1,     1,      1,       1,        1,       1,      3),
           lwd=c(     NA,     3,     2,     2,      2,       2,        2,       2,      1),
           bg="white",
           cex=if(separate) 1.25 else 0.75,
           ncol=if(separate)   2 else 1
           )
    if(separate)
      layout( 1 )
    title(paste("Spline method comparison\n", if(logx)"(log scale)"else"(natural scale)", "\n"))
    
  }


###################################################
### code chunk number 26: natural_scale_one_plot
###################################################
doplot(FALSE, FALSE)


###################################################
### code chunk number 27: natural_scale_separate_plots
###################################################
doplot(FALSE, TRUE)


###################################################
### code chunk number 28: log_scale_one_plot
###################################################
doplot(TRUE, FALSE)


###################################################
### code chunk number 29: log_scale_separate_plots
###################################################
doplot(TRUE, TRUE)


###################################################
### code chunk number 30: SII:::sii.excel
###################################################
args(SII:::sii.excel)


###################################################
### code chunk number 31: sii
###################################################
args(sii)


###################################################
### code chunk number 32: sii.print
###################################################



###################################################
### code chunk number 33: show.reload.constants
###################################################
SII:::reload.constants 


###################################################
### code chunk number 34: run.reload.constants (eval = FALSE)
###################################################
## SII:::reload.constants(xls.path="./SII/extdata")


###################################################
### code chunk number 35: ComparisonTable1
###################################################
sii.left <- sii(
                speech="raised",
                threshold=c(25,25,30,35,45,45,55,60),
                freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                importance="NU6",
                interpolate=TRUE
                )


## comparison of our interpolation and matlab's
tab <- cbind(
             matlab=SII:::sii.constants[,"Ti'.THDN"],
             R     =sii.left$table[,"T'i"],
             delta =NA
             )
tab[,3] <- tab[,1] - tab[,2]
rownames(tab) <- SII:::sii.constants[,".NFreqLin"]
t(round(tab,2))


###################################################
### code chunk number 36: <plot.sii
###################################################
compare.plot <- function(x, matlab, title)
  {
    plot(x)
    lines(SII:::sii.constants[,".NFreqLin"],
          matlab,
          type="l", col="red", lwd=3)

    legend("topleft",
           legend=c(
             "Measured data",
             "Matlab INTERP1(X, Y, XI, 'spline', 'extrap')",
             "R's approx(X,Y, XI, xout=XI,\n method='linear', rule=2)"  
             ),
           col=c("black", "red","blue","green","orange","magenta"),
           pch=c(      1,    NA,    2,     NA,      NA,       NA),
           lty=c(     NA,     1,    1,      1,       1,        1),
           lwd=c(     NA,     3,    2,      2,       2,        2),
           bg="white"
           )
    title(title)
  }


###################################################
### code chunk number 37: ComparisonFigure1
###################################################

compare.plot(sii.left, matlab=SII:::sii.constants[, "Ti'.THDN"], title="Spline method comparison, Left Ear")
      


###################################################
### code chunk number 38: ComparisonTable2
###################################################
# comparison of our interpolation and matlab's

matlab <- c(15.00, 15.00, 13.34, 14.27, 16.06, 17.82, 19.18, 20.00,
            20.26, 20.54, 21.44, 23.36, 26.93, 31.38, 34.63, 35.13, 
            35.00, 38.29, 43.98, 48.80, 49.58)

sii.right <- sii( 
                 speech="raised",
                 threshold=c(15,15,20,25,35,35,45,50), 
                 freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 importance="NU6",
                 interpolate=TRUE
                 )

tab <- cbind(
             matlab=matlab,
             R     =sii.right$table    [,"T'i"],
             delta =NA
             )
tab[,3] <- tab[,1] - tab[,2]
rownames(tab) <- SII:::sii.constants[,".NFreqLin"]
t(round(tab,2))


###################################################
### code chunk number 39: ComparisonFigure2
###################################################
compare.plot(sii.right, matlab=matlab, title="Spline method comparison, Right Ear")


###################################################
### code chunk number 40: <sii.left.old
###################################################
SII:::sii.excel( 
          c(25,25,30,35,45,45,55,60),
          c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000) 
          )


###################################################
### code chunk number 41: sii.left
###################################################
sii.left <- sii(
                speech="raised",
                threshold=c(25,25,30,35,45,45,55,60),
                freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                importance="NU6",
                interpolate=TRUE
                )
sii.left


###################################################
### code chunk number 42: <sii.right.old
###################################################
SII:::sii.excel( 
          c(15,15,20,25,35,35,45,50), 
          c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000) 
          )


###################################################
### code chunk number 43: sii.right
###################################################
sii.right <- sii( 
                 speech="raised",
                 threshold=c(15,15,20,25,35,35,45,50), 
                 freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 importance="NU6",
                 interpolate=TRUE
                 )
sii.right


###################################################
### code chunk number 44: <sii.best.old
###################################################
SII:::sii.excel( 
          rep(0,8),
          freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000) 
          )


###################################################
### code chunk number 45: sii.best
###################################################
sii.best <- sii( 
                threshold=rep(0,8),
                freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                interpolate=TRUE 
                )
sii.best


###################################################
### code chunk number 46: <sii.worst.old
###################################################
SII:::sii.excel( 
          rep(100,8),
          c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000) 
          )


###################################################
### code chunk number 47: sii.worst
###################################################
sii.worst <- sii( 
                 threshold=rep(100,8),
                 freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 interpolate=TRUE
                )
sii.worst


###################################################
### code chunk number 48: sii.missing
###################################################
sii.worst <- sii( 
                 threshold=c(NA, rep(100,7)),
                 freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 interpolate=TRUE
                )
sii.worst


###################################################
### code chunk number 49: sii.missing
###################################################
sii.right <- sii( 
                 speech="raised",
                 threshold=c(0,15,15,20,25,35,35,45,50), 
                 freq=c(NA, 250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 importance="NU6",
                 interpolate=TRUE
                 )
sii.right


###################################################
### code chunk number 50: sii.all.missing
###################################################
## This should fail, because there is no data!
sii.NONE <- try(
             sii( 
                 threshold=rep(NA,8),
                 freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 interpolate=TRUE
                )
            )
sii.NONE


###################################################
### code chunk number 51: sii.missing
###################################################
sii.right <- sii( 
                 speech="raised",
                 threshold=c(15,15,20,NA,35,35,45,50), 
                 freq=c(250, 500, 1000, 2000, 3000, 4000, 6000, 8000),
                 importance="NU6",
                 interpolate=TRUE
                 )
sii.right


###################################################
### code chunk number 52: missing.C.1
###################################################
sii.C1.NA <- sii(
              speech   = c(50.0, 40.0, 40.0, NA,   20.0,  0.0),
              noise    = c(70.0, 65.0, 45.0, 25.0,  1.0,-15.0),
              threshold= c( 0.0,  0.0,  0.0,  0.0,   NA,  0.0),
              freq     = c( 250,  500, 1000, 2000, 4000, 8000),
              method="octave",
              importance="CST",
              interpolate=TRUE
              )
sii.C1
sii.C1.NA


###################################################
### code chunk number 53: save data
###################################################
save.image("SII-Code.Rda")

Try the SII package in your browser

Any scripts or data that you put into this service are public.

SII documentation built on May 2, 2019, 6:47 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close