+
−% A test script for ao/whiten2D
+
−% 
+
−% L. Ferraioli 18-11-08
+
−% 
+
−% $Id: test_ao_whiten2D.m,v 1.6 2009/04/20 14:31:14 luigi Exp $
+
−% 
+
−
+
−%%
+
−
+
−clear all
+
−
+
−%% CSD Noise models
+
−
+
−% Mdc2 noise model
+
−% mod = ltpnoise();
+
−
+
−load mpsd.mat % load mpsd.mat first column is f then csd11, csd12 and csd22
+
−
+
−f = mpsd(:,1); % frequency vector
+
−csd11 = mpsd(:,2);
+
−csd12 = mpsd(:,3);
+
−csd22 = mpsd(:,4);
+
−fs = 10; % sampling frequency
+
−
+
−% Building AOs
+
−ltp_csd11 = ao(fsdata(f,csd11,fs));
+
−ltp_csd12 = ao(fsdata(f,csd12,fs));
+
−ltp_csd21 = ao(fsdata(f,conj(csd12),fs));
+
−ltp_csd22 = ao(fsdata(f,csd22,fs));
+
−
+
−%% 
+
−
+
−iplot(ltp_csd11,ltp_csd12,ltp_csd22)
+
−
+
−%% Colored Noise generation
+
−
+
−% Make white noise
+
−
+
−a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 10000));
+
−a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 10000));
+
−a = [a1 a2];
+
−
+
−% Generation from white noise
+
−
+
−pl = plist(...
+
−    'csd11', ltp_csd11,...
+
−    'csd12', ltp_csd12,...
+
−    'csd21', ltp_csd21,...
+
−    'csd22', ltp_csd22,...
+
−    'MaxIter', 70, ...
+
−    'PoleType', 2, ...
+
−    'MinOrder', 20, ...
+
−    'MaxOrder', 30, ...
+
−    'Weights', 2, ...
+
−    'Plot', false,...
+
−    'MSEVARTOL', 1e-2,...
+
−    'FITTOL', 1e-3,...
+
−    'UseSym', 0,...
+
−    'Disp', false);
+
−
+
−ac = noisegen2D(a, pl);
+
−
+
−%% Checking generated noise
+
−
+
−acxx1 = ac(1).psd;
+
−acxx2 = ac(2).psd;
+
−acxx12 = ac.cpsd;
+
−
+
−iplot(abs(acxx1),abs(ltp_csd11))
+
−iplot(abs(acxx12),abs(ltp_csd12))
+
−iplot(abs(acxx2),abs(ltp_csd22))
+
−
+
−%% Noise whitening
+
−
+
−pl = plist(...
+
−    'csd11', ltp_csd11, ...
+
−    'csd12', ltp_csd12, ...
+
−    'csd21', ltp_csd21, ...
+
−    'csd22', ltp_csd22, ...
+
−    'MaxIter', 75, ...
+
−    'PoleType', 3, ...
+
−    'MinOrder', 20, ...
+
−    'MaxOrder', 40, ...
+
−    'Weights', 2, ...
+
−    'Plot', false,...
+
−    'MSEVARTOL', 1e-2,...
+
−    'FITTOL', 1e-3,...
+
−    'UseSym', 0,...
+
−    'Disp', false);
+
−
+
−aw = whiten2D(ac, pl);
+
−
+
−%% Checking whitened noise
+
−
+
−awxx1 = aw(1).psd;
+
−awxx2 = aw(2).psd;
+
−
+
−%%
+
−iplot(awxx1,awxx2)
+
−
+
−%% Checking
+
−
+
−aw1 = aw(1);
+
−aw2 = aw(2);
+
−awxx1 = awxx(1,1);
+
−awxx2 = awxx(2,2);
+
−
+
−mw1 = aw1.mean; % 0 expected
+
−vw1 = aw1.var; % 1 expected
+
−
+
−mw2 = aw2.mean; % 0 expected
+
−vw2 = aw2.var; % 1 expected
+
−
+
−mwxx1 = awxx1.mean; % 0.2 expected
+
−vwxx1 = awxx1.var;
+
−
+
−mwxx2 = awxx2.mean; % 0.2 expected
+
−vwxx2 = awxx2.var;
+
−
+
−% correlation coefficients
+
−% diagonal terms are aw1 and aw2 variances
+
−% cross terms are data series covariances, a value close to 0 means no
+
−% correlation between the data series
+
−crc = corrcoef([aw1.data.y.' aw2.data.y.']);
+
−
+
−
+
−% END