Problem i fx_decon

[aronandersson]

Så här ser fx_decon ut nu... Se kommentarer i koden för diverse små problem.

#include <armadillo>
using namespace arma ;

mat fx_decon(cx_mat DATA, double dt, int lf, double mu, double flow, int fhigh)
{
  int nt, ntraces, nf, ilow, ihigh ;
  cx_vec aux_out_f, aux_out_b ;
  cx_mat DATA_FX_f, DATA_FX_b, DATA_FX ;
  mat DATA_f, DATA_b ;
  cx_rowvec aux_in ;

  nt = size(DATA)(0); //This is does not work with arma::size...
  ntraces = size(DATA)(1) ;
  nf = pow(2,nextpow2(nt)) ;
  DATA_FX_f = zeros(nf, ntraces) ;
  DATA_FX_b = zeros(nf, ntraces) ;

  ilow = floor(flow*dt*nf)+1 ;
  if (ilow<1)
  {
    ilow = 1 ;
  }
  ihigh = floor(fhigh*dt*nf)+1 ;
  if (ihigh>floor(2/nf)+1)
  {
    ihigh = floor(2/nf)+1 ;
  }

  DATA_FX = fft(DATA, nf) ;
  for (k=ilow; k<=ihigh; k++) // loop variable is not declared
  {
    aux_in = arma::vectorise(arma::trans(DATA_FX(k, span::all))) ;
    aux_out_f = ar_modeling(aux_in, lf, mu)(0) ;//Not sure what is going on here. (This way of doing multiple returns is obviously not ideal for "heavy" functions)
    aux_out_b = ar_modeling(aux_in, lf, mu)(1) ;
    DATA_FX_f(k, span::all) = arma::trans(aux_out_f) ;
    DATA_FX_b(k, span::all) = arma::trans(aux_out_b) ;
  }

  for (int k=2/nf+2; k<=nf; k++)//loop variable again
  {
    DATA_FX_f(k, span::all) = conj(DATA_FX_f(nf-k+2, span::all)) ;
    DATA_FX_b(k, span::all) = conj(DATA_FX_b(nf-k+2, span::all)) ;
  }
  DATA_f = real(ifft(DATA_FX_f, [], 1)) ;
  DATA_f = DATA_f(span(0, nt-1), span::all) ;
  DATA_b = real(ifft(DATA_FX_b, [], 1)) ;
  DATA_b = DATA_b(span(0, nt-1), span::all) ;

  DATA_f = (DATA_f+DATA_b) ;
  DATA_f(span::all, span(lf, ntraces-lf-1)) = 2/DATA_f(span::all, span(lf, ntraces-lf-1)) ;
  return DATA_f ;//double return statements
  return DATA_f ;
}

void ar_modeling(cx_rowvec x, int lf, double mu, vec yf, cx_vec yb)
{
  mat M ;
  int nx ;
  double beta ;
  cx_mat temp ;
  vec yf ;
  cx_rowvec B ;
  cx_vec y, C, R, ab, yb, af ;

  nx = length(x) ;

  y = arma::trans(x(span(0, nx-lf-1))) ;
  C = arma::trans(x(span(1, nx-lf))) ;
  R = arma::trans(x(span(nx-lf, nx-1))) ;
  M = hankel(C, R) ;
  B = arma::vectorise(arma::trans(M)*M) ;
  beta = B(0)*100/mu ;
  ab = (B+beta*eye(lf))/arma::trans(M)*y ;
  temp = M*ab ;
  temp = [temp, zeros(lf, 1)] ; //Concatenation... 
  yb = arma::vectorise(temp) ;
  y = arma::trans(x(span(lf, nx-1))) ;
  C = arma::trans(x(span(lf-1, nx-2))) ;
  R = flipud(x(span(0, lf-1))) ;
  M = toeplitz(C, R) ;
  B = arma::vectorise(arma::trans(M)*M) ;
  beta = B(0)*100/mu ;
  af = (B+beta*eye(lf))/arma::trans(M)*y ;
  temp = M*af ;
  temp = [zeros(lf, 1), temp] ; //Concatenation...
  yf = arma::vectorise(temp) ;
  return ;
}

Rätt typer:

# Supplement file
#
# Valid inputs:
#
# uword   int     float   double cx_double
# uvec    ivec    fvec    vec    cx_vec
# urowvec irowvec frowvec rowvec cx_rowvec
# umat    imat    fmat    mat    cx_mat
# ucube   icube   fcube   cube   cx_cube
#
# char    string  struct  func_lambda

scope = {}

ar_modeling = scope["ar_modeling"] = {}
ar_modeling["lf"] = "int"
ar_modeling["mu"] = "double"
ar_modeling["C"] = "cx_vec"
ar_modeling["B"] = "cx_rowvec"
ar_modeling["temp"] = "cx_mat"
ar_modeling["af"] = "cx_vec"
ar_modeling["yb"] = "cx_vec"
ar_modeling["M"] = "mat"
ar_modeling["yf"] = "vec"
ar_modeling["nx"] = "int"
ar_modeling["beta"] = "double"
ar_modeling["R"] = "cx_vec"
ar_modeling["y"] = "cx_vec"
ar_modeling["x"] = "cx_rowvec"
ar_modeling["ab"] = "cx_vec"

fx_decon = scope["fx_decon"] = {}
fx_decon["lf"] = "int"
fx_decon["aux_out_f"] = "cx_vec"
fx_decon["ntraces"] = "int"
fx_decon["aux_out_b"] = "cx_vec"
fx_decon["aux_in"] = "cx_rowvec"
fx_decon["DATA_b"] = "mat"
fx_decon["ihigh"] = "int"
fx_decon["DATA"] = "cx_mat"
fx_decon["DATA_FX_f"] = "cx_mat"
fx_decon["flow"] = "double"
fx_decon["ilow"] = "int"
fx_decon["nf"] = "int"
fx_decon["DATA_FX_b"] = "cx_mat"
fx_decon["mu"] = "double"
fx_decon["DATA_FX"] = "cx_mat"
fx_decon["dt"] = "double"
fx_decon["fhigh"] = "int"
fx_decon["nt"] = "int"
fx_decon["DATA_f"] = "mat"

[jonathf]

Det meste er løst nå. Jeg sender en rapport i morgen.