options(repr.plot.width=3.5, repr.plot.height=3.5)
options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"

library(Matrix)
library(tuneR)
library(colorRamps)
library(audio)

# Importing the libraries
for (f in list.files(path="nt_toolbox/toolbox_general/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_general/", f, sep=""))
}
for (f in list.files(path="nt_toolbox/toolbox_signal/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_signal/", f, sep=""))
}

n = 1024*16
s = 3 #number of sounds
p = 2 #number of micros

x = Matrix(0,nrow=n, ncol=3)
x[,1] = as.vector(load_sound("nt_toolbox/data/bird.wav", n))
x[,2] = as.vector(load_sound("nt_toolbox/data/female.wav", n))
x[,3] = as.vector(load_sound("nt_toolbox/data/male.wav", n))

x = x/matrix(rep(c(sd(x[,1]),sd(x[,2]),sd(x[,3])),each=n),nrow=n)

theta = seq(from=0, to =pi, by=pi/s)[-(s+1)]
theta[1] = 0.2
M = rbind(cos(theta),sin(theta))

y = x %*% t(M)

options(repr.plot.width=6, repr.plot.height=3.5)
for (i in (1:s)){
    plot(x[, i], main=paste("Source", i), type="l", ylab="", xlab="", col="blue")
    }

options(repr.plot.width=6, repr.plot.height=3.5)
for (i in (1:p)){
    plot(y[, i], main=paste("Micro", i), type="l", ylab="", xlab="", col="blue")
    }

w = 128   #size of the window
q = floor(w/4)  #overlap of the window

X = replicate(s, Matrix(0,nrow=w,ncol=4*w+1), simplify=FALSE)
Y = replicate(p, Matrix(0,nrow=w,ncol=4*w+1), simplify=FALSE)

for (i in (1:s))
{
    X[[i]] = perform_stft(x[,i], w,q,n)
    plot_spectogram(X[[i]], paste("Source", i))
}

source("nt_solutions/audio_separation/exo1.R")

## Insert your code here.

mf = dim(Y[[1]])[1]
mt = dim(Y[[1]])[2]
P = matrix(0, nrow=mt*mf, ncol=p)
c = c()
d = c()
for (j in (1:mf))
{
    c = c(c,Y[[1]][j,])
    d = c(d,Y[[2]][j,])
}
P[,1] = c
P[,2] = d
P = rbind(Re(P), Im(P))

npts = 6000

options(repr.plot.width=5, repr.plot.height=4)

sel = sample(n)

sel = sel[1:npts]

plot(y[sel,1], y[sel,2], type="p", pch=19, col = "blue", cex=0.05, main="Time domain", xlab="", ylab="")

source("nt_solutions/audio_separation/exo2.R")

## Insert your code here.

nrow = dim(P)[1]
Theta = c()
for (i in (1:nrow))
{
    Theta[i] = atan2(P[i,2], P[i,1])%%pi
}

options(repr.plot.width=4, repr.plot.height=4)

nbins = 100
t = seq(from=pi/200, pi, (pi-pi/200)/(nbins-1))
hist = hist(Theta, xlim=c(min(Theta),max(Theta)), main="", breaks=nbins, plot=FALSE)
h = hist$counts/sum(hist$counts)
barplot(h, xlab="Theta", col ="DarkBlue" , tick=TRUE)

source("nt_solutions/audio_separation/exo3.R")

## Insert your code here.

source("nt_solutions/audio_separation/exo4.R")

## Insert your code here.

A = matrix(0, nrow=mt*mf, ncol=p)
c = c()
d = c()
for (j in (1:mf))
{
    c = c(c,Y[[1]][j,])
    d = c(d,Y[[2]][j,])
}
A[,1] = c
A[,2] = d

C = abs(t(M) %*% t(A))

vec = c(1:dim(C)[2])
tmp = apply(C[,vec],2,max)
I = max.col(t(C))
I = t(matrix(I, nrow=mt, ncol=mf))

T = .05
D = sqrt(abs(Y[[1]])**2+abs(Y[[2]])**2)
I = I*(D > T)

options(repr.plot.width=8, repr.plot.height=2.5)
options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"

image(t(I[1:floor(mf/2),]),col=matlab.like2(50),interpolation_type="nearest")

Proj = t(M) %*% t(A)
Xr = replicate(s, matrix(0,nrow=w,ncol=4*w+1), simplify=FALSE)
for (i in (1:s))
{
    Xr[[i]] = t(matrix(Proj[i,], nrow=mt, ncol=mf))*((I+1)==i)
}

xr = Matrix(0, nrow=n, ncol=s)
for (i in (1:s))
{
    xr[,i] = perform_stft(Re(Xr[[i]]), w,q,n)
}


options(repr.plot.width=5, repr.plot.height=3)
for (i in (1:s))
{
    plot(xr[,i], main=paste("Estimated source #", i), type="l", ylab="", xlab="", col="blue")
}

i = 3
play(x[,i], rate=15000)

play(xr[,i] * 1e4, rate=15000)

<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>