1  % 'FFT2_detrend': calculate the 2D spectrum of the input scalar after removing the linear trend


2  %


3  %%%% Use the general syntax for transform fields with a single input %%%%


4  % OUTPUT:


5  % DataOut: output field structure


6  %


7  %INPUT:


8  % DataIn: first input field structure


9 


10  %=======================================================================


11  % Copyright 20082016, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


12  % http://www.legi.grenobleinp.fr


13  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


14  %


15  % This file is part of the toolbox UVMAT.


16  %


17  % UVMAT is free software; you can redistribute it and/or modify


18  % it under the terms of the GNU General Public License as published


19  % by the Free Software Foundation; either version 2 of the license,


20  % or (at your option) any later version.


21  %


22  % UVMAT is distributed in the hope that it will be useful,


23  % but WITHOUT ANY WARRANTY; without even the implied warranty of


24  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


25  % GNU General Public License (see LICENSE.txt) for more details.


26  %=======================================================================


27 


28  function DataOut=FFT2_detrend(DataIn)


29  %


30  DataOut=[];


31  if strcmp(DataIn,'*')


32  DataOut.InputFieldType='scalar';


33  return


34  end


35  %%%%%%%%%%%%%%%%%%%%%%%%%


36  [CellInfo,NbDim,errormsg]=find_field_cells(DataIn);


37  if ~isempty(errormsg)


38  DataOut.Txt=errormsg;


39  return


40  end


41  DataOut.ListVarName={};


42  DataOut.VarDimName={};


43  for ilist=1:numel(CellInfo)


44  if NbDim(ilist)==2 && numel(CellInfo{ilist}.CoordIndex)==2 % field with structured coordinates


45  %process coordinates


46  CoordName=DataIn.ListVarName(CellInfo{ilist}.CoordIndex);


47  x1 = DataIn.(CoordName{2}); y1 = DataIn.(CoordName{1});


48  [x y] = meshgrid(x1,y1);


49  coeff(1,1) = sum(sum(x.^2)); coeff(1,2) = sum(sum(x.*y)); coeff(1,3) = sum(sum(x));


50  coeff(2,1) = sum(sum(x.*y)); coeff(2,2) = sum(sum(y.^2)); coeff(2,3) = sum(sum(y));


51  coeff(3,1) = sum(sum(x)); coeff(3,2) = sum(sum(y)); coeff(3,3) = length(x1)*length(y1);


52  delta_x = x1(2)  x1(1); delta_y = y1(2)  y1(1);


53  Nx = length(x1); Ny = length(y1);


54  Nxa = 1:Nx; Nya = 1:Ny;


55  ssx = find(Nxa<Nx/2); ssy = find(Nya<Ny/2);


56  Nxa(Nxssx+1) = Nxa(ssx)+1; Nya(Nyssy+1) = Nya(ssy)+1;


57  [Nxa Ix] = sort(Nxa); [Nya Iy] = sort(Nya);


58  kx1 = (2*pi/delta_x/Nx)*(Nxa1); ky1 = (2*pi/delta_y/Ny)*(Nya1);


59  ss = find(ky1>=0); ky1 = ky1(ss);


60  [kx ky] = meshgrid(kx1,ky1);


61  DataOut.(CoordName{2}) = kx1; DataOut.(CoordName{1}) = ky1;


62  if isfield(DataIn,'CoordUnit')


63  DataOut.CoordUnit=[DataIn.CoordUnit '^{1}'];


64  DataOut.ListGlobalAttribute={'CoordUnit'};


65  end


66  %process scalar


67  ivar=CellInfo{ilist}.VarIndex_scalar(1);


68  VarName=DataIn.ListVarName{ivar};


69  z=DataIn.(VarName);


70  rhs(1) = sum(sum(x.*z)); rhs(2) = sum(sum(y.*z)); rhs(3) = sum(sum(z));


71  lin_coeff = inv(coeff)*rhs';


72  lin_trend = lin_coeff(1)*x + lin_coeff(2)*y + lin_coeff(3);


73  z2 = z  lin_trend;


74  spec2 = abs(fft2(z2)).^2;


75  spec2 = spec2(Iy,Ix);


76  spec2 = spec2(ss,:);


77  %DataOut.(VarName) = log(spec2);


78  DataOut.(VarName) = spec2;


79  spec_sum=sum(sum(spec2));


80  kx_mean=sum(sum(spec2.*kx))/spec_sum;


81  ky_mean=sum(sum(spec2.*ky))/spec_sum;


82  theta=atand(ky_mean/kx_mean);


83  lambda=2*pi/(sqrt(kx_mean*kx_mean+ky_mean*ky_mean));


84 


85  DataOut.ListVarName=[CoordName {VarName}];%list of variables


86  DataOut.VarDimName=[CoordName {CoordName}];%list of dimensions for variables


87  break


88  end


89  end


90 


91 


92 

