#!/usr/bin/ruby

=begin
Allan分散を求めるためのコード
=end

class Array
  def sample(n)
    res = []
    0.step(self.size - 1, n){|i| res << self[i]}
    return res
  end
  
  def sum
    self.inject(0.0){|temp, a| temp + a}
  end
  
  def mean
    self.sum / self.size
  end
  
  def sum_sum2
    self.inject([0.0, 0.0]){|temp, a|
      temp[0] += a
      temp[1] += a ** 2
      temp
    }
  end
  
  def mean_std
    sum, sum2 = self.sum_sum2
    mean = sum / self.size
    std = Math::sqrt((sum2 - sum * mean) / (self.size - 1))
    return [mean, std]
  end
  
  def std
    self.mean_std.last
  end
  
  def diff
    res = [self.first]
    self[1..-2].each{|item|
      res[-1] = item - res[-1]
      res << item
    }
    res[-1] = self.last - res[-1]
    return res
  end
end

=begin
移植元のMatlabコード

function [sig,sig2,osig,msig,tsig,tau]=avar(y,tau0)

s=[];
x=[];
n=length(y);
jj=floor( log((n-1)/3)/log(2) );

for j=0:jj
  
  m=2^j;
    
  tau(j+1)=m*tau0;
  D=zeros(1,n-m+1);
  for i=1:n-m+1
    D(i)=sum(y(i:i+m-1))/m;
  end
  
  %N sample 
  sig(j+1)=std(D(1:m:n-m+1));
  
  %AVAR 
  sig2(j+1)=sqrt(0.5*mean((diff(D(1:m:n-m+1)).^2)));
  
  %OVERAVAR
  z1=D(m+1:n+1-m);
  z2=D(1:n+1-2*m);
  u=sum((z1-z2).^2);
  osig(j+1)=sqrt(u/(n+1-2*m)/2);
  
  %MVAR
  u=zeros(1,n+2-3*m);
  for L=0:n+1-3*m
      z1=D(1+L:m+L);
      z2=D(1+m+L:2*m+L);
      u(L+1)=(sum(z2-z1))^2;
  end
  
  uu=mean(u);
  msig(j+1)=sqrt(uu/2)/m;
  
  %TVAR
  tsig(j+1)=tau(j+1)*msig(j+1)/sqrt(3);
end

=end

class Allan
  def Allan.avar(y, tau0)
    s = []
    x = []
    n = y.size
    
    sig = []
    sig2 = []
    osig = []
    msig = []
    tsig = []
    tau = []
    
    ((Math::log((n-1)/3) / Math::log(2)).floor + 1).times{|j|
      m = 2 ** j
    
      tau << m * tau0
      d = []
      (n-m+1).times{|i| d << y[i...(i+m)].sum / m}
      
      #N sample 
      sig << d[0..(n-m)].sample(m).std

      #AVAR 
      sig2 << Math::sqrt(d[0..(n-m)].sample(m).diff.collect!{|item|
          item ** 2
        }.mean / 2)

      #OVERAVAR
      z1 = d[m..(n-m)]
      z2 = d[0..(n-(2*m))]
      u = z1.zip(z2).collect!{|item|
          (item.first - item.last) ** 2
        }.sum
      osig << Math::sqrt(u / (n+1-(2*m)) / 2)

      #MVAR
      u = []
      (n+2-3*m).times{|l|
        z1 = d[l...(m+l)]
        z2 = d[(m+l)...((2*m)+l)]
        u << z2.zip(z1).collect!{|item|
            item.first - item.last
          }.sum ** 2
      }
      
      uu = u.mean
      msig << Math::sqrt(uu/2) / m;

      #TVAR
      tsig << tau.last * msig.last / Math::sqrt(3);
    }
    
    return [sig,sig2,osig,msig,tsig,tau]
  end
end

if $0 == __FILE__ then
  require 'test/unit'

  class TC_Allan < Test::Unit::TestCase
    def setup
      @sample_data = [
        0.918624068347636, 0.209475798478216, 0.692552834577734, 0.91048674362905, 
        0.662206614063334, 0.148297596073455, 0.938753014152008, 0.902320152806411, 
        0.119996139906719, 0.0947379081387505, 0.0982231953355901, 0.00660114349832464, 
        0.0339708079482285, 0.17183127621723, 0.892318418281166, 0.411675756564584, 
        0.138062931231765, 0.578273649321607, 0.129059928893944, 0.0213952352727891, 
        0.418986210364156, 0.430381019071647, 0.344252565078658, 0.197027736384632, 
        0.879710213130768, 0.783776835129653, 0.530268578010031, 0.499293191661211, 
        0.661504352713999, 0.855895806614516, 0.760074671439117, 0.902842235920218]
      @sample_tau0 = 0.1
    end
  
    def teardown
    end
  
    def test_avar
      res = Allan.avar(@sample_data, @sample_tau0)
      p res
      expected = [
          [0.3310,    0.3006,    0.2587,    0.2608], # sig
          [0.2673,    0.2732,    0.2058,    0.2596], # sig2
          [0.2673,    0.2299,    0.1888,    0.1667], # osig
          [0.2673,    0.1820,    0.1612,    0.0900], # msig
          [0.0154,    0.0210,    0.0372,    0.0416], # tsig
          [0.1000,    0.2000,    0.4000,    0.8000], # tau
        ]
      
      expected.zip(res).each{|v1s, v2s|
        v1s.zip(v2s).each{|v1, v2|
          assert_in_delta(v1, v2, 0.0005, message="")
        }
      }
    end
  end
end

=begin
a = [0.918624068347636, 0.209475798478216, 0.692552834577734, 0.91048674362905,...
0.662206614063334, 0.148297596073455, 0.938753014152008, 0.902320152806411,...
0.119996139906719, 0.0947379081387505, 0.0982231953355901, 0.00660114349832464,...
0.0339708079482285, 0.17183127621723, 0.892318418281166, 0.411675756564584,...
0.138062931231765, 0.578273649321607, 0.129059928893944, 0.0213952352727891,...
0.418986210364156, 0.430381019071647, 0.344252565078658, 0.197027736384632,...
0.879710213130768, 0.783776835129653, 0.530268578010031, 0.499293191661211,...
0.661504352713999, 0.855895806614516, 0.760074671439117, 0.902842235920218]'

[sig,sig2,osig,msig,tsig,tau] = avar(a, 0.1)

% sig  = 0.3310    0.3006    0.2587    0.2608
% sig2 = 0.2673    0.2732    0.2058    0.2596
% osig = 0.2673    0.2299    0.1888    0.1667
% msig = 0.2673    0.1820    0.1612    0.0900
% tsig = 0.0154    0.0210    0.0372    0.0416
% tau  = 0.1000    0.2000    0.4000    0.8000
=end