%
%   vibrationdata_wavelet_decomposition_engine_q_alt.m  ver 1.2  by Tom Irvine
%
%
function[acceleration,wavelet_table,wavelet_low,wavelet_up,alpha,beta]=...
  vibrationdata_wavelet_decomposition_engine_q_alt(t,accel,dt,first,...
                                                   ffmax,nt,nfr,start_time)
%
tp=2*pi;
%
residual=accel;
%
num2=max(size(accel));
%
fprintf(' number of input points= %d \n',num2);
%
duration=t(num2)-t(1);
%
sr=1/dt;
%
fprintf(' sample rate = %10.4g \n\n',sr);
%
fl=3/duration;
fu=sr/10;
%
fprintf(' fl = %10.4g Hz   fu = %10.4g Hz \n\n',fl,fu);
%
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);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   
        [ax1r,ax2r,ax3r,ax4r,alphar,betar,wlr,wur]=...
        decomposition_wgen_q_nw(num2,t,residual,duration,fl,fu,nt,ffmax,first,sr,start_time);    
    
        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);
% 
end
%
fprintf('%d.  ave=%12.4g  sd=%12.4g \n',ie,ave,sd);
pause(0.1);
progressbar(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');

wavelet_table = [(1:nfr)', x2r./tp, x1r, x3r, x4r];

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));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

aaa=zeros(num2,1);

progressbar;
for k=1:num2
%
        progressbar(k/num2);
%
        tt=t(k);
%  
        for j = 1:nfr
            if x2r(j) > 1e-10
                t1 = x4r(j) + t(1);
                t2 = pi * (x3r(j) / x2r(j)) + t1;
        
                if tt >= t1 && tt <= t2
                    arg = x2r(j) * (tt - t1);
                    aaa(k) = aaa(k) + x1r(j) * sin(arg / x3r(j)) * sin(arg);
                end
            end
        end
%
end
%
%
pause(0.1);
progressbar(1);
%
t=fix_size(t);
aaa=fix_size(aaa);

acceleration=[t aaa];