#=======================
drive <- "C:"
code.dir <- paste(drive, "DataMining", "Code", sep="/")
data.dir <- paste(drive, "DataMining", "Data", sep="/")
#
source(paste(code.dir, "borderhist.r", sep="/"))
source(paste(code.dir, "waveio.r", sep="/"))
   #
library(fastICA)
   #
   # Create and display two signals
S.1 <- sin((1:1000)/20) 
S.2 <- rep((((1:200)-100)/100), 5)
S <- cbind(S.1, S.2)
plot(S.1)
plot(S.2)

   #========================================
   # Mix the two signals.
   # Normally at random, but for this example rotate.
a <- pi/4
A <- matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
X <- S%*%A
plot(X[,1])
plot(X[,2])

border.hist(S.1, S.2)
border.hist(X[,1], X[,2])

   # Look at what happens as we rotate back.
b <- pi/36
W <- matrix(c(cos(b), -sin(b), sin(b), cos(b)), 2, 2)
   #
XX <- X
   #
graphics.off()
for (i in 1:9) {
   XX <- XX%*%W
   border.hist(XX[,1], XX[,2])
   readline("Press Enter...")
}


plot(XX[,1])
plot(XX[,2])

     
   #========================================
S.1 <- sin((1:1000)/20) 
S.2 <- rep((((1:200)-100)/100), 5)
S.3 <- rep(c(exp(seq(0,.99,.01))-1.845617, -exp(seq(0,.99,.01))+1.845617), 5)
   #
S <- cbind(S.1, S.2, S.3)
   #
A <- matrix(runif(9), 3, 3)        # Set a random mixing
X <- S%*%A
   #
a <- fastICA(X, 3, alg.typ = "parallel", fun = "logcosh", alpha = 1, 
             method = "R", row.norm = FALSE, maxit = 200, 
             tol = 0.0001, verbose = TRUE)
oldpar <- par(mfcol = c(3, 3), mar=c(2, 2, 2, 1))
   #
plot(1:1000, S[,1], type = "l", main = "Original Signals", xlab = "", ylab = "")
for (i in 2:3) {
   plot(1:1000, S[,i ], type = "l", xlab = "", ylab = "")
}
   #
plot(1:1000, X[,1 ], type = "l", main = "Mixed Signals", xlab = "", ylab = "")
for (i in 2:3) {
   plot(1:1000, X[,i], type = "l", xlab = "", ylab = "")
}   
   #
plot(1:1000, a$S[,1 ], type = "l", main = "ICA source estimates", xlab = "", ylab = "")
for (i in 2:3) {
   plot(1:1000, a$S[,i], type = "l", xlab = "", ylab = "")
}   
par(oldpar)

   #==============================
S.1 <- sin((1:1000)/20) 
S.2 <- rep((((1:200)-100)/100), 5)
s.3 <- tan(seq(-pi/2+.1,pi/2-.1,.0118))
S.3 <- rep(s.3, 4)
S.4 <- rep(c(exp(seq(0,.99,.01))-1.845617, -exp(seq(0,.99,.01))+1.845617), 5)
   #
S <- cbind(S.1, S.2, S.3, S.4)
   #
(A <- matrix(runif(16), 4, 4))     # Mix
X <- S%*%A
   #
a <- fastICA(X, 4, alg.typ = "parallel", fun = "logcosh", alpha = 1, 
             method = "R", row.norm = FALSE, maxit = 200, 
             tol = 0.0001, verbose = TRUE)

oldpar <- par(mfcol = c(4, 3), mar=c(2, 2, 2, 1))
plot(1:1000, S[,1], type = "l", main = "Original Signals", xlab = "", ylab = "")
for (i in 2:4) {
   plot(1:1000, S[,i ], type = "l", xlab = "", ylab = "")
}
   #
plot(1:1000, X[,1 ], type = "l", main = "Mixed Signals", xlab = "", ylab = "")
for (i in 2:4) {
   plot(1:1000, X[,i], type = "l", xlab = "", ylab = "")
}   
   #
plot(1:1000, a$S[,1 ], type = "l", main = "ICA source estimates", xlab = "", ylab = "")
for (i in 2:4) {
   plot(1:1000, a$S[,i], type = "l", xlab = "", ylab = "")
}   
par(oldpar)

   #========== Too few mixed ============
(A <- matrix(runif(12), 4, 3))
X <- S%*%A
   #
a <- fastICA(X, 4, alg.typ = "parallel", fun = "logcosh", alpha = 1, 
             method = "R", row.norm = FALSE, maxit = 200, 
             tol = 0.0001, verbose = TRUE)

oldpar <- par(mfcol = c(4, 3), mar=c(2, 2, 2, 1))
plot(1:1000, S[,1], type = "l", main = "Original Signals", xlab = "", ylab = "")
for (i in 2:4) {
   plot(1:1000, S[,i ], type = "l", xlab = "", ylab = "")
}
   #
plot(1:1000, X[,1 ], type = "l", main = "Mixed Signals", xlab = "", ylab = "")
for (i in 2:3) {
   plot(1:1000, X[,i], type = "l", xlab = "", ylab = "")
}   
plot(0, type="n")      # Dummy to fill
plot(1:1000, a$S[,1 ], type = "l", main = "ICA source estimates", xlab = "", ylab = "")
for (i in 2:3) {
   plot(1:1000, a$S[,i], type = "l", xlab = "", ylab = "")
}   
plot(0, type="n")      # Dummy to fill
par(oldpar)

   #========= Cocktail Party ==================
