Mercurial > hg > ltpda
view m-toolbox/test/test_ao_psd_variance_montecarlo.m @ 1:2014ba5b353a database-connection-manager
Remove old code
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Sat, 03 Dec 2011 18:13:55 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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 (%)')