#### ---- All files in the web directory: use the website path + filename


#### ---- load in audio file
library(audio) ; source("audioFunctions.R") 

S<-load.wave("audio2019-07-22-04:00:02.wav") 

S<-memsProcess(S) 

play(S) 



#### ---- STFT 

A<-(stft(S[1,])$A + stft(S[2,])$A )/2 ; A<-A[2:(nrow(A)/2),]

Y<-log(A) 

par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
image(t(Y)) 



#### ---- SVD fit  

sY<-svd(Y) 

U<-sY$u
V<-sY$v
d<-sY$d 

par(mfrow=c(1,3),mar=c(3,3,1,1),mgp=c(1.75,.75,0)) 

plot(d) 
plot(U[,1]) 
plot(V[,1]) 

Ysvd2<-U[,1:2]%*%diag(d[1:2])%*%t(V[,1:2]) 
sum( (Y-Ysvd2)^2 ) 


#### ---- additive fit 
mu<-mean(Y)
a<-rowMeans(Y)-mu 
b<-colMeans(Y)-mu 

Yadd<-mu+outer(a,b,"+")
sum( (Y-Yadd)^2 )

sum(dim(Y))

image(t(Y)) ; image(t(Ysvd2)) ; image(t(Yadd)) 



#### ---- AMMI fit 
R<-Y-Yadd
sR<-svd(R)

U<-sR$u
V<-sR$v
d<-sR$d 

plot(d) 

par(mfrow=c(2,2)) 
plot(U[,1]) 
plot(V[,1]) 
plot(U[,2])    
plot(V[,2])   

Yammi2<-Yadd + U[,1:2]%*%diag(d[1:2])%*%t(V[,1:2])

sum( (Y-Yammi2)^2 )

par(mfrow=c(2,2)) 

image(t(Y)) ; image(t(Ysvd2)) ; image(t(Yadd)) ; image(t(Yammi2) )



#### ----  comparisons 

## rank 1 vs additive 
sum( (Y- sY$d[1]*outer(sY$u[,1],sY$v[,1]) )^2 ) 

sum( (Y-Yadd)^2 )

rssSVD<-(sum(Y^2) - cumsum(sY$d^2))[1:10] 

rssAMMI<-(sum(R^2) - c(0,cumsum(sR$d^2)))[1:10] 

plot(rssSVD,ylim=range(c(rssSVD,rssAMMI)),col="blue",pch=16,cex=2) 
points(rssAMMI,col="green",pch=16) 



#### ---- SVD not always better
n<-10 ; p<-8

mu<-rnorm(1) 
a<-rnorm(n) 
b<-rnorm(p) 

M<-mu + outer(a,b,"+") 

Y<-M+matrix(rnorm(n*p),n,p)/2

## additive 
muhat<-mean(Y) 
ahat<-rowMeans(Y)-muhat
bhat<-colMeans(Y)-muhat

Madd<-muhat+outer(ahat,bhat,"+") 

sum( (Y-Madd)^2 ) 

## best rank 1  
sY<-svd(Y)
Msvd<-sY$d[1]*outer(sY$u[,1],sY$v[,1],"*") 

sum( (Y-Msvd)^2 )

plot(Y,Msvd,ylim=range(c(Msvd,Madd))) ; abline(0,1) 
points(Y,Madd,col="green",pch=16 ) 


## best rank 2 
Msvd<-sY$u[,1:2]%*%diag(sY$d[1:2])%*%t(sY$v[,1:2]) 
sum( (Y-Msvd)^2 )