Prima versione ufficiale powdist_f.m, ripulita e sistemata

F0 = F1(1:floor(N/2)+1);

nblog = 100;

Flog = logspace(log10(F1(2)/2),log10(F0(end)+0.5/Tseg),nblog);

Fqlog = zeros(size(Flog),'like',Flog);

for i = 1:length(Flog)-1

    [~,minind] = min(abs(F0(:)-0.5/Tseg-(Flog(i))));

    Fqlog(i) = F0(minind)-0.5/Tseg;

Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);

figure('Units','centimeters','Position', [0.75 0.75 21 29.7])

t = tiledlayout(3,1,'Padding','tight');

% t.Units = 'centimeters';

% t.OuterPosition = [0.75 0.75 29.7 21];

figure('Units','centimeters','Position', [0.75 0.75 19 20.5])

% set(gcf, 'DefaultAxesFontName', 'Futura Book');

t = tiledlayout(3,1);

nexttile

semilogx(F0,Y0)

loglog(F0,Y0)

title('Power spectrum (PS) up to Nyq. freq.')

xlabel('Frequency [Hz]')

ylabel('Leahy-normalized power (no corr.)')

nexttile

yme = mean(Y0(round(Hfreq*Tseg)+1:end));

semilogx(Fqlmid,Yqlbin)

loglog(Fqlmid,Yqlbin)

xlim([10 NyqFreq+1000])

ylim([yme*0.8 yme*1.2])

yline(yme,'--',['Avg. pow. for f > ',num2str(fix(Hfreq)),' Hz'])

maxcount = max(xtemp);

datadist = zeros(1,maxcount);

for bu = 1:maxcount

    datadist(bu) = length(xtemp(xtemp == bu))/N;

    datadist(bu) = length(xtemp(xtemp == bu));

end

ptot0 = length(xtemp(xtemp == 0))/N;

mi = -log(ptot0);

% mivalues(i,m) = mi;

extrpoiss = poisspdf([1:maxcount],mi);

extrpoiss = poisspdf(1:maxcount,mi);

epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

meanpow = mean(Y0(round(Hfreq*Tseg)+1:end))/2;

vsuedat = var(xtemp)/mean(xtemp);

% vsuedat = var(xtemp)/mean(xtemp);

ncr = 8;

% epsvalues(i,m) = epscr;

% epsmedio = epsmedio + epscr/M;

pcr = 1- (1-epscr)^(1/ncr);

% expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

qcr = (1-epscr)^(1/ncr);

p8_1 = zeros(1,maxcount);

p8_1(1) = (1-epscr);

end

expec1 = sum((1:maxcount).*p8_1(1:maxcount));

var1 = epscr + (epscr^2)*(7-9/ncr)/2;

% enf = 1 + var1/(expec1^2);

expectot = expec1*mi;

vartot = var1*mi + mi*expec1^2;

vsuetot(1) = vartot/expectot;

% varvalues(i,m,1) = meanpow/vsuetot;

% meanpow/vsuedat

extrtot8 = ptotk(p8_1,extrpoiss);

extrtot8 = N.*ptotk(p8_1,extrpoiss);

ncr = 4;

pcr = 1- (1-epscr)^(1/ncr);

% expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

var1 = epscr + (epscr^2)*(7-9/ncr)/2;

% enf = 1 + var1/(expec1^2);

qcr = (1-epscr)^(1/ncr);

p4_1 = zeros(1,maxcount);

p4_1(1) = (1-epscr);

expectot = expec1*mi;

vartot = var1*mi + mi*expec1^2;

vsuetot(2) = vartot/expectot;

% varvalues(i,m,2) = meanpow/vsuetot;

% meanpow/vsuedat

extrtot4 = ptotk(p4_1,extrpoiss);

extrtot4 = N.*ptotk(p4_1,extrpoiss);

% The case for n = 1 is especially simple: it's a Polya–Aeppli dist.;

% [see Univ. Disc. Dist. - Johnson, Kemp, and Kotz, sec 9.7, same

% notation, while in their chap. 5 their p is our q and vice versa]

ncr = 1;

% ncr = 1; % unneeded

pcr = epscr;

% enf = 1 + var1/(expec1^2);

qcr = (1-epscr);

p1_1 = zeros(1,maxcount);

p1_1(1:maxcount) = (pcr.^(0:maxcount-1)).*(1-pcr);

expectot = expec1*mi;

vartot = var1*mi + mi*expec1^2;

vsuetot(3) = vartot/expectot;

% varvalues(i,m,3) = meanpow/vsuetot;

% meanpow/vsuedat

extrtot1 = ptotk(p1_1,extrpoiss);

extrtot1 = N.*ptotk(p1_1,extrpoiss);

blindpoiss = poisspdf((1:maxcount),mean(xtemp));

blindpoiss = N.*poisspdf((1:maxcount),mean(xtemp));

t = tiledlayout(2,1);

nexttile

semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],blindpoiss,[1:maxcount],datadist,'+')

legend('1-pixel crosstalk','4-pixel crosstalk','8-pixel crosstalk','Pure Poissonian','Data distr.')

semilogy(1:maxcount,extrtot1,1:maxcount,extrtot4,1:maxcount,extrtot8,1:maxcount,blindpoiss,1:maxcount,datadist,'+')

legend({'1-pixel crosstalk','4-pixel crosstalk','8-pixel crosstalk','Pure Poissonian','Data distribution'},'Location','southwest')

xlabel('Counts')

ylabel('Number of bins with x counts')

title('Extrapolated gener. Poiss. dist. vs data')

babbo = nexttile;

semilogy([1:maxcount],abs(extrtot1-datadist)./datadist,'x',[1:maxcount],abs(extrtot4-datadist)./datadist,'o',[1:maxcount],abs(extrtot8-datadist)./datadist,'*',[1:maxcount],abs(blindpoiss-datadist)./datadist,'square')

semilogy(1:maxcount,abs(extrtot1-datadist)./datadist,'x',1:maxcount,abs(extrtot4-datadist)./datadist,'o',1:maxcount,abs(extrtot8-datadist)./datadist,'*',1:maxcount,abs(blindpoiss-datadist)./datadist,'square')

ylim([0.0001 2])

% legend('abs(extrtot1-datadist)/datadist','abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

yline([0.01 0.1 0.5 1.0],'--')