Home > atmlab > randomize > iaaft > iaaft_loop_3d.m

iaaft_loop_3d

PURPOSE ^

INPUT:

SYNOPSIS ^

function [surrogate, error_amplitude, error_spec] = iaaft_loop_3d(coeff_3d,sorted_values_prof)

DESCRIPTION ^

 INPUT: 
 coeff_3d:           The 3 dimensional Fourier coefficients that describe the
                     structure and implicitely pass the size of the matrix
 sorted_values_prof: A vector with the wanted amplitudes (e.g. LWC of LWP
                     values) sorted in acending order for each height level.
 OUTPUT:
 y:                  The 3D IAAFT surrogate field
 error_amplitude:    The amount of addaption that was made in the last
                     amplitude addaption relative to the total standard deviation.
 error_spec:         The amont of addaption that was made in the last fourier
                     coefficient addaption relative to the total standard deviation

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

iaaft_loop_3d.m

SOURCE CODE ^

0001 function [surrogate, error_amplitude, error_spec] = iaaft_loop_3d(coeff_3d,sorted_values_prof)
0002 % INPUT:
0003 % coeff_3d:           The 3 dimensional Fourier coefficients that describe the
0004 %                     structure and implicitely pass the size of the matrix
0005 % sorted_values_prof: A vector with the wanted amplitudes (e.g. LWC of LWP
0006 %                     values) sorted in acending order for each height level.
0007 % OUTPUT:
0008 % y:                  The 3D IAAFT surrogate field
0009 % error_amplitude:    The amount of addaption that was made in the last
0010 %                     amplitude addaption relative to the total standard deviation.
0011 % error_spec:         The amont of addaption that was made in the last fourier
0012 %                     coefficient addaption relative to the total standard deviation
0013 
0014 % settings
0015 errorThresshold = 1e-4;    % The addaptation made in the last iterative step relative to the total standard deviation
0016 timeThresshold  = inf;     %14*3600; % CPU-time in seconds that the iteration maximally uses
0017 speedThresshold = 1e-4;    % Minimal convergence speed in the maximum error.
0018 
0019 % Initialise function
0020 [no_values_y, no_values_x, no_values_z] = size(coeff_3d);
0021 error_amplitude = 1;
0022 error_spec      = 1;
0023 oldTotalError   = 10;
0024 speed = 1;
0025 standard_deviation = std(sorted_values_prof(:));
0026 
0027 % The method starts with a randomized uncorrelated time series y with the pdf of
0028 % sorted_values
0029 surrogate = zeros(no_values_y*no_values_x,no_values_z);
0030 for j=1:no_values_z
0031     [dummy,index] = sort(rand(size(sorted_values_prof(:,j))));
0032     surrogate(index,j) = sorted_values_prof(:,j);
0033 end
0034 surrogate=reshape(surrogate,no_values_y,no_values_x,no_values_z);
0035 
0036 % Main intative loop
0037 counter = 0;
0038 at=tic;
0039 while ( (error_amplitude > errorThresshold | error_spec > errorThresshold) & (toc(at) < timeThresshold) & (speed > speedThresshold) ) %0.0001 300
0040     % addapt the power spectrum
0041     old_surrogate = surrogate;    
0042     x=ifftn(double(surrogate));
0043     phase = angle(x);
0044 
0045     x = double(coeff_3d) .* exp(i*phase);
0046     surrogate = fftn(double(x));
0047     difference = mean(mean(mean(abs(double(real(surrogate))-double(real(old_surrogate))))));
0048     error_spec = difference/standard_deviation;
0049     
0050     % addept the amplitude distribution
0051     old_surrogate = surrogate;
0052     surrogate = reshape(surrogate, no_values_y*no_values_x, no_values_z);
0053     for j=1:no_values_z
0054         [dummy, index] = sort(real(surrogate(:,j)));
0055         surrogate(index,j) = sorted_values_prof(:,j);
0056     end
0057     surrogate = reshape(surrogate, no_values_y,no_values_x, no_values_z);
0058     difference = mean(mean(mean(abs(double(real(surrogate))-double(real(old_surrogate))))));
0059     error_amplitude = difference/standard_deviation;
0060     counter = counter + 1;
0061 
0062     TotalError = error_spec + error_amplitude;
0063     speed = (oldTotalError - TotalError) / TotalError;
0064     oldTotalError = TotalError;
0065 end
0066 
0067 error_spec;
0068 error_amplitude;
0069 surrogate = real(surrogate);

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