%--------------------------------------------------------------------------
% DWT-QIC estimator
%--------------------------------------------------------------------------
% Copyright (C) 2008  Giada Giorgi, University of Padova, Italy.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
%USAGE
%  function[temp]=dwtquantile(x,J)
%INPUT
%  x = time series samples
%  J = maximum time scale
%OUTPUT
%  temp.ca = approximation coefficients
%  temp.cd = detail coefficients 
%  temp.va = variances of approx. coeff. at different time scales
%  temp.vd = variances of detail coeff. at different time scales
%  temp.qa = quantiles of approx. coeff. at different time scales
%  temp.qd = quantiles of detail coeff. at different time scales
%AUTHOR
%-------------------------------------------------------------------------


function[temp]=dwtquantile(x,J)
p=[0.1 0.5 1 5 10 20 30 40 50 60 70 80 90 95 99 99.5 99.9];
Np=length(p);
j=1;
[cA,cD]=dwt(x,'haar');
va(1)=var(cA);  vd(1)=var(cD);
qa(1,1:Np)=prctile(cA,p);
qd(1,1:Np)=prctile(cD,p);
ca{1}=cA;       cd{1}=cD;
clear cA; clear cD;
for j=2:J
    [cA,cD]=dwt(ca{j-1},'haar');
    va(j)=var(cA);  vd(j)=var(cD);
    qa(j,1:Np)=prctile(cA,p);
    qd(j,1:Np)=prctile(cD,p);
    ca{j}=cA; cd{j}=cD;
    clear cA; clear cD;
end
temp.ca=ca;
temp.cd=cd;
temp.va=va;
temp.vd=vd;
temp.qa=qa;
temp.qd=qd;

figure(21); hold;
plot((1:J),log2(temp.qd(:,17)),'-kp','MarkerSize',6,'MarkerFacecolor','k'); %99.9
plot((1:J),log2(temp.qd(:,16)),'-k<','MarkerSize',6,'MarkerFacecolor','w'); %99.5
plot((1:J),log2(temp.qd(:,15)),'-k^','MarkerSize',6,'MarkerFacecolor','k'); %99
plot((1:J),log2(temp.qd(:,14)),'-kv','MarkerSize',6,'MarkerFacecolor','w'); %95
plot((1:J),log2(temp.qd(:,13)),'-kd','MarkerSize',6,'MarkerFacecolor','k'); %90
plot((1:J),log2(temp.qd(:,12)),'-k*','MarkerSize',6); %80
plot((1:J),log2(temp.qd(:,11)),'-ko','MarkerSize',6,'MarkerFacecolor','k'); %70
plot((1:J),log2(temp.qd(:,10)),'-ks','MarkerSize',6,'MarkerFacecolor','w'); %60
axis([0 J -inf inf]);
title('\fontsize{16}wavelet spectrum: quantiles');
xlabel('\fontsize{16}scales j');
ylabel('\fontsize{16}log_2[quantile]');
set(gca,'FontSize',14);
legend('\fontsize{12}99.9%', '\fontsize{12}99.5%', '\fontsize{12}99%',...
        '\fontsize{12}95%', '\fontsize{12}90%', '\fontsize{12}80%',...
        '\fontsize{12}70%', '\fontsize{12}60%','Location','NorthWest');
grid;
fprintf(1,'-------------------------------------------------------\n');
reply=input('Estimate H? [YES]=1 [NO]=another\n');
while(reply==1)
    nq=input('Choose the quantile:\n[1]60% [2]70% [3]80% [4]90% [5]95% [6]99% [7]99.5% [8]99.9%\n');
    j1=input('Insert j1:\n');
    j2=input('Insert j2:\n');
    coeff=polyfit((j1:1:j2),log2(temp.qd(j1:j2,nq+9)'),1);
    figure(21); plot((j1:1:j2),coeff(1)*(j1:1:j2)+coeff(2),'--r','LineWidth',2);
    H=coeff(1)+ 0.5;
    fprintf(1,'H=%d\n',H);
    fprintf(1,'-------------------------------------------------------\n');
    reply=input('Estimate H? [YES]=1 [NO]=another\n');
end

%Quantili coefficienti di approssimazione
figure(22); hold;
plot((1:J),log2(temp.qa(:,17)-temp.qa(:,9)),'-kp','MarkerSize',6,'MarkerFacecolor','k'); %99.9
plot((1:J),log2(temp.qa(:,16)-temp.qa(:,9)),'-k<','MarkerSize',6,'MarkerFacecolor','w'); %99.5
plot((1:J),log2(temp.qa(:,15)-temp.qa(:,9)),'-k^','MarkerSize',6,'MarkerFacecolor','k'); %99
plot((1:J),log2(temp.qa(:,14)-temp.qa(:,9)),'-kv','MarkerSize',6,'MarkerFacecolor','w'); %95
plot((1:J),log2(temp.qa(:,13)-temp.qa(:,9)),'-kd','MarkerSize',6,'MarkerFacecolor','k'); %90
plot((1:J),log2(temp.qa(:,12)-temp.qa(:,9)),'-k*','MarkerSize',6); %80
plot((1:J),log2(temp.qa(:,11)-temp.qa(:,9)),'-ko','MarkerSize',6,'MarkerFacecolor','k'); %70
plot((1:J),log2(temp.qa(:,10)-temp.qa(:,9)),'-ks','MarkerSize',6,'MarkerFacecolor','w'); %60
axis([0 J -inf inf]);
title('\fontsize{16}wavelet spectrum: quantiles');
xlabel('\fontsize{16}scales j');
ylabel('\fontsize{16}log_2[quantile]');
set(gca,'FontSize',14);
legend('\fontsize{12}99.9%', '\fontsize{12}99.5%', '\fontsize{12}99%',...
        '\fontsize{12}95%', '\fontsize{12}90%', '\fontsize{12}80%',...
        '\fontsize{12}70%', '\fontsize{12}60%','Location','NorthWest');
grid;
fprintf(1,'-------------------------------------------------------\n');
reply=input('Estimate H? [YES]=1 [NO]=another\n');
while(reply==1)
    nq=input('Choose the quantile:\n[1]60% [2]70% [3]80% [4]90% [5]95% [6]99% [7]99.5% [8]99.9%\n');
    j1=input('Insert j1:\n');
    j2=input('Insert j2:\n');
    coeff=polyfit((j1:1:j2),log2(temp.qa(j1:j2,nq+9)'-temp.qa(j1:j2,9)'),1);
    figure(22); plot((j1:1:j2),coeff(1)*(j1:1:j2)+coeff(2),'--r','LineWidth',2);
    H=coeff(1)+ 0.5;
    fprintf(1,'H=%d\n',H);
    fprintf(1,'-------------------------------------------------------\n');
    reply=input('Estimate H? [YES]=1 [NO]=another\n');
end

%Spettro wavelet coefficienti di dettaglio
figure(23); hold;
plot((1:J),log2(temp.vd),'-k*','LineWidth',1,'MarkerSize',6);
title('\fontsize{16}wavelet spectrum: variance');
xlabel('\fontsize{16}scales j');
ylabel('\fontsize{16}log_2[var]');
set(gca,'FontSize',14);
grid;
axis([0 J -inf inf]);
fprintf(1,'-------------------------------------------------------\n');
reply=input('Estimate H? [YES]=1 [NO]=altro\n');
while(reply==1)
    j1=input('Insert j1:\n');
    j2=input('Insert j2:\n');
    coeff=polyfit((j1:1:j2),log2(temp.vd(j1:j2)),1);
    figure(23); plot((j1:1:j2),coeff(1)*(j1:1:j2)+coeff(2),'--r','LineWidth',2);
    H=(coeff(1)+ 1)/2;
    fprintf(1,'H=%d\n',H);
    fprintf(1,'-------------------------------------------------------\n');
    reply=input('Estimate H? [YES]=1 [NO]=altro\n');
end
