#definitions and utilities for ex. on fts
 1    ;
function d=delta(n,i) 
  d=zeros(1,n);
  d(i)=1;
end

function s=sqr(x) 
  s=x .* x;
end

function m=wmean(y,w)
  m=sum(y .* w)/ sum(w);
end

function view_ft(N,dspan,dsigma)
  sup=linspace(0,dspan*(N-1)/N,N);
  ispan=N/dspan;
  center=dspan/2;
  filter=menu("scegli una distribuzione","gaussiana","lorenziana","esponenziale","gaussiana in 0");
  switch(filter)
    case 1
      dir_y=normpdf(sup,center,dsigma);
    case 2
      dir_y=dsigma*dsigma./(dsigma*dsigma+(sup-center).*(sup-center));
    case 3
      dir_y=exp(-abs(sup-center)/dsigma);
    case 4
      dir_y=normpdf(sup,0,dsigma)+normpdf(sup,dspan,dsigma);
  endswitch

  subplot(2,2,1);
  plot(sup,dir_y,';diretta;');
  isup=linspace(0,ispan*(N-1)/N,N);
  inv_y=fft(dir_y);
  subplot(2,2,2);
  plot(isup,inv_y,';inversa;',isup,imag(inv_y),';im inversa;');
  xlim([0,ispan]);  
  subplot(2,2,3);
  dir_y2=dir_y .* conj(dir_y);
  inv_y2=inv_y .* conj(inv_y);
  plot(sup,dir_y2,';diretta^2;');
  subplot(2,2,4);
  plot(isup(1:N/2),inv_y2(1:N/2),isup(N/2+1:N)-ispan,inv_y2(N/2+1:N),
       ';inversa^2;');
  xlim([-ispan/2,ispan/2]);
  m_y=wmean(sup,dir_y2);
  s_y=wmean(sqr(sup-m_y),dir_y2);
  printf("la media e la sigma di y^2 : %f, %f\n", m_y, sqrt(s_y));
  isup(N/2+1:N)=isup(N/2+1:N)-ispan;
  im_y=wmean(isup,inv_y2);
  is_y=wmean(sqr(isup-im_y),inv_y2);
  printf("la media e la sigma di ~y^2 : %f, %f\n", im_y, sqrt(is_y));
end

global speclen=8192;

function [spec,x]=sim_spectrum(npeak,width)
  global speclen
  spec=zeros(1,speclen);
  x=linspace(0,speclen-1,speclen);
  left=speclen*0.2;
  right=speclen*0.8;
  if (npeak==1)
    spacing=right-left;
  else
    spacing=(right-left)/(npeak-1);
  endif
#  width=min(spacing,left);
  for i=0:npeak-1
    spec=spec+cauchy_pdf(x,left+spacing*i,width*0.3);
  endfor
endfunction

function [sp,xx]=ex_sample()
  [spec,x]=sim_spectrum(40,3);
  speclen=length(x);
  plot(x,spec);
  input("cont","s");
  spect=fftshift(fft(spec));
  plot(x,spect);
  specu=spect(floor(speclen*0.1):ceil(speclen*0.9));
  xu=x(floor(speclen*0.1):ceil(speclen*0.9));
  input("cont","s");
  plot(xu,specu);
  input("cont","s");
  speca=ifft(fftshift(specu));
  plot(xu,speca);
  filter=menu("scegli un filtro","gstrong","gweak","triang","raised_cos");
  c_xu=mean(xu);
  switch(filter)
    case 1
      specw=specu.*normpdf(xu,c_xu,speclen*0.1);
    case 2
      specw=specu.*normpdf(xu,c_xu,speclen*0.5);
    case 3
      specw=specu.*(1-abs(xu/c_xu-1));
    case 4
      scale=pi/(c_xu-xu(1));
      specw=specu.*(cos((xu-c_xu)*scale)+1);
  endswitch
  plot(xu,specw);
  input("cont","s");
  speca=ifft(fftshift(specw));
  plot(xu,speca);
  sp=speca;
  xx=xu;
endfunction

#[sp1,xx]=ex_sample();
#[sp2,~]=ex_sample();
#plot(xx,sp1/max(sp1),';sp1;',xx,sp2/max(sp2),';sp2;');

function [sp,xx]=make_spectrum(specname)
  fid=fopen(specname);
  spec=textscan(fid,'%f %f');
  spec=cell2struct(spec,{'x';'y'},2);
  lx=max(spec.x)-min(spec.x);
  minx=min(spec.x)-lx*0.1;
  lx=lx*1.2
  dx=lx/8192
  xx=(0:8191)*dx+minx;
  sp=zeros(1,8192);
  for i = 1:size(spec.x,1)
    ind=round((spec.x(i)-minx)/dx);
    sp(ind)=spec.y(i);
  endfor
end

function show_filter_effect(spec,x)
  speclen=length(x);
  input("lo spettro","s");
  plot(x,spec);
  speca=fft(spec);
  lx=speclen*(x(2)-x(1));
  maxk=pi/lx;
  xv=(0:8191)*maxk/4096;
  xu=fftshift(xv);
  xu(1:4096)=xu(1:4096)-2*maxk;
  #xu(1:128:8192)
  input("its fourier transform","s");
  plot(xv,speca);

  specu=fftshift(speca);
  input("its fourier transform recentered","s");
  plot(xu,specu);
  
  filter=menu("scegli un filtro","gstrong","gweak","triang","raised_cos");
  switch(filter)
    case 1
      wf=normpdf(xu,0,maxk*0.1);
    case 2
      wf=normpdf(xu,0,maxk*0.5);
    case 3
      wf=(1-abs(xu)/maxk);
    case 4
      scale=pi/maxk;
      wf=(cos(xu*scale)+1);
  endswitch
  specw=specu.*wf*8192/sum(wf); #normalization?
  input("filtered ft","s");
  plot(xu,specw);
  input("the original spectrum and the convoluted","s");
  specf=ifft(fftshift(specw));
  plot(x,abs(specf),x,spec);
endfunction

function [sp,xx,xx2]=make_spectrum_c13(specname)
  fid=fopen(specname);
  spec=textscan(fid,'%f %f');
  spec=cell2struct(spec,{'x';'y'},2);
  lx=max(spec.x)-min(spec.x);
  minx=min(spec.x)-lx*0.1;
  lx=lx*1.2;
  dx=lx/8192;
  xx=(0:8191)*dx+minx;
  sp=zeros(1,8192);
  for i = 1:size(spec.x,1)
    ind=round((spec.x(i)-minx)/dx);
    sp(ind)=spec.y(i);
  endfor
  mass12=12*16/(12+16)
  mass13=13*16/(13+16)
  cb=0;
  dp=spec.x(2)-spec.x(1);
  #search the band center
  for i=2:size(spec.x,1)
    dpt=spec.x(i)-spec.x(i-1);
    if (dpt<dp*1.6)
      dp=dpt;
    else
      cb=(spec.x(i)+spec.x(i-1))*0.5
      break
    endif
  endfor
  cb13=cb*sqrt(mass12/mass13);
  xx2=(xx-cb)*(mass12/mass13)+cb13;
end

  