numb.source <- 9
in.file <- matrix(0, numb.source, 1)
mix.file <- matrix(0, numb.source, 1)
out.file <- matrix(0, numb.source, 1)
for (i in 1:numb.source) {
   in.file[i,] <-  paste(data.dir, "/source",i,".wav", sep="")
   mix.file[i,] <- paste(data.dir, "/m",i,".wav", sep="")
   out.file[i,] <- paste(data.dir, "/s",i,".wav", sep="")
}   

in.wav <- {}
for (m in 1:numb.source) {
   in.wav <- c(in.wav, list(read.wav(in.file[m,])))
}

wav.char <- function(wav) {
   cat("RIFF               = ", wav$RIFF, "\n")
   cat("Length             = ", wav$File.Len, "\n")
   cat("Wave               = ", wav$WAVE, "\n")
   cat("Format             = ", wav$format, "\n")
   cat("Format Length      = ", wav$len.of.format, "\n")
   cat("One                = ", wav$f.one, "\n")
   cat("Number of Channels = ", wav$Channel.numbs, "\n")
   cat("Sample Rate        = ", wav$Sample.Rate, "\n")
   cat("Bytes/Sec          = ", wav$Bytes.P.Sec, "\n")
   cat("Bytes/Sample       = ", wav$Bytes.P.Sample, "\n")
   cat("Bits/Sample        = ", wav$Bits.P.Sample, "\n")
   cat("Data               = ", wav$DATA, "\n")
   cat("Data Length        = ", wav$data.len, "\n")
}

wav.char(in.wav[[1]])
   # =============== Mix them =================
(A <- matrix(runif(numb.source*numb.source),numb.source,numb.source))
   #
mixed <- {}
for (i in 1:numb.source) {
   mixed <- cbind(mixed, in.wav[[i]]$data)
}
mixed <- mixed%*%A

old.par <- par(mfcol = c(numb.source, 1))
par(mar=c(2,2,1,2)+0.1)
plot(mixed[,1], type="l", main="Mixed")
for (m in 2:numb.source) {
   plot(mixed[,m], type="l")
}   
#if (dev.cur()[[1]]!=1)  bringToTop(which=dev.cur())
par(old.par)

   # =================================
mix.wav <- {}
for (m in 1:numb.source) {
   mix.wav <- c(mix.wav, list(in.wav[[m]]))
}
   #
for (m in 1:numb.source) {
   mix.wav[[m]]$data <- mixed[,m]
   write.wav(mix.file[m,], mix.wav[[m]])
}
   # =============== Play them =================
library(audio)
play(load.wave(mix.file[1,]))
play(load.wave(mix.file[2,]))
play(load.wave(mix.file[3,]))
play(load.wave(mix.file[4,]))
play(load.wave(mix.file[5,]))
play(load.wave(mix.file[6,]))
play(load.wave(mix.file[7,]))
play(load.wave(mix.file[8,]))
play(load.wave(mix.file[9,]))


   # =============== Unmix them =================
mixed.all <- {}
for (i in 1:numb.source) {
   mixed.all <- cbind(mixed.all, mixed[,i])
}
   #
ICA.wavs <- fastICA(mixed.all,
             numb.source, alg.typ = "parallel", fun = "logcosh", alpha = 1, 
             method = "R", row.norm = FALSE, maxit = 200, 
             tol = 0.0001, verbose = TRUE)

   # =============== Save them =================
new.wav <- {}
for (m in 1:numb.source) {
   new.wav <- c(new.wav, list(in.wav[[m]]))
}
   #
for (m in 1:numb.source) {
   new.wav[[m]]$data <- 5*ICA.wavs$S[,m]
   write.wav(out.file[m,], new.wav[[m]])
}

   # =============== Play them =================
play(load.wave(out.file[1,]))
play(load.wave(out.file[2,]))
play(load.wave(out.file[3,]))
play(load.wave(out.file[4,]))
play(load.wave(out.file[5,]))
play(load.wave(out.file[6,]))
play(load.wave(out.file[7,]))
play(load.wave(out.file[8,]))
play(load.wave(out.file[9,]))

   # =============== Plot them =================
old.par <- par(mfcol = c(numb.source, 3))
par(mar=c(2,2,1,1)+0.1)
   #
plot(in.wav[[1]]$data, type="l", main="Original")
for (m in 2:numb.source) {
   plot(in.wav[[m]]$data, type="l")
}
plot(mixed[,1], type="l", main="Mixed")
for (m in 2:numb.source) {
   plot(mixed[,m], type="l")
}   
plot(ICA.wavs$S[,1], type="l", )
if (dev.cur()[[1]]!=1)  bringToTop(which=dev.cur())
for (m in 2:numb.source) {
   plot(ICA.wavs$S[,m], type="l")
}   
par(old.par)

   # =============== Original - play =================
play(load.wave(in.file[1,]))
play(load.wave(in.file[2,]))
play(load.wave(in.file[3,]))
play(load.wave(in.file[4,]))
play(load.wave(in.file[5,]))
play(load.wave(in.file[6,]))
play(load.wave(in.file[7,]))
play(load.wave(in.file[8,]))
play(load.wave(in.file[9,]))


   #=========== End ===========