view m-toolbox/test/test_ao_tfe_variance_montecarlo.m @ 17:7afc99ec5f04
database-connection-manager
Update ao_model_retrieve_in_timespan
author |
Daniele Nicolodi <nicolodi@science.unitn.it> |
date |
Mon, 05 Dec 2011 16:20:06 +0100 (2011-12-05) |
parents |
f0afece42f48 |
children |
|
line source
% test_ao_tfe_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_tfe_variance_montecarlo.m,v 1.1 2009/08/11 14:20:10 miquel Exp $
% function test_ao_tfe_variance_montecarlo()
clear
% data
nsecs = 200;
fs = 10;
pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'sin(2*pi*2*t) + randn(size(t))');
% Make a filter
f1 = miir(plist('type', 'highpass', 'fc', 4, 'fs', fs));
% Window
Nfft = 100;
win = specwin('Hanning', Nfft);
pl2 = plist('Nfft',Nfft, 'win',win,'Olap',-1)
% loop
for i = 1:100
a1 = ao(pl);
a2 = filter(a1,plist('filter', f1));
t1(i) = tfe(a1,a2,pl2);
% matlab's
[txy, f] = tfestimate(a1.data.y, a2.data.y, win.win, Nfft/2, Nfft, a1.data.fs);
t2(i) = ao(fsdata(f.', txy.'));
end
%% mean
index = 6;
% compare mean
mn = [mean(t1(:).y(index)) mean(t2(:).y(index))]
% error
err = std(t1(:).y(index))
% compare standard deviation
clear rel
for i =1:len(t1(1))
mn(i) = [mean(t1(:).y(i))]; % both means are equal
rel(:,i) = [std(t1(:).y(i)) mean(t1(:).dy(i))]/abs(mn(i));
end
figure
loglog(t1(1).x,rel')
figure
loglog(t1(1).x,100*abs(rel(2,:)-rel(1,:)))
ylabel('difference (%)')