Home > atmlab > randomize > iaaft > iaaft_loop_2d_horizontal.m

iaaft_loop_2d_horizontal

PURPOSE ^

This function is the main loop of the Iterative Amplitude Adapted Fourier

SYNOPSIS ^

function [surrogate, error_amplitude, error_spec] = iaaft_loop_2d(coeff_2d,sorted_values)

DESCRIPTION ^

 This function is the main loop of the Iterative Amplitude Adapted Fourier
 Transform  (IAAFT) method to make surrogate fields. The method was developped by
 Schreiber and Schmitz, see e.g. Phys. Rev Lett. 77, pp. 635-, 1996, for 
 statistical non-linearity tests for time series.
 This method makes fields that have a specified amplitude distribution and
 power spectral coefficients. This can be seen as a field made with linear
 dynamics as seen through a static non-linear filter. It works by
 iteratively addaption the amplitude distribution and the Fourier
 coefficients (the phases are not changed in this step).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

iaaft_loop_2d_horizontal.m

SOURCE CODE ^

0001 function [surrogate, error_amplitude, error_spec] = iaaft_loop_2d(coeff_2d,sorted_values)
0002 
0003 % This function is the main loop of the Iterative Amplitude Adapted Fourier
0004 % Transform  (IAAFT) method to make surrogate fields. The method was developped by
0005 % Schreiber and Schmitz, see e.g. Phys. Rev Lett. 77, pp. 635-, 1996, for
0006 % statistical non-linearity tests for time series.
0007 % This method makes fields that have a specified amplitude distribution and
0008 % power spectral coefficients. This can be seen as a field made with linear
0009 % dynamics as seen through a static non-linear filter. It works by
0010 % iteratively addaption the amplitude distribution and the Fourier
0011 % coefficients (the phases are not changed in this step).
0012 
0013 % This Matlab version was written by Victor Venema,
0014 % Victor.Venema@uni-bonn.de, http:\\www.meteo.uni-bonn.de\victor, for the
0015 % generation of surrogate cloud fields. First version: May 2003.
0016 
0017 % INPUT:
0018 % coeff_2d:      The 2 dimensional Fourier coefficients that describe the
0019 %                structure and implicitely pass the size of the matrix
0020 % sorted_values: A vector with all the wanted amplitudes (e.g. LWC of LWP
0021 %                values) sorted in acending order.
0022 % OUTPUT:
0023 % y:               The 2D IAAFT surrogate field
0024 % error_amplitude: The amount of addaption that was made in the last
0025 %                  amplitude addaption relative to the total standard deviation.
0026 % error_spec:      The amont of addaption that was made in the last fourier
0027 %                  coefficient addaption relative to the total standard deviation
0028 
0029 % settings
0030 errorThresshold = 0.005   % Maximum error level that is allowed for both error measures; is probably too high for non-linearity testing.
0031 timeThresshold  = 6000;   % Time in seconds that this loop maximally uses in cpu time. Meant as protection against an infinite loop.
0032 speedThresshold = 1e-5;   % Minimal convergence speed in the maximum error.
0033 
0034 % Initialse function
0035 [no_values_y, no_values_x] = size(coeff_2d);
0036 oldTotalError = 1e6;
0037 speed = 1;
0038 error_amplitude=1;
0039 error_spec=1;
0040 standard_deviation=std(sorted_values);
0041 
0042 % The method starts with a randomized uncorrelated time series y with the pdf of
0043 % sorted_values
0044 
0045 if 0
0046    load 2dtest.mat 
0047    x1=setdiff(1:no_values_y^2,xi); 
0048    x0=find(sorted_values<=0.2147-1e-3);
0049    [dummy,index]=sort(rand(size(x0)));
0050    surrogate(xi)=0.2147;
0051    surrogate(x1(index)) = sorted_values(x0);
0052 end
0053 
0054 if 1
0055    [dummy,index]=sort(rand(size(sorted_values)));
0056    surrogate(index) = sorted_values;
0057 end
0058 
0059 surrogate=reshape(surrogate,no_values_y,no_values_x);
0060 
0061 
0062 % Main intative loop
0063 
0064 tic;
0065 while ( (error_amplitude > errorThresshold | error_spec > errorThresshold) & (toc < timeThresshold) & (speed > speedThresshold) ) %0.0001 300
0066     % addapt the power spectrum
0067     old_surrogate = real(surrogate);    
0068     x=ifft2(surrogate);
0069     phase = angle(x);
0070     x = coeff_2d .* exp(i*phase);
0071     surrogate = fft2(x);
0072     difference=mean(mean(abs(real(surrogate)-old_surrogate)));
0073     error_spec = difference/standard_deviation
0074     % addept the amplitude distribution
0075     old_surrogate = real(surrogate);
0076     surrogate=reshape(surrogate,1,no_values_y*no_values_x);
0077     
0078     if 1 
0079        [dummy,index]=sort(real(surrogate));
0080        surrogate(index)=sorted_values;
0081     end     
0082 
0083     if 0
0084        [dummy,index]=sort(real(surrogate(x1)));
0085        surrogate(x1(index))=sorted_values(x0);
0086        surrogate(xi)=0.2147;
0087     end
0088 
0089     surrogate=reshape(surrogate,no_values_y,no_values_x);
0090     
0091     difference=mean(mean(abs(real(surrogate)-old_surrogate)));
0092     error_amplitude = difference/standard_deviation
0093     
0094 
0095     totalError = error_spec + error_amplitude
0096     speed = (oldTotalError - totalError) / totalError;
0097     oldTotalError = totalError;
0098 end
0099 
0100 surrogate = real(surrogate);
0101 figure  
0102 imagesc(real(surrogate)) 
0103 
0104 %keyboard
0105

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005