Problem with identifying the entropy stabilization point in SSA for accurate signal decomposition.

In Singular Spectrum Analysis (SSA), the goal is to decompose a signal into three main components:

1. Trend: The overall trend of the signal, which represents long-term changes.
2. Oscillations: Parts of the signal that change periodically, typically corresponding to medium-term variations.
3. Noise: Parts of the signal that change randomly, usually representing noise or disturbances in the signal.

In this process, the signal is first transformed into a trajectory matrix, and then its components are extracted via Singular Value Decomposition (SVD).

Main Challenge:

• After performing SVD and calculating the entropy of each component, we are looking for a point where the entropy stabilizes. This point is considered where additional components no longer contribute meaningful information to the model.
• The problem is that the point of entropy stabilization is not being correctly identified, which results in an improper decomposition of the signal into the three components (trend, oscillations, and noise).

Question: Is there a method in MATLAB that can more accurately identify the point of entropy stabilization and correctly decompose the signal into trend, oscillations, and noise?

clc;

clear;

close all;

% Load the CSV file

[file, path] = uigetfile('*.csv', 'Please upload your CSV file:');

if isequal(file, 0)

error('No file selected.');

end

data = readtable(fullfile(path, file));

% Display available columns for the user to choose from

disp('Available columns in your dataset:');

