Generates surrogate data for real and complex signals Surrogate data generation using iterated amplitude adjusted fourier transform (iAAFT) method for real-valued series, and complex iAAFT method for complex series. USAGE: Xs = surrogate(X, Ns) INPUTS: X: Input time series (can be real-valued or complex) Ns: No. of surrogates to be generated OUTPUTS: Xs: Generated surrogates A Delay Vector Variance (DVV) toolbox for MATLAB (c) Copyright Danilo P. Mandic 2008 http://www.commsp.ee.ic.ac.uk/~mandic/dvv.htm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You can obtain a copy of the GNU General Public License from http://www.gnu.org/copyleft/gpl.html or by writing to Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function Xs = surrogate(X, Ns) 0002 % Generates surrogate data for real and complex signals 0003 % 0004 % Surrogate data generation using iterated amplitude adjusted fourier 0005 % transform (iAAFT) method for real-valued series, and complex iAAFT method 0006 % for complex series. 0007 % 0008 % 0009 % USAGE: Xs = surrogate(X, Ns) 0010 % 0011 % INPUTS: 0012 % X: Input time series (can be real-valued or complex) 0013 % Ns: No. of surrogates to be generated 0014 % 0015 % OUTPUTS: 0016 % Xs: Generated surrogates 0017 % 0018 % 0019 % A Delay Vector Variance (DVV) toolbox for MATLAB 0020 % (c) Copyright Danilo P. Mandic 2008 0021 % http://www.commsp.ee.ic.ac.uk/~mandic/dvv.htm 0022 % 0023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0024 % This program is free software; you can redistribute it and/or modify 0025 % it under the terms of the GNU General Public License as published by 0026 % the Free Software Foundation; either version 2 of the License, or 0027 % (at your option) any later version. 0028 % 0029 % This program is distributed in the hope that it will be useful, 0030 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0031 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0032 % GNU General Public License for more details. 0033 % 0034 % You can obtain a copy of the GNU General Public License from 0035 % http://www.gnu.org/copyleft/gpl.html or by writing to 0036 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0038 0039 % Default Parameters 0040 if (nargin<1) 0041 error('Not enough Input arguments'); 0042 end 0043 if (nargin<2) 0044 Ns = 25; 0045 end 0046 0047 % Initial conditions and parameter initializations 0048 iter = 0; 0049 max_it = 1000; % Maximum iterations 0050 error_threshold = 1e-5; 0051 MSE_start = 100; 0052 MSE = 1000; 0053 0054 % Makes input vector X a column vector 0055 if (size(X,2) > size(X,1)) 0056 X = X'; 0057 end 0058 0059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0060 % Implements univariate iAAFT algorithm for real-valued signals 0061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0062 if (isreal(X)) 0063 0064 Xs = zeros(length(X),Ns); 0065 0066 % Stores the Amplitude vector of the original time series 0067 X_amp = abs(fft(X)); 0068 0069 % Stores the sorted version(ascending) of the original time series 0070 X_sorted = sort(X); 0071 0072 % Iteration process for multiple surrogate time series generation 0073 for a = 1:Ns 0074 0075 % Random permutation of the original time series 0076 temp = randperm(length(X)); 0077 X_random = X(temp); 0078 0079 % Initializations 0080 r_prev = X_random; 0081 MSE_prev = MSE_start; 0082 0083 % Iterate untill convergence condition is met or max iterations are reached 0084 while (abs(MSE-MSE_prev) > error_threshold && iter < max_it) 0085 0086 MSE_prev = MSE; 0087 0088 % Amplitude spectrum matching 0089 ang_r_prev = angle(fft(r_prev)); 0090 s = ifft(X_amp .* exp(ang_r_prev.*sqrt(-1))); 0091 0092 % Rank ordering in order to scale to original signal distribution 0093 [s_sort, Ind] = sort(s); 0094 r(Ind,:) = X_sorted; 0095 0096 % Convergence Metric calculation 0097 MSE = mean(abs(X_amp - abs(fft(r)))); 0098 0099 r_prev = r; 0100 iter = iter+1; 0101 end 0102 0103 % iter 0104 Xs(:,a) = r; 0105 iter = 0; 0106 end 0107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0108 % Implements Complex iAAFT Algorithm for complex signals 0109 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0110 else 0111 0112 % Initial conditions 0113 Xs = zeros(length(X),Ns); 0114 r_real = zeros(length(X),1); 0115 r_imag = zeros(length(X),1); 0116 r = zeros(length(X),1); 0117 r_prev = zeros(length(X),1); 0118 0119 % Stores the Amplitude vector of the original time series 0120 X_amp = abs(fft(X)); 0121 0122 % Stores the sorted version(modulus, ascending), of the original time series 0123 C_sorted = sort(abs(X)); 0124 0125 % Sorted versions(real and imaginary), of original time series 0126 C_real_sorted = sort(real(X)); 0127 C_imag_sorted = sort(imag(X)); 0128 0129 % Iteration process for multiple surrogate time series generation 0130 for a = 1:Ns 0131 0132 % Random permutation of the original time series 0133 r_prev = X(randperm(length(X))); 0134 MSE_prev = MSE_start; 0135 0136 % Iterate untill convergence condition is met or max iterations are reached 0137 while( abs(MSE-MSE_prev) > error_threshold && iter < max_it) 0138 0139 MSE_prev = MSE; 0140 0141 % Amplitude spectrum matching 0142 ang_r_prev = angle(fft(r_prev)); 0143 s = ifft(X_amp .* exp(ang_r_prev.*sqrt(-1))); 0144 0145 % Rank ordering real and imaginary parts of s to match that of original signal 0146 [temp , Ind_real] = sort(real(s)); 0147 [temp , Ind_imag] = sort(imag(s)); 0148 r_real(Ind_real) = C_real_sorted; 0149 r_imag(Ind_imag) = C_imag_sorted; 0150 0151 r = complex(r_real, r_imag); 0152 0153 % Rank ordering in order to scale to original signal distribution 0154 [temp, Ind_abs] = sort(abs(r)); 0155 AUX1 = abs(temp); 0156 r(Ind_abs) = r(Ind_abs).*(abs(C_sorted)./(AUX1+(AUX1==0))); 0157 0158 % Convergence Metric 0159 MSE = mean(abs(X_amp - abs(fft(r)))); 0160 r_prev = r; 0161 iter = iter + 1; 0162 end 0163 0164 % iter 0165 Xs(:,a) = r; 0166 iter = 0; 0167 end 0168 0169 end