%
%   vibrationdata_wavelet_decomposition_engine_q.m  ver 1.2  by Tom Irvine
%
%
function[acceleration,velocity,displacement,wavelet_table,...
                                     wavelet_low,wavelet_up,alpha,beta]=...
  vibrationdata_wavelet_decomposition_engine_q(t,accel,dt,first,...
                                                   ffmax,nt,nfr,start_time)
%
tp=2*pi;
%
residual=accel;
%
num2=max(size(accel));
%
out1=sprintf(' number of input points= %d ',num2);
disp(out1);
%
duration=t(num2)-t(1);
%
sr=1/dt;
%
out1=sprintf(' sample rate = %10.4g \n',sr);
disp(out1);
%
fl=3/duration;
fu=sr/10;
%
out1=sprintf(' fl = %10.4g Hz   fu = %10.4g Hz \n',fl,fu);
disp(out1);
%
clear y;
%
progressbar;
%
x1r=zeros(nfr,1);
x2r=zeros(nfr,1);
x3r=zeros(nfr,1);
x4r=zeros(nfr,1);
%
alpha=zeros(nfr,1);
beta=zeros(nfr,1);
wavelet_up=zeros(nfr,1);
wavelet_low=zeros(nfr,1);
%
for ie=1:nfr
%
        progressbar(ie/nfr);
%
        fprintf(' frequency case %d \n',ie);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
        % disp(' ref 1');
   
        [ax1r,ax2r,ax3r,ax4r,alphar,betar,wlr,wur]=...
        decomposition_wgen_q_nw(num2,t,residual,duration,fl,fu,nt,ffmax,first,sr,start_time);    
    
    
        % disp(' ref 2');
        x1r(ie)=ax1r;
        x2r(ie)=ax2r;
        x3r(ie)=ax3r;
        x4r(ie)=ax4r;
        alpha(ie)=alphar;
        beta(ie)=betar;
        wavelet_low(ie)=wlr;
        wavelet_up(ie)=wur;        
     
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

        % Pre-calculate constants outside the loop for efficiency
        t1 = x4r(ie) + t(1);
        t2=t1 + tp*x3r(ie)/(2.*x2r(ie));
           
        % Logical indexing for efficient condition check
        mask = t >= t1 & t <= t2;
        
        % Vectorized calculations within the loop
        arg = x2r(ie) .* (t(mask) - t1);
        y = x1r(ie) .* sin(arg ./ double(x3r(ie))) .* sin(arg);
        
        % Update residual using logical indexing
        residual(mask) = residual(mask) - y;

        ave=mean(residual); 
        sd=std(residual);
%
        fprintf(' ave=%12.4g  sd=%12.4g \n\n',ave,sd);
% 
end
%
pause(0.4);
progressbar(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear aaa;
clear vvv;
clear ddd;
clear mm;
%
aaa=zeros(num2,1); 
vvv=zeros(num2,1); 
ddd=zeros(num2,1); 
%
disp(' ');
wavelet_table=zeros(nfr,5);
%
for i=1:nfr
%
        fprintf(' amp=%10.4f   freq=%10.3f Hz   nhs=%d   delay=%10.4f \n',x1r(i),x2r(i)/tp,x3r(i),x4r(i));
%
        wavelet_table(i,1)=i;
        wavelet_table(i,2)=x2r(i)/tp;   
        wavelet_table(i,3)=x1r(i);
        wavelet_table(i,4)=x3r(i); 
        wavelet_table(i,5)=x4r(i);  
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
vlast=zeros(num2,1);
%
iscale=1;
%
progressbar;
for k=1:num2
%
        progressbar(k/num2);
%
        tt=t(k);
%  
        for j=1:nfr
%
            w=0.;
            v=0.;
            d=0.;
%
            t1=x4r(j)+t(1);
%
            if(x2r(j)>=1.0e-10)
%
            else
                x2r(j)=fl*tp;
                x1r(j)=1.0e-20;
                x3r(j)=3;
                x4r(j)=0;
                t1=x4r(j)+t(1);  
            end
            t2=tp*x3r(j)/(2*x2r(j))+t1; 
%
            if( tt>=t1  && tt <=t2  )
%
                arg=x2r(j)*(tt-t1);  
%
                w=  x1r(j)*sin(arg/double(x3r(j)))*sin(arg);
%
                aa=x2r(j)/double(x3r(j));
                bb=x2r(j);
%
                te=tt-t1;
%
                alpha1=aa+bb;
                alpha2=aa-bb;
%
                alpha1te=alpha1*te;
                alpha2te=alpha2*te;   
%
                v1= -sin(alpha1te)/(2*alpha1);
                v2= +sin(alpha2te)/(2*alpha2);
%
                d1= +(cos(alpha1te)-1)/(2*(alpha1^2));
                d2= -(cos(alpha2te)-1)/(2*(alpha2^2));
%
                v=(v2+v1)*iscale*x1r(j);
                d=(d2+d1)*iscale*x1r(j);
%
                vlast(j)=v;
%
            end
%
            aaa(k)=aaa(k)+w; 
            vvv(k)=vvv(k)+v;
            ddd(k)=ddd(k)+d;
%
            if(x3r(j)<1)
                printf(' error x3r ');
                break;
            end
%
        end
%
end
%
%
pause(0.5);
progressbar(1);
%
t=fix_size(t);
aaa=fix_size(aaa);
vvv=fix_size(vvv);
ddd=fix_size(ddd);
%
acceleration=[t aaa];
    velocity=[t vvv];
displacement=[t ddd];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%