source: trunk/src/transform_field/FFT2_detrend.m

Last change on this file was 1127, checked in by g7moreau, 9 months ago

Update Joel email

File size: 4.0 KB
RevLine 
[574]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 %%%%
[510]4% OUTPUT:
[574]5% DataOut:   output field structure
[810]6%
[510]7%INPUT:
[574]8% DataIn:  first input field structure
[810]9
10%=======================================================================
[1126]11% Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
[810]12%   http://www.legi.grenoble-inp.fr
[1127]13%   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
[810]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
[510]28function DataOut=FFT2_detrend(DataIn)
[574]29%------------------------------------------------------------------------
[510]30DataOut=[];
31if strcmp(DataIn,'*')   
32    DataOut.InputFieldType='scalar';
33    return
34end
35%%%%%%%%%%%%%%%%%%%%%%%%%
[532]36[CellInfo,NbDim,errormsg]=find_field_cells(DataIn);
[510]37if ~isempty(errormsg)
38    DataOut.Txt=errormsg;
39    return
40end
41DataOut.ListVarName={};
42DataOut.VarDimName={};
[532]43for ilist=1:numel(CellInfo)
44    if NbDim(ilist)==2 && numel(CellInfo{ilist}.CoordIndex)==2 % field with structured coordinates
[510]45        %process coordinates
[532]46        CoordName=DataIn.ListVarName(CellInfo{ilist}.CoordIndex);
[510]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
[532]67        ivar=CellInfo{ilist}.VarIndex_scalar(1);
[510]68        VarName=DataIn.ListVarName{ivar};
69        z=DataIn.(VarName);
[1123]70        z(isnan(z))=0;% set to 0 NaN values
[510]71        rhs(1) = sum(sum(x.*z)); rhs(2) = sum(sum(y.*z)); rhs(3) = sum(sum(z));
72        lin_coeff = inv(coeff)*rhs';
73        lin_trend = lin_coeff(1)*x + lin_coeff(2)*y + lin_coeff(3);
74        z2 = z - lin_trend;
75        spec2 = abs(fft2(z2)).^2;
76        spec2 = spec2(Iy,Ix);
77        spec2 = spec2(ss,:);
[1123]78        DataOut.(VarName) = log10(spec2);
79        %DataOut.(VarName) = spec2;
[510]80        spec_sum=sum(sum(spec2));
[532]81        kx_mean=sum(sum(spec2.*kx))/spec_sum;
82        ky_mean=sum(sum(spec2.*ky))/spec_sum;
83        theta=atand(ky_mean/kx_mean);
84        lambda=2*pi/(sqrt(kx_mean*kx_mean+ky_mean*ky_mean));
[510]85       
86        DataOut.ListVarName=[CoordName {VarName}];%list of variables
87        DataOut.VarDimName=[CoordName {CoordName}];%list of dimensions for variables
[1113]88        DataOut.VarAttribute{1}.Role='coord_y';
89        DataOut.VarAttribute{2}.Role='coord_x';
90        DataOut.VarAttribute{3}.Role='scalar';
[510]91        break
92    end
93end
94 
95
96
Note: See TracBrowser for help on using the repository browser.