function [norm_E1, norm_E] = wavelet_high_low_frequencies_norm(temp,fs)

% WAVELET_HIGH_LOW_FREQUENCIES_NORM Computes the normalized energy distribution 
% of a signal using Discrete Wavelet Transform (DWT) decomposition.
%
% Inputs:
%   temp - Input time-series signal (e.g., Center of Pressure data in ML or AP)
%   fs   - Sampling frequency (Hz)
%
% Outputs:
%   norm_E1 - 1x2 vector containing dynamically pooled energy percentages:
%             [1] High-frequency component (dynamically bound between (approximately) 10 Hz and 1 Hz)
%             [2] Low-frequency component (from the sub-1 Hz level down to the approximation)
%   norm_E  - 1x(levels+1) vector containing normalized energy percentages of all 
%             individual detail levels and the final approximation level.

%% Parameters & Initialization

N = length(temp);         % Number of data samples
time = 0 :1/fs:(N-1)/fs;  % Time vector corresponding to the data
levels = floor(log2(N));  % Maximum allowable wavelet decomposition levels
wname = 'db2';            % Daubechies wavelet family type

%% Frequency Band Estimation

% Tracks cutoff frequency boundaries for each decomposition level
frequencies = zeros(levels,2); % Col 1: lower-bound (HP); Col 2: upper-bound (LP)
frequencies(1,1) = fs/4;
frequencies(1,2) = fs/2;
for i = 2 : levels
    frequencies(i,2) = frequencies(i-1,1);
    frequencies(i,1) = frequencies(1,2)/(2^i);
end

%% Wavelet Decomposition

% Perform discrete wavelet decomposition
[C,L] = wavedec(temp,levels,wname);

% Reconstruct the final approximation signal from decomposition coefficients
A = wrcoef('a', C, L, wname);

% Reconstruct specific detail signals for each decomposition level
D = {};
for i = 1 : levels
    D(i) = {wrcoef('d',C,L,wname,i)};
end

%% Energy Content Calculation

% Compute absolute energy (sum of squares) for each detail level
E_detail = zeros(1, levels);
for k = 1:levels
    E_detail(k) = sum(D{k}.^2);
end

% Compute absolute energy of the approximation signal
E_appr = sum(A.^2);

%% Total Energy Calculation

% Calculate the cumulative energy across all reconstructed components
Et = E_appr + sum(E_detail);

%% Normalized Energy Distribution (%)

% Normalize energy levels as percentages of the total signal energy
norm_E_detail = (E_detail./Et)*100;
norm_E_appr = (E_appr./Et)*100;

norm_E = [norm_E_detail norm_E_appr];

%% Dynamic Band Boundary Identification

% Find detail levels for the High-Frequency band (frequency < 10 Hz and >= 1 Hz)
high_freq_indices = find(frequencies(:,1) < 10 & frequencies(:,2) >= 1);
start_high = high_freq_indices(1);
end_high = high_freq_indices(end);

% The Low-Frequency band starts right after the high-frequency band 
% and spans until the final index of norm_E (which includes the approximation level)
start_low = end_high + 1;
end_low = length(norm_E);

%% Pooled Band Energy (High vs. Low Frequencies)

% Sum the percentage energy contributions across the dynamically found groups
norm_E1(1) = sum(norm_E(start_high : end_high)); % Sum of the High-frequency components
norm_E1(2) = sum(norm_E(start_low : end_low));   % Sum of the Low-frequency components

end