view m-toolbox/test/test_ao_psd_variance_montecarlo.m @ 48:16aa66670d74
database-connection-manager
Fix LTPDA Preferences tooltip
author |
Daniele Nicolodi <nicolodi@science.unitn.it> |
date |
Tue, 06 Dec 2011 19:07:27 +0100 (2011-12-06) |
parents |
f0afece42f48 |
children |
|
line source
% test_ao_psd_variance_montecarlo
%
% Tests that the standard deviation returned by ao.dy in one
% frequency bin is equivalent to the matlab's std taking
% considering all realisations
%
% M Nofrarias 22-07-09
%
% $Id: test_ao_psd_variance_montecarlo.m,v 1.2 2009/08/11 14:20:10 miquel Exp $
% function test_ao_psd_variance_montecarlo()
clear
% data
nsecs = 500;
fs = 5;
pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
% Window
Nfft = 100;
win = specwin('Hanning', Nfft);
pl2 = plist('Nfft',Nfft, 'win',win,'Olap',-1,'scale','PSD')
% loop
for i = 1:100
a(i) = ao(pl);
b1(i) = psd(a(i),pl2);
% matlab's
[txy, f] = pwelch(a(i).data.y, win.win, Nfft/2, Nfft, a(i).data.fs);
b2(i) = ao(fsdata(f.', txy.'));
end
%% mean
index = 6;
% compare mean
mn = [mean(b1(:).y(index)) mean(b2(:).y(index))]
% error
err = std(b1(:).y(index))
% compare standard deviation
clear rel
for i =1:len(b1(1))
mn(i) = [mean(b1(:).y(i))]; % both means are equal
rel(:,i) = [std(b1(:).y(i)) mean(b1(:).dy(i))]/abs(mn(i));
end
figure
loglog(b1(1).x,rel')
figure
loglog(b1(1).x,rel(2,:)-rel(1,:))
ylabel('difference (%)')