| 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(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
|
|---|
| 72 | end
|
|---|
| 73 |
|
|---|
| 74 |
|
|---|
| 75 |
|
|---|
