function [MAW0,Fmax,a,MAW1,MAW2] = MeanAtrialWave(X,Fs)

ff = (0:length(X)-1)*Fs/length(X);
YY = abs(fft(X));
ii = find(ff>=3 & ff<=12);
[~,pos] = max(YY(ii));
Fmax = ff(pos+ii(1)-1);
% DFAcendo = mean(Fmax);

if length(X) > 2304
        B = fir1(768,[(Fmax-1.5)/(Fs/2) (Fmax+1.5)/(Fs/2)],'bandpass',chebwin(768+1,40));
        MAW0 = filtfilt(B,1,X);
    
        B = fir1(768,[(2*Fmax-1.5)/(Fs/2) (2*Fmax+1.5)/(Fs/2)],'bandpass',chebwin(768+1,40));
        MAW1 = filtfilt(B,1,X);
    
        B = fir1(768,[(3*Fmax-1.5)/(Fs/2) (3*Fmax+1.5)/(Fs/2)],'bandpass',chebwin(768+1,40));
        MAW2 = filtfilt(B,1,X);
else
    filteroreder = floor(length(X)/3)-3;
        B = fir1(filteroreder,[(Fmax-1.5)/(Fs/2) (Fmax+1.5)/(Fs/2)],'bandpass',chebwin(filteroreder+1,40));
        MAW0 = filtfilt(B,1,X);
    
        B = fir1(filteroreder,[(2*Fmax-1.5)/(Fs/2) (2*Fmax+1.5)/(Fs/2)],'bandpass',chebwin(filteroreder+1,40));
        MAW1= filtfilt(B,1,X);
    
        B = fir1(filteroreder,[(3*Fmax-1.5)/(Fs/2) (3*Fmax+1.5)/(Fs/2)],'bandpass',chebwin(filteroreder+1,40));
        MAW2 = filtfilt(B,1,X);
end

% relative subband energy of the first and second harmonics

% a = (sum((MAW1.^2)')+sum((MAW2.^2)'))./sum((MAW0.^2)');
a = (sum((MAW1.^2)')+sum((MAW2.^2)'))./(sum((MAW1.^2)')+sum((MAW2.^2)')+sum((MAW0.^2)'));
