// calcula espectrograma de seņal s tomando logaritmo de espectro de potencia
// de corto plazo sobre ventanas de hamming de ancho M separadas por d
function x = xfft(s,M,d)
  [lhs,rhs] = argn(0)
  if rhs < 2 then M = 512; end
  if rhs < 3 then d = M/2; end
  L = size(s,'c')
  H = .54-.46*cos(2*%pi*(0:M-1)/(M-1))
  x = zeros(M,floor((L-M)/d)+1)
  for i=1:size(x,'c')
    x(:,i) = log10(abs(fft(s((1:M)+i*d-d).*H,-1)).^2)'
  end

// procesa archivo de sonido f (wav)
function [x,s] = procwav(f)
  s = loadwave(f)
  x = xfft(s,256,128)
  x = x(1:128,:)

// filtro exponencial con factor a sobre secuencia x
function x = fltexp(x,a)
  v = x(:,1)
  for i=1:size(x,'c')
    v(:) = a*x(:,i)+(1-a)*v
    x(:,i) = v
  end

// filtro exponencial con factor a sobre secuencia x en orden inverso
function x = fltexpr(x,a)
  v = x(:,$)
  for i=size(x,'c'):-1:1
    v(:) = a*x(:,i)+(1-a)*v
    x(:,i) = v
  end

// filtro gaussiano con ventana de ancho 2*l+1 sobre secuencia x
function y = fltgauss(x,l)
  [f,c] = size(x)
  if f > 1
    for i=1:f; x(i,:) = fltgauss(x(i,:),l); end
    y = x
    return
  end
  h = 3*(-l:l)/l
  h = exp(-h.*h)
  h = h/sum(h)
  y = convol(h,[x(1)*ones(1,l) x x($)*ones(1,l)])
  y = y(2*l+1:$-2*l)

// normalizacion de promedio sobre secuencia de vectores x
function x = normp(x)
  x = x - sum(x,'c')/size(x,'c')*ones(x(1,:))

// normalizacion de longitud sobre secuencia de vectores x
function x = norml(x)
  x = x./(ones(x(:,1))*sqrt(sum(x.*x,'r')))

// transformacion de Karhunen-Loewe (componentes proncipales)
// x es una secuencia de vectores (matriz en que cada columna es un vector)
// t es la matriz de proyeccion.
// t*x es la proyeccion
// t(1:2,:)*x es una reduccion dimensional (a dos dimensiones)
// antes de calcular esta transformacion es conveniente
// normalizar x haciendo x = normp(x)
function t = klt(x)
  [t,v] = spec(x*x')
  t = t(:,$:-1:1)'

// calcula un histograma bidimensional para puntos (x,y)
function r = hist2d(x,y,m)
  p = []; for i=x; p = [p, [i*ones(y);y]]; end
  u = ones(p(1,:))
  r = zeros(u)
  for c=m; r = r + exp(-sum((p-c*u).^2,'r')*30); end
  r = matrix(r,length(y),-1)

// grafica curvas (filas) de x
function plt(x)
  plot2d((1:size(x,'c'))',x',rect=[0,min(x),size(x,'c')+1 max(x)])

// grafica una matriz con niveles de gris
function mplt(m)
  v = [0 1]
  [f,c] = size(m)
  m = m - min(m)
  m = floor(m/max(m)*31.9999)
  m = m + 33
  oldmap = xget("colormap")
  if size(oldmap,'r')==32
    cmap = graycolormap(32); cmap = cmap($:-1:1,:)
    //cmap = hotcolormap(32)
    xset("colormap",[oldmap;cmap])
  end
  plot2d(0,0,rect=[0 0 c+1 1],axesflag=0)
  Matplot1(m($:-1:1,:),[0.5 v(1) c+.5 v(2)])

// graba una matriz m en archivo de nombre f (string)
function msave(f,m)
  MAT = m
  save(f,MAT)

// lee matriz almacenada en archivo f
function m = mload(f)
  load(f)
  m = MAT

// funcion sigmoide logistica
function l = logit(v)
  l = (1+exp(-v)).^(-1)

// crea una nueva ventana grafica
function w = nw()
  l = winsid()
  if l==[]
    w = 0
  else
    w = (0:max(l)+1)*0
    w(l+1) = 1
    w = find(w==0)-1
    w = w(1)
  end
  xset("window",w)

//ejemplo de lo que se puede hacer
function ej()
  [x,s] = procwav("aeiou.wav")
  n = size(s,'c')
  nw()
  subplot(4,1,1); plt(s)
  subplot(4,1,2); mplt(x)
  e = sum(x,'r')
  x = norml(x)
  subplot(4,1,3); mplt(x)
  subplot(4,1,4); plt([e;-555*ones(e)])
  xx = normp(x(:,find(e>-555)))
  t = klt(xx)
  nw()
  subplot(4,1,1); mplt(xx)
  y = t(1:4,:)*xx
  subplot(4,1,2); plt(y)
  subplot(4,1,3); mplt(t(1:2,:)'*t(1:2,:)*xx)
  subplot(4,1,4); mplt(t(1:4,:)'*t(1:4,:)*xx)
  c = zeros(x(1,:));
  c(17:34) = 1; c(57:76) = 2; c(98:118) = 3; c(145:160) = 4; c(184:197) = 5
  z = t*x
  nw()
  for i=1:size(z,'c')
    plot2d(z(1,i),z(2,i),style=-c(i))
  end
  plot2d(0,zeros(1,5),frameflag=0,style=-[1:5],leg="a@e@i@o@u")
  [x,s,t,y,z,c] = return(x,s,t,y,z,c)

function x = flt2d(x)
  x = fltgauss(x',10)'-fltgauss(x',20)'
  x = fltgauss(x,10)