%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %
%%           synthArVH.m               %
%%                                     %
%%   P. Chainais & R. Riedi & P. Abry  %
%%                                     %
%%             march 2004              %
%%                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Synthesizes infinitely divisible cascading motion A_r(t),
%                                        and random walk VH(t),
% associated to a multifractal (power law scaling) Compound Poisson Cascade 
% or Log Normal cascade. 
% 
% See references : ( www.isima.fr/~chainais/PUB/publications.html ) 
%
% On non scale invariant infinitely divisible cascades 
% LIMOS/RR-04-06 - Auteurs : P. Chainais1, R. Riedi, P. Abry . 
% at    http://www.isima.fr/limos/publications.htm 
%
% On non scale invariant infinitely divisible cascades ,
% P. Chainais, R. Riedi, et P. Abry,Ê in Actes du GRETSI, Paris, 2003. 
% 
% Scale invariant infinitely divisible  cascades ,
% P. Chainais, R. Riedi, et P. Abry,Ê in Proceedings of  PSIP'03, Grenoble, 2003. 
% 
% Compound  Poisson Cascades ,
% P. Chainais, R. Riedi, et P. Abry,Ê Annales de l'UniversitŽ Blaise Pascal 
% Workshop Self-Similarity and Applications, Clermont-Ferrand 2002. 
%
%%%%%%%%%%%%%%%%%
%
%  Uses   dilatrandn.m
%         synthfbmcircul.m
%         
% ----------------------------        
% Input : 
%       . H   : Hurst parameter of underlying fractional Brownian motion,
%       . Qr : IDC noise,
%       . Dt : time sampling rate (ex Dt = 2.rmin).
%       . JL    : number of points of the signal
%       . JS    : rmin=2^(-JS),
%       . oversamprate : over sampling rate for injection of Ar in fBm to create VH,
%       . printout : 1 if plot needed.
%       
% Output :
%       . Ar : IDC motion,
%       . VH : IDC-fractional Brownian motion,
%       . BH : corresponding realization of fBm,
%        
%%%%%%%%%%%%%%%%%%%
% Example : 
%  
% H=1/3 ; oversamprate = 16 ; printout = 0;
% [Ar,VH,BH] = synthArVH(H,Qr,Dt,oversamprate,printout) ;
%
%%%%%%%%%%%%%%%%%%%%%%%%%

function [Ar,VH,BH] = synthArVH(H,Qr,Dt,oversamprate,printout) ;

if nargin>2
    
%--- Provides Ar
Ar = cumsum(Qr)*Dt ;

if nargout>1  % allows to get Ar only if needed
    
%--- provides fBm on [0, VHduration]
VHduration = Ar(end) ;
nfbm = 1+2^ceil(log2(VHduration/Dt*oversamprate)) ;
if (H-1/2)
    [BH] = synthfbmcircul(nfbm,H) ;   % fbmcircul.m plus rapide pour N=2^k+1
else
    [BH] = cumsum(randn(1,nfbm)) ;
end

BH = VHduration^H*BH ;       % BH was synthesized on [0 1] -> [0 VHduration] 
% [BH] = fbmcircul(2^(ceil(log2(nfbm)))+1,H) ; 

%--- Creates IDCRW
ind = round(Ar/Dt*oversamprate) ;
%ind=ind(find(ind));
ind=max(1,ind);
VH = BH(ind) ;
%nVH = length(VH) ;
%BH = BH(1:oversamprate:nVH*oversamprate) ;
BH = BH(1:oversamprate:end) ;

%--- Affiche
if printout == 1
	%time = 1:1:nVH ;
    time = (0:(length(VH)-1))*Dt ;
	figure
	plot(linspace(0,VHduration,length(BH)),BH,'r',time,VH,'b')
	legend('FBM','IDCRW')
	grid
end

end % if nargout

end % if nargout