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


2 


3  %


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


5  % OUTPUT:


6  % DataOut: output field structure


7 


8  %INPUT:


9  % DataIn: first input field structure


10  %


11  function DataOut=FFT2_detrend(DataIn)


12  %


13  DataOut=[];


14  if strcmp(DataIn,'*')


15  DataOut.InputFieldType='scalar';


16  return


17  end


18  %%%%%%%%%%%%%%%%%%%%%%%%%


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


20  if ~isempty(errormsg)


21  DataOut.Txt=errormsg;


22  return


23  end


24  DataOut.ListVarName={};


25  DataOut.VarDimName={};


26  for ilist=1:numel(CellInfo)


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


28  %process coordinates


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


45  if isfield(DataIn,'CoordUnit')


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


47  DataOut.ListGlobalAttribute={'CoordUnit'};


48  end


49  %process scalar


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


51  VarName=DataIn.ListVarName{ivar};


52  z=DataIn.(VarName);


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


54  lin_coeff = inv(coeff)*rhs';


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


56  z2 = z  lin_trend;


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


58  spec2 = spec2(Iy,Ix);


59  spec2 = spec2(ss,:);


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


61  DataOut.(VarName) = spec2;


62  spec_sum=sum(sum(spec2));


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


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


65  theta=atand(ky_mean/kx_mean);


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


67 


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


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


70  break


71  end


72  end


73 


74 


75 

