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

Last change on this file since 574 was 574, checked in by sommeria, 11 years ago

clean the transform fcts

File size: 3.0 KB
Line
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%------------------------------------------------------------------------
11function DataOut=FFT2_detrend(DataIn)
12%------------------------------------------------------------------------
13DataOut=[];
14if strcmp(DataIn,'*')
15    DataOut.InputFieldType='scalar';
16    return
17end
18%%%%%%%%%%%%%%%%%%%%%%%%%
19[CellInfo,NbDim,errormsg]=find_field_cells(DataIn);
20if ~isempty(errormsg)
21    DataOut.Txt=errormsg;
22    return
23end
24DataOut.ListVarName={};
25DataOut.VarDimName={};
26for 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(Nx-ssx+1) = -Nxa(ssx)+1; Nya(Ny-ssy+1) = -Nya(ssy)+1;
40        [Nxa Ix] = sort(Nxa); [Nya Iy] = sort(Nya);
41        kx1 = (2*pi/delta_x/Nx)*(Nxa-1); ky1 = (2*pi/delta_y/Ny)*(Nya-1);
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
72end
73
74
75
Note: See TracBrowser for help on using the repository browser.