Home > atmlab > randomize > iaaft > siaaft_1d.m

siaaft_1d

PURPOSE ^

INPUT:

SYNOPSIS ^

function [bestSurrogate, spectralDiffBest] = siaaft_loop_1d_experimental(fourierCoeff, sortedValues, counterThresshold, iterativeCounterAmplitude, iterativeDenomAmplitude, algorithmVersion)

DESCRIPTION ^

 INPUT: 
 fourierCoeff:              The 1 dimensional Fourier coefficients that describe the
                            structure and implicitely pass the size of the matrix.
 sortedValues:              A vector with all the wanted amplitudes (e.g. LWC of LWP
                            values) sorted in acending order.
 counterThresshold:         The number of unsuccesful iterations before
                            the algorithm stops.
 iterativeCounterAmplitude: The fraction of amplitude adaptations is give
                            by iterativeCounterAmplitude/iterativeDenomAmplitude.
                            For some algorithm versions 2/4 can be a bit
                            different that 1/2. Please look in the code.
 iterativeDenomAmplitude:   See variable iterativeCounterAmplitude above.
 algorithmVersion:          There are 4 different algorithm versions. 
                            1: Partially stochastic version (default)
                            2: Deterministic version
                            3: Fully stochastic version (where the exact
                               fraction of amplitude adaptations is
                               performed).
                            4: Fully stochastic version (more efficient,
                               but not exactly the right fraction of
                               amplitude adaptations.
 
 OUTPUT:
 bestSurrogate:    The 1D SIAAFT surrogate time series
 spectralDiffBest: The RMS difference between the specified Fourier
                   spectrum (from the template) and the spectrum of the
                   surrogate.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

siaaft_1d.m

SOURCE CODE ^

0001 function [bestSurrogate, spectralDiffBest] = siaaft_loop_1d_experimental(fourierCoeff, sortedValues, counterThresshold, iterativeCounterAmplitude, iterativeDenomAmplitude, algorithmVersion)
0002 % INPUT:
0003 % fourierCoeff:              The 1 dimensional Fourier coefficients that describe the
0004 %                            structure and implicitely pass the size of the matrix.
0005 % sortedValues:              A vector with all the wanted amplitudes (e.g. LWC of LWP
0006 %                            values) sorted in acending order.
0007 % counterThresshold:         The number of unsuccesful iterations before
0008 %                            the algorithm stops.
0009 % iterativeCounterAmplitude: The fraction of amplitude adaptations is give
0010 %                            by iterativeCounterAmplitude/iterativeDenomAmplitude.
0011 %                            For some algorithm versions 2/4 can be a bit
0012 %                            different that 1/2. Please look in the code.
0013 % iterativeDenomAmplitude:   See variable iterativeCounterAmplitude above.
0014 % algorithmVersion:          There are 4 different algorithm versions.
0015 %                            1: Partially stochastic version (default)
0016 %                            2: Deterministic version
0017 %                            3: Fully stochastic version (where the exact
0018 %                               fraction of amplitude adaptations is
0019 %                               performed).
0020 %                            4: Fully stochastic version (more efficient,
0021 %                               but not exactly the right fraction of
0022 %                               amplitude adaptations.
0023 %
0024 % OUTPUT:
0025 % bestSurrogate:    The 1D SIAAFT surrogate time series
0026 % spectralDiffBest: The RMS difference between the specified Fourier
0027 %                   spectrum (from the template) and the spectrum of the
0028 %                   surrogate.
0029 
0030 % input checking
0031 if (nargin < 3)
0032     counterThresshold = 1e4; % Number of iterative steps with no approvement before the algorithm stops.
0033 end
0034 if (nargin < 5)
0035     iterativeCounterAmplitude = 1;
0036     iterativeDenomAmplitude   = 5;
0037 end
0038 if (nargin < 6)
0039     algorithmVersion = 1;
0040 end
0041 
0042 % Settings
0043 verbose = 0;   % Boolean indicating whether to write to the command window or not.
0044 makePlots = 0; % Best used together during debugging (together with verbose).
0045 
0046 % Initialise function
0047 noValues = length(fourierCoeff);
0048 % amplitudeDiff = 1;
0049 spectralDiff = 1;
0050 spectralDiffBest = 100;
0051 standardDeviation = std(sortedValues);
0052 noAdaptedValues = floor(noValues*iterativeCounterAmplitude/iterativeDenomAmplitude);
0053 
0054 indices = init_regular_indices(iterativeCounterAmplitude, iterativeDenomAmplitude, noValues);
0055 noIndices = size(indices,2);
0056 
0057 % Set limits for plots.
0058 if makePlots
0059     xmin = 1;
0060     xmax = noValues;
0061     ymin = sortedValues(1);
0062     ymax = sortedValues(end);
0063 end
0064 
0065 % The method starts with a randomized uncorrelated surrogate time series
0066 % with the amplitudes (PDF) of sorted_values.
0067 [dummy,index] = sort(rand(size(sortedValues)));
0068 surrogate(index) = sortedValues;
0069 
0070 % Main intative loop
0071 counter = 1; % Counts the number of iterations in every stage.
0072 number = 1;  % The index to the vector to use from the matrix indices (for algorithmVersion 1 and 2).
0073 % counterSpec = 1;
0074 % maxCounter = 1e4;
0075 
0076 for stageCounter=1:2
0077     if (stageCounter == 2)
0078         % In this second stage all amplitudes and coefficients are used
0079         % in the iteration, just like in the original iaaft algorithm.
0080         disp('Second stage')
0081         iterativeCounterAmplitude = iterativeDenomAmplitude;
0082         noAdaptedValues = noValues;        
0083         % use the best surrogate as the new one.
0084         surrogate = bestSurrogate;   % The best surrogate of the first phase is the starting surrogate of the second phase
0085         spectralDiffBest = 100;      % make the second phase use the old surrogate as its best one up to now.
0086         counter = 1;
0087         indices = (1:noValues)';     % use all indices, not just a fraction, in the second stage.
0088         noIndices = 1;
0089         number = 1;
0090         counterThresshold = 100; % In the second round with the IAAFT algorithm, it doesn't help to iterate after convergence has stopped.
0091     end
0092     
0093     bestCounter = 0; % This counter counts the number of iteration, and is set to zero each time a new better fitting surrogate is found.
0094     while ( bestCounter < counterThresshold )
0095         % addapt the power spectrum
0096         x=ifft(surrogate);
0097         phase = angle(x);
0098         difference = sqrt( mean( (abs(x)-abs(fourierCoeff)).^2  ) );
0099         spectralDiff = difference/standardDeviation;
0100 %         spectralDiffg(counterSpec) = spectralDiff;
0101 %         counterSpec = counterSpec + 1;
0102 
0103         % If the calculation of the power spectrum reveals that the
0104         % surrogate time series from the previous iteration with perfect
0105         % amplitudes had the best power spectrum up to now, save this
0106         % surrogate and its accuracy and reset bestCounter.
0107         if ( (spectralDiff < spectralDiffBest) & (counter > 1) )
0108             fprintf(1, 'SpectralDiffBest: %e, bestCounter:, %d\n', spectralDiff, bestCounter);
0109             spectralDiffBest  = spectralDiff;
0110 %             amplitudeDiffBest = amplitudeDiff;
0111             bestSurrogate     = surrogate;
0112             bestCounter = 0;
0113 %            phaseCounter;
0114         else
0115             bestCounter = bestCounter + 1;
0116         end
0117            
0118         x = fourierCoeff .* exp(i*phase);
0119         surrogate = fft(x);
0120     
0121         % Plot surrogate after spectral adaptation.
0122         if ( verbose ), spectralDiff, end        
0123         if ( makePlots )
0124             subplot(3,1,1);
0125             plot(real(surrogate))
0126             title('Surrogate after spectal adaptation')
0127             axis([xmin xmax ymin ymax])        
0128             drawnow
0129         end
0130             
0131         % Adapt the amplitude distribution
0132         [dummy, index] = sort(real(surrogate)); % We need only the indices. The first value of index points to the highest value, etc.
0133  
0134         switch algorithmVersion
0135             case 1
0136                 % random selection of noIndices vectors with regularly ordered indices
0137                 number = ceil(rand(1)*noIndices);
0138                 surrogate(index(indices(:,number)))=sortedValues(indices(:,number));                            
0139             case 2
0140                 % deterministic selection of noIndices vectors with regularly ordered indices
0141                 number = increment(number, 1, 1, noIndices);
0142                 surrogate(index(indices(:,number)))=sortedValues(indices(:,number));            
0143             case 3
0144                 % Fully random indices, with exactly the right fraction of amplitude adaptations.
0145                 if (iterativeCounterAmplitude == iterativeDenomAmplitude)
0146                     surrogate(index) = sortedValues;
0147                 else
0148                     [dummy, indices] = sort(rand(size(surrogate)));
0149                     indices = indices(1:noAdaptedValues);
0150                     surrogate(index(indices))=sortedValues(indices);
0151                 end
0152             case 4
0153                 % Fully random indices, more efficient version of algorithmVersion 3.
0154                 if (iterativeCounterAmplitude == iterativeDenomAmplitude)
0155                     surrogate(index) = sortedValues;
0156                 else                
0157                     indices = ceil(rand(noAdaptedValues,1)*noValues);
0158                     surrogate(index(indices))=sortedValues(indices);
0159                 end
0160         end
0161         
0162         % Plot surrogate after amplitude adaptation.
0163         if ( verbose ), amplitudeDiff, end       
0164         if ( makePlots )
0165             subplot(3,1,2);
0166             plot(real(surrogate))
0167             title('Surrogate after amplitude adaptation')
0168             axis([xmin xmax ymin ymax])        
0169             drawnow
0170         end
0171       
0172         counter=counter+1;
0173 %        bestCounter;
0174         drawnow % Allows escaping the while loop with ^c.
0175     end % while not conveged
0176 end % for the 2 stages.
0177 
0178 surrogate = double(surrogate);

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