source: trunk/src/transform_field/FFT2_detrend.m @ 951

Last change on this file since 951 was 924, checked in by g7moreau, 8 years ago
  • Update Copyright Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
File size: 3.8 KB
Line 
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 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
12%   http://www.legi.grenoble-inp.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
28function DataOut=FFT2_detrend(DataIn)
29%------------------------------------------------------------------------
30DataOut=[];
31if strcmp(DataIn,'*')   
32    DataOut.InputFieldType='scalar';
33    return
34end
35%%%%%%%%%%%%%%%%%%%%%%%%%
36[CellInfo,NbDim,errormsg]=find_field_cells(DataIn);
37if ~isempty(errormsg)
38    DataOut.Txt=errormsg;
39    return
40end
41DataOut.ListVarName={};
42DataOut.VarDimName={};
43for 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(Nx-ssx+1) = -Nxa(ssx)+1; Nya(Ny-ssy+1) = -Nya(ssy)+1;
57        [Nxa Ix] = sort(Nxa); [Nya Iy] = sort(Nya);
58        kx1 = (2*pi/delta_x/Nx)*(Nxa-1); ky1 = (2*pi/delta_y/Ny)*(Nya-1);
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
89end
90 
91
92
Note: See TracBrowser for help on using the repository browser.