disp(data.Properties.VariableNames');

column_name = input('Enter the column name you want to analyze: ', 's');

% Extract the signal (time series) from the chosen column

if ~any(strcmp(data.Properties.VariableNames, column_name))

error(['Column "', column_name, '" not found in the dataset.']);

end

signal = data.(column_name);

signal = signal(~isnan(signal)); % Drop NaN values

% SSA Functions

function X = trajectory_matrix(signal, window_length)

N = length(signal);

K = N - window_length + 1;

X = zeros(window_length, K);

for i = 1:K

X(:, i) = signal(i:i + window_length - 1);

end

end

function entropy = singular_entropy(S)

P = S / sum(S);

entropy = -sum(P .* log(P));

end

function reconstructed_signal = reconstruct_signal(U, S, V, indices, original_length)

reconstructed_trajectory = zeros(size(U, 1), size(V, 1));

for i = indices

reconstructed_trajectory = reconstructed_trajectory + S(i) * (U(:, i) * V(i, :));

end

[rows, cols] = size(reconstructed_trajectory);

signal_length = rows + cols - 1;

reconstructed_signal = zeros(signal_length, 1);

counts = zeros(signal_length, 1);

% Extracting anti-diagonals

for k = 1:(rows + cols - 1)

anti_diag_elements = [];

for i = 1:rows

j = k - i + 1;

if j >= 1 && j <= cols

anti_diag_elements = [anti_diag_elements; reconstructed_trajectory(i, j)];

end

end

reconstructed_signal(k) = mean(anti_diag_elements);

counts(k) = length(anti_diag_elements);

end

reconstructed_signal = reconstructed_signal ./ counts; % Normalize the signal

reconstructed_signal = reconstructed_signal(1:original_length);

end

% Perform SSA

window_length = floor(length(signal) / 3); % Adjust window length

X = trajectory_matrix(signal, window_length);

% Perform SVD decomposition

[U, S_mat, V] = svd(X, 'econ');

S = diag(S_mat);

% Calculate Singular Entropy for each component cumulatively

cumulative_entropy = zeros(length(S), 1);

for i = 1:length(S)

cumulative_entropy(i) = singular_entropy(S(1:i));

end

% Plot Singular Entropy to observe stabilization point

figure;

plot(1:length(S), cumulative_entropy, '-o', 'LineWidth', 1.5);

xlabel('Component Index');

ylabel('Cumulative Singular Entropy');

title('Cumulative Singular Entropy vs. Component Index');

grid on;

% Calculate the point where entropy stabilizes

entropy_diff = diff(cumulative_entropy); % Difference between consecutive entropy values

threshold = 0.01; % Threshold to detect stabilization

stable_point = find(entropy_diff < threshold, 1);

% Mark the stabilization point

hold on;

xline(stable_point + 1, '--r', ['Stabilization Point: ', num2str(stable_point + 1)]);

legend('Cumulative Singular Entropy', 'Stabilization Point');

% Number of components based on entropy stabilization

num_components = stable_point + 1;

disp(['Number of components based on entropy stabilization: ', num2str(num_components)]);

% Reconstruct the signal by categorizing components

trend_indices = 1; % Trend corresponds to the first component

oscillation_indices = 2:num_components; % Oscillations from second component to stabilization point

noise_indices = (num_components + 1):length(S); % Noise corresponds to components after stabilization

% Reconstruct signals based on selected indices

trend_signal = reconstruct_signal(U, S, V', trend_indices, length(signal));

oscillation_signal = reconstruct_signal(U, S, V', oscillation_indices, length(signal));

noise_signal = reconstruct_signal(U, S, V', noise_indices, length(signal));

% Plot the original signal, Trend, Oscillations, and Noise

figure;

plot(signal, 'b', 'DisplayName', 'Original Signal', 'LineWidth', 1.5);

hold on;

plot(trend_signal, 'g', 'DisplayName', 'Trend', 'LineWidth', 1.5);

plot(oscillation_signal, 'orange', 'DisplayName', 'Oscillations', 'LineWidth', 1.5);

plot(noise_signal, 'r', 'DisplayName', 'Noise', 'LineWidth', 1.5);

title('Signal Decomposition');

xlabel('Time');

ylabel('Amplitude');

legend;

grid on;

Answer 1

[William Rose]

@Trais

You say SSA decomposes a signal into three components: trend, oscillations, and noise. Shouldn't the sum of the three components equal the signal? Even if the point of entropy stabilization is not found correctly? In your second plot, the three components do not add up to the orignal signal. I would fix that before working on a better way to find the point of entropy stabilization.

When I run your code with your data, I get an error - see below.

I notice that reconstruct_signal is defined

function reconstructed_signal = reconstruct_signal(U, S, V, indices, original_length)

and that you pass V', not V, when you call it. But changing V' to V in the call did not fix the error. I would check vector and array dimensions in reconstruct_signal(). Then we'll get to the point of entropy stabilization.

More later.

% Load the CSV file
data = importdata('VAL320.csv');
signal=data.data;
signal = signal(~isnan(signal)); % Drop NaN values
% SSA Functions
function X = trajectory_matrix(signal, window_length)
N = length(signal);
K = N - window_length + 1;
X = zeros(window_length, K);
for i = 1:K
    X(:, i) = signal(i:i + window_length - 1);
end
end
function entropy = singular_entropy(S)
P = S / sum(S);
entropy = -sum(P .* log(P));
end

function reconstructed_signal = reconstruct_signal(U, S, V, indices, original_length)
reconstructed_trajectory = zeros(size(U, 1), size(V, 1));
for i = indices
    reconstructed_trajectory = reconstructed_trajectory + S(i) * (U(:, i) * V(i, :));
end
[rows, cols] = size(reconstructed_trajectory);
    signal_length = rows + cols - 1;
    reconstructed_signal = zeros(signal_length, 1);
    counts = zeros(signal_length, 1);
    % Extracting anti-diagonals
    for k = 1:(rows + cols - 1)
        anti_diag_elements = [];
        for i = 1:rows
            j = k - i + 1;
            if j >= 1 && j <= cols
                anti_diag_elements = [anti_diag_elements; reconstructed_trajectory(i, j)];
            end
        end
        reconstructed_signal(k) = mean(anti_diag_elements);
        counts(k) = length(anti_diag_elements);
    end
    reconstructed_signal = reconstructed_signal ./ counts; % Normalize the signal
    reconstructed_signal = reconstructed_signal(1:original_length);
end

% Perform SSA
window_length = floor(length(signal) / 3); % Adjust window length
X = trajectory_matrix(signal, window_length);
% Perform SVD decomposition
[U, S_mat, V] = svd(X, 'econ');
S = diag(S_mat);
% Calculate Singular Entropy for each component cumulatively
cumulative_entropy = zeros(length(S), 1);
for i = 1:length(S)
    cumulative_entropy(i) = singular_entropy(S(1:i));
end
% Plot Singular Entropy to observe stabilization point
figure;
plot(1:length(S), cumulative_entropy, '-o', 'LineWidth', 1.5);
xlabel('Component Index');
ylabel('Cumulative Singular Entropy');
title('Cumulative Singular Entropy vs. Component Index');
grid on;
% Calculate the point where entropy stabilizes
entropy_diff = diff(cumulative_entropy); % Difference between consecutive entropy values
threshold = 0.01; % Threshold to detect stabilization
stable_point = find(entropy_diff < threshold, 1);
% Mark the stabilization point
hold on;
xline(stable_point + 1, '--r', ['Stabilization Point: ', num2str(stable_point + 1)]);
legend('Cumulative Singular Entropy', 'Stabilization Point');

% Number of components based on entropy stabilization
num_components = stable_point + 1;
disp(['Number of components based on entropy stabilization: ', num2str(num_components)]);

Number of components based on entropy stabilization: 64

% Reconstruct the signal by categorizing components
trend_indices = 1; % Trend corresponds to the first component
oscillation_indices = 2:num_components; % Oscillations from second component to stabilization point
noise_indices = (num_components + 1):length(S); % Noise corresponds to components after stabilization
% Reconstruct signals based on selected indices
trend_signal = reconstruct_signal(U, S, V', trend_indices, length(signal));

Arrays have incompatible sizes for this operation.

Error in solution>reconstruct_signal (line 22)
reconstructed_trajectory = reconstructed_trajectory + S(i) * (U(:, i) * V(i, :));

oscillation_signal = reconstruct_signal(U, S, V', oscillation_indices, length(signal));
noise_signal = reconstruct_signal(U, S, V', noise_indices, length(signal));
% Plot the original signal, Trend, Oscillations, and Noise
figure;
plot(signal, 'b', 'DisplayName', 'Original Signal', 'LineWidth', 1.5);
hold on;
plot(trend_signal, 'g', 'DisplayName', 'Trend', 'LineWidth', 1.5);
plot(oscillation_signal, 'orange', 'DisplayName', 'Oscillations', 'LineWidth', 1.5);
plot(noise_signal, 'r', 'DisplayName', 'Noise', 'LineWidth', 1.5);
title('Signal Decomposition');
xlabel('Time');
ylabel('Amplitude');
legend;
grid on;

OK

[Trais]

Here are the steps of Singular Spectrum Analysis (SSA) in a structured and concise manner:

Steps of Singular Spectrum Analysis (SSA)

Embedding:

• Transform the original time series X={x1,x2,…,xN} into a trajectory matrix.
• Choose a window length L such that 2≤L≤N.
• The trajectory matrix X\mathbf{X}X is constructed as follows: \mathbf{X} = \begin{bmatrix} x_1 & x_2 & \cdots & x_K \\ x_2 & x_3 & \cdots & x_{K+1} \\ \vdots & \vdots & \ddots & \vdots \\ x_L & x_{L+1} & \cdots & x_N \end{bmatrix}​​​
• where K=N−L+1.

Decomposition:

• Perform Singular Value Decomposition (SVD) on the trajectory matrix: \mathbf{X} = \mathbf{U} \mathbf{S} \mathbf{V}^\top
• U: Left singular vectors.
• S: Diagonal matrix of singular values.
• \mathbf{V}: Right singular vectors.

Reconstruction:

• Reconstruct the signal by grouping singular values and corresponding components: \mathbf{X}_i = s_i \cdot \mathbf{U}_i \mathbf{V}_i^\top​
• Each \mathbf{X}_i​ corresponds to a component of the signal (trend, oscillations, or noise).

Interpretation:

• Analyze the reconstructed components:
• Trend: Represents the long-term behavior of the signal.
• Oscillations: Captures cyclical patterns.
• Noise: Contains random or less significant components.

[Trais]

Here is the Python code:

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from google.colab import files

# Upload the file in Google Colab

print("Please upload your CSV file:")

uploaded = files.upload()

# Load the CSV file into a pandas dataframe

file_name = list(uploaded.keys())[0]

data = pd.read_csv(file_name)

# Display available columns for the user to choose from

print("\nAvailable columns in your dataset:")

print(data.columns.tolist())

column_name = input("\nEnter the column name you want to analyze: ")

# Extract the signal (time series) from the chosen column

if column_name not in data.columns:

raise ValueError(f"Column '{column_name}' not found in the dataset.")

signal = data[column_name].dropna().values # Drop NaN values

# SSA Functions

def trajectory_matrix(signal, window_length):

N = len(signal)

K = N - window_length + 1

return np.array([signal[i:i + window_length] for i in range(K)]).T

def singular_entropy(S):

P = S / np.sum(S)

return -np.sum(P * np.log(P))

def reconstruct_signal(U, S, Vt, indices, original_length):

reconstructed_trajectory = np.zeros_like(U @ np.diag(S) @ Vt)

for i in indices:

reconstructed_trajectory += np.outer(U[:, i], S[i] * Vt[i, :])

rows, cols = reconstructed_trajectory.shape

signal_length = rows + cols - 1

reconstructed_signal = np.zeros(signal_length)

counts = np.zeros(signal_length)

# Extracting anti-diagonals (from top-right to bottom-left)

for k in range(rows + cols - 1):

anti_diag_elements = []

for i in range(rows):

j = k - i

if 0 <= j < cols:

anti_diag_elements.append(reconstructed_trajectory[i, j])

# Average of the anti-diagonal elements

reconstructed_signal[k] = np.mean(anti_diag_elements)

counts[k] = len(anti_diag_elements)

reconstructed_signal /= counts # Normalize the signal

return reconstructed_signal[:original_length]

# Perform SSA

window_length = int(len(signal) / 3) # Adjust window length

X = trajectory_matrix(signal, window_length)

# Perform SVD decomposition

U, S, Vt = np.linalg.svd(X, full_matrices=False)

# Calculate Singular Entropy for each component cumulatively

cumulative_entropy = []

for i in range(1, len(S) + 1):

entropy = singular_entropy(S[:i])

cumulative_entropy.append(entropy)

# Plot Singular Entropy to observe stabilization point

plt.figure(figsize=(10, 6))

plt.plot(range(1, len(S) + 1), cumulative_entropy, marker='o', label='Cumulative Singular Entropy')

plt.xlabel('Component Index')

plt.ylabel('Cumulative Singular Entropy')

plt.title('Cumulative Singular Entropy vs. Component Index')

plt.legend()

plt.show()

# Calculate the point where entropy stabilizes

entropy_diff = np.diff(cumulative_entropy) # Calculate difference between consecutive entropy values

threshold = 0.01 # Threshold to detect stabilization

stable_point = np.argmax(entropy_diff < threshold) # Find the point where entropy stabilizes

# Mark the stabilization point

plt.figure(figsize=(10, 6))

plt.plot(range(1, len(S) + 1), cumulative_entropy, marker='o', label='Cumulative Singular Entropy')

plt.axvline(x=stable_point + 1, color='r', linestyle='--', label=f'Stabilization Point: {stable_point + 1}')

plt.xlabel('Component Index')

plt.ylabel('Cumulative Singular Entropy')

plt.title('Cumulative Singular Entropy vs. Component Index')

plt.legend()

plt.show()

# Number of components based on entropy stabilization

num_components = stable_point + 1

print(f"Number of components based on entropy stabilization: {num_components}")

# Reconstruct the signal by categorizing components

trend_indices = [0] # Trend corresponds to the first component

oscillation_indices = list(range(1, num_components)) # Oscillations are from second component to stabilization point

noise_indices = list(range(num_components, len(S))) # Noise corresponds to the components after stabilization

# Reconstruct signals based on selected indices

trend_signal = reconstruct_signal(U, S, Vt, trend_indices, len(signal))

oscillation_signal = reconstruct_signal(U, S, Vt, oscillation_indices, len(signal))

noise_signal = reconstruct_signal(U, S, Vt, noise_indices, len(signal))

# Plot the original signal, Trend, Oscillations, and Noise

plt.figure(figsize=(12, 8))

plt.plot(signal, label='Original Signal', color='blue')

plt.plot(trend_signal, label='Trend', color='green')

plt.plot(oscillation_signal, label='Oscillations', color='orange')

plt.plot(noise_signal, label='Noise', color='red')

plt.title('Signal Decomposition')

plt.xlabel('Time')

plt.ylabel('Amplitude')

plt.legend()

plt.show()

[William Rose]

@Trais,

I can't fix the error at

reconstructed_trajectory=reconstructed_trajectory+...;

in the main program until I know exactly what you are trying to do there. (reconstructed_trajectory=... in the function works fine.) You don't use the result anyway, so I commented it out in the main program.

There's still the issue of whether to use V or V', and getting the dimensions right for the array multiplication. You originally had

function reconstructed_signal = reconstruct_signal(U, S, V, indices, original_length)

reconstructed_trajectory = zeros(size(U, 1), size(V, 1));

for i = indices

reconstructed_trajectory = reconstructed_trajectory + S(i) * (U(:, i) * V(i,:));

end

...

end

and then in your main program you call reconstruct_signal() with V', not V. There are array size errors. One way to fix it is to change the function to

function reconstructed_signal = reconstruct_signal(U, S, V, indices, original_length)

reconstructed_trajectory = zeros(size(U, 1),size(V, 1));

for i = indices

reconstructed_trajectory = reconstructed_trajectory + S(i)*(V(i,:)*U(:, i));

end

...

end

and pass V, not V', each time you call reconstruct_signal(). This cures the error, but gives weird results. See below. How do you want to fix it?

% Load the CSV file
data = importdata('VAL320.csv');
signal=data.data;
signal = signal(~isnan(signal)); % Drop NaN values
% SSA Functions
function X = trajectory_matrix(signal, window_length)
N = length(signal);
K = N - window_length + 1;
X = zeros(window_length, K);
for i = 1:K
    X(:, i) = signal(i:i + window_length - 1);
end
end

function entropy = singular_entropy(S)
P = S / sum(S);
entropy = -sum(P .* log(P));
end

function reconstructed_signal = reconstruct_signal(U, S, V, indices, original_length)
reconstructed_trajectory = zeros(size(U, 1),size(V, 1));
for i = indices
    reconstructed_trajectory = reconstructed_trajectory + S(i)*(V(i, :)*U(:, i));
end
[rows, cols] = size(reconstructed_trajectory);
signal_length = rows + cols - 1;
reconstructed_signal = zeros(signal_length, 1);
counts = zeros(signal_length, 1);
% Extracting anti-diagonals
for k = 1:(rows + cols - 1)
    anti_diag_elements = [];
    for i = 1:rows
        j = k - i + 1;
        if j >= 1 && j <= cols
            anti_diag_elements = [anti_diag_elements; reconstructed_trajectory(i, j)];
        end
    end
    reconstructed_signal(k) = mean(anti_diag_elements);
    counts(k) = length(anti_diag_elements);
end
reconstructed_signal = reconstructed_signal ./ counts; % Normalize the signal
reconstructed_signal = reconstructed_signal(1:original_length);
end

% Perform SSA
window_length = floor(length(signal) / 3); % Adjust window length
X = trajectory_matrix(signal, window_length);
% Perform SVD decomposition
[U, S_mat, V] = svd(X, 'econ');
S = diag(S_mat);
% Calculate Singular Entropy for each component cumulatively
cumulative_entropy = zeros(length(S), 1);
for i = 1:length(S)
    cumulative_entropy(i) = singular_entropy(S(1:i));
end
% Plot Singular Entropy to observe stabilization point
figure;
plot(1:length(S), cumulative_entropy, '-o', 'LineWidth', 1.5);
xlabel('Component Index');
ylabel('Cumulative Singular Entropy');
title('Cumulative Singular Entropy vs. Component Index');
grid on;
% Calculate the point where entropy stabilizes
entropy_diff = diff(cumulative_entropy); % Difference between consecutive entropy values
threshold = 0.01; % Threshold to detect stabilization
stable_point = find(entropy_diff < threshold, 1);
% Mark the stabilization point
hold on;
xline(stable_point + 1, '--r', ['Stabilization Point: ', num2str(stable_point + 1)]);
legend('Cumulative Singular Entropy', 'Stabilization Point');

% Number of components based on entropy stabilization
num_components = stable_point + 1;
disp(['Number of components based on entropy stabilization: ', num2str(num_components)]);

Number of components based on entropy stabilization: 64

% Reconstruct the signal by categorizing components
trend_indices = 1; % Trend corresponds to the first component
oscillation_indices = 2:num_components; % Oscillations from second component to stabilization point
noise_indices = (num_components + 1):length(S); % Noise corresponds to components after stabilization
% Reconstruct signals based on selected indices
trend_signal = reconstruct_signal(U, S, V, trend_indices, length(signal));
% reconstructed_trajectory = reconstructed_trajectory + S(i) * (U(:, i) * V(i, :));
oscillation_signal = reconstruct_signal(U, S, V, oscillation_indices, length(signal));
noise_signal = reconstruct_signal(U, S, V, noise_indices, length(signal));
% Plot the original signal, Trend, Oscillations, and Noise
figure;
plot(signal, 'b', 'DisplayName', 'Original Signal', 'LineWidth', 1.5);
hold on;
plot(trend_signal, 'g', 'DisplayName', 'Trend', 'LineWidth', 1.5);
plot(oscillation_signal, 'm', 'DisplayName', 'Oscillations', 'LineWidth', 1.5);
plot(noise_signal, 'r', 'DisplayName', 'Noise', 'LineWidth', 1.5);
title('Signal Decomposition');
xlabel('Time');
ylabel('Amplitude');
legend;
grid on;

[William Rose]

@Trais,

[Edit: I deleted a comment below ("I would use the transpose..."). This comment reflected a lack of understanding.]

From reading your code I am starting to understand SSA.

SSA does a kind of pricipal component analysis on signal Y. It does this by creating a matrix X, whose columns are windowed chunks of Y. Let N be the length of Y. Define window length L=floor(N/3). This seems a bit arbitrary, but it is what it is. Matrix X has size L x (N+1-L). Each column of X is a windowed chunk of Y: X(:,1)=Y(1:L), X(:,2)=Y(2:L+1), etc. (X defined this way, with L~=N/3, will have approximately twice as many columns as rows.)

Do SVD on X: [U,S,V]=svd(X,'econ'). Then S is LxL, U is LxL, and V is (N+1-L)xL. This is kind of like doing principal component analysis (PCA) on X, but in PCA, the row vector mu=mean(X) is computed and is subtracted form every row of X, before doing SVD. Another practical distinction is, when doing PCA on X in Matlab, each row of X represents one observation of the variable or process in question. @Trais, your matrix X has the observations (the windowed chunks of the signal) as columns, not rows.

An entropy-like quantity is computed for each "principal component". (In quotes because it is not comuted in the usual way.) The stabilization point is the index of the last "principal component" which improves the cumulative entropy by more than threshold. (threshold equals 0.01. Why this value?)

The first PC is assumed to represent the trend. PCs 2 through the stabilization point are assumed to represent "oscillations" in the signal. The remaining PCs are assumed to correspond to noise.

The code attempts to use the PCs to reconstruct the signal, but it is not working right.

I would do the reconstruction differently.

I would make simulated data with clearly defined trend and oscillations, to test the algorithm.

[William Rose]

@Trais,

Attached code does SSA. The plot of the signal and its components and the reconstruction look reasonable. I have not tried to figure out exactly why your code does not produce a reasonable looking result.

My script does not find the stabilization point. But you can easily modify the script to do so, since you have all the singular values, and you can compute the entropies* of the components.

The plot above shows the contribution to variance of the principal components found by SSA. The y-axis scale in the lower plot above shows that the first component accounts for 95% of the variance in the signal. The large fraction of "variance" accounted for by component 1 is due in part to the fact that the signal is not zero mean. (Which is why I put variance in quotes. Maybe "power" would be a more accurate term for this non-zero-mean signal.) The plot of cumulative variance seems to have a shoulder, or breakpoint, after 6 or 7 components.

If I were trying to find the stabilization point, I would remove the mean value before doing the SSA, so that the first component is not so huge. Then the remaining components would be larger, and it might make it easier to separate oscillation components from noise components.

The next plot shows the signal and its components.

The plot above shows the signal x, the exact reconstructed signal xre (which uses all the principal components), and the first six individual components. The first component, in red, appears to be a reasonable estimate of trend. The plot shows that x and xre are exactly the same, as we expect.

* The entropy in a dataset, is , i.e. the entropy is the sum of all the individual components. It is not accurate to refer to the "entropy" of individual components. But I understand why you do, and maybe the quantity p*log(p) will be a good way to find the stabilization point.

[Trais]

Thank you so much for the time you've dedicated and the helpful explanations you've provided. Your insights have been invaluable in better understanding SSA and improving the code for better results. I truly appreciate your time and attention.

[Torsten]

@Trais

If the answer could help you, you should accept it to show your appreciation.

[William Rose]

@Trais, you're welcome.

[William Rose]

@Trais,

Do you have a reference explaining the use and justification for the entropy stabilization point in SSA? I am asking because the book Singular Spectrum Analysis for Time Series, 2nd edition (2020) by Golyandina and Zhigljavsky, does not mention the entropy stabilization point. This book, or its first edition, is a commonly cited authority for SSA. The authors do mentin other methods, whcih I was also thinking about before I got the book: use Akaike information criterion, or Bayes information criterion (AIC or BIC), to decide where to draw the line for noise. I used AIC to select the number of model components in this publication.

One could also take a graphical approach to decide which components are oscillation and which are noise. See image below, which uses real experimental data.

Source: William Rose. Biomedical Signal Processing Lecture Notes, 2024. Concept from D. Winter, Biomechanics, 3rd edition, 2004.