pi=%pi
e=%e
i=%i

// dot product for functions
// f and g must be of the same size
function out=func_dot(f,g)
    //TODO check sizes
    out = f*g'
endfunction

// returns the norm of f
function out=func_norm(f,dot)
    out = sqrt(dot(f,f));
endfunction

function out = project(f,g,dot)
    out = dot(f,g)/dot(g,g)*g
endfunction

function out = reject(f,g,dot)
    out = f-project(f,g,dot)
endfunction

function out = func_rotate(func,f,g,dot,theta)
    ff = f/func_norm(f,dot)
    g = reject(g,f)
    gg = g/func_norm(g,dot)
    // TODO : check these next lines
    out = reject(reject(func,ff,dot),gg,dot) +...
          dot(func,ff)*( cos(theta).*ff + sin(theta).*gg )+...
          dot(func,gg)*(-sin(theta).*ff + cos(theta).*gg )
endfunction

snd1 = loadwave("/home/doek/Downloads/Roland-JV-2080-PizzicatoStr-C5.wav")
snd2 = loadwave("/home/doek/Downloads/Roland-JV-2080-SynVox-C4.wav")
snd3 = loadwave("/home/doek/Downloads/Alesis-Fusion-Bright-Acoustic-Piano-C4.wav")
snd4 = loadwave("/home/doek/Downloads/a.wav")
snd5 = loadwave("/home/doek/Downloads/o.wav")

// one channel, same length
snd1 = snd1(1,1:50000)
snd2 = snd2(1,1:50000)
snd3 = snd3(1,1:50000)
snd4 = snd4(1,1:50000)
//snd5 = snd5(1,1:50000)

//playsnd(snd1)
//playsnd(snd2)
//playsnd(snd3)
disp(func_dot(snd1,snd2))
disp(func_dot(snd2,snd3))
disp(func_dot(snd1,snd3))

//out = func_rotate(snd1,snd1,snd2,func_dot,(1:50000)/50000*2*2*pi)
//out = func_rotate(snd1,snd2,snd3,func_dot,pi/4)
//playsnd(out)

//rate = 44100
//snd1_fft = fft(snd1);
//snd1_t = (0:length(snd1)-1)*rate/length(snd1) - rate/2;
//clf()
//plot(snd1_t,abs(fftshift(snd1_fft)))
//title('Magnitude')

//samp=2/44100
//t=0:samp:1;
//f=sin(2*pi*440*t)
//playsnd(f)
//g=exp(t)/25
//clf()
//plot(t,f)
//plot(t,g)

samp=2/44100
t=0:samp:1
f1=sin(2*pi*4*t)
f2=exp(t)
f3=sin(1./(t+0.00001))*10
clf()
//out = func_rotate(f3,f1,f2,func_dot,t*2*2*pi)
//plot(t,out)

for tt = 0:0.01:1
    clf()
    out = func_rotate(f3,f1,f2,func_dot,tt*2*2*pi)
    plot(t,out)
end
