In this study, a native clinical database was used, comprising EEG signals recorded from 200 subjects (100 males and 100 females) with a mean age of 32±832\pm 832±8 years. The subjects were divided into four groups as follows:

I. Healthy group (50 subjects): no history of psychological disorders.

II. Anxiety group (50 subjects): diagnosed with Generalized Anxiety Disorder (GAD) according to DSM-5 criteria III. Chronic stress group (50 subjects): individuals experiencing chronic stress for at least six months IV. Depression group (50 subjects): diagnosed with Major Depressive Disorder (MDD) EEG signals were recorded using a 19-channel system based on the international 10–20 standard, with a sampling frequency of 256 Hz. Each recording session consisted of 10 minutes of resting-state EEG with eyes closed. After removing noise-contaminated segments, a total of 4,071 five-second signal segments were used for training and testing the model. Table 1 presents the final data distribution.

As can be seen, the data distribution is nearly balanced, which helps prevent bias during model training.

The preprocessing of EEG signals includes the following steps:

I. Filtering: the signals were first filtered using a 4th-order Butterworth filter in the frequency range of 0.5–45 Hz to remove low- and high-frequency noise.

II. Artifact removal: to eliminate ocular and muscular artifacts, Independent Component Analysis (ICA) was applied. The components related to artifacts were manually identified and removed.

III. Normalization: each channel was independently normalized using the z-score method:

z = σ . (1)

As can be seen, the data distribution is approximately balanced, which helps prevent bias during model training.

where x denotes the sample value, μ represents the mean, and σ is the standard deviation of the channel.

I. Segmentation: the continuous signals were segmented into 5-second windows (1280 samples) with 50% overlap.

II. Data augmentation: to prevent overfitting and to balance the classes, data augmentation techniques were applied, including the addition of Gaussian noise (SNR = 20 dB), time scaling (±10%), and temporal shifting.

The proposed CNN architecture is designed to extract deep features from one-dimensional EEG signals. This network consists of the following components:

Input layer: a matrix of size 19×128019 (19 channels × 1280 time samples).

Convolutional block 1:

I. Conv1D layer: 64 filters, kernel size = 7, ReLU activation II. Batch Normalization layer III. MaxPooling1D layer with pool size = 2 IV. Dropout layer with rate = 0.3 Convolutional block 2:

I. Conv1D layer: 128 filters, kernel size = 5, ReLU activation II. Batch Normalization layer III. MaxPooling1D layer with pool size = 2 IV. Dropout layer with rate = 0.4 Convolutional block 3:

I. Conv1D layer: 256 filters, kernel size = 3, ReLU activation II. Batch normalization layer III. MaxPooling1D layer with pool size = 2 IV. Dropout layer with rate = 0.5 Fully connected layers:

I. Flatten II. Dense layer: 512 neurons, ReLU activation, Dropout = 0.5 III. Dense layer: 128 neurons, ReLU activation, Dropout = 0.3 IV. Output layer: 4 neurons with Softmax activation The output of this network is a probability vector P = [p1, p2, p3, p4]], which represents the probability that the input signal belongs to each of the four classes (healthy, anxiety, stress, depression). Fig. 1 illustrates the detailed architecture of the proposed CNN.

For each of the four disease classes, a separate HMM is designed. Each HMM is characterized by the following parameters:

λ = (A, B, π), (2)

where:

I. A: state transition probability matrix between hidden states II. B: emission (observation) probability matrix III. π: initial state probability distribution Number of hidden states: for each disorder, five hidden states S={s1,s2,s3,s4,s5} are considered, representing different stages of disease progression. The number of states was determined based on empirical analysis and the Bayesian Information Criterion (BIC) to achieve a balance between model complexity and data fitting.

Transition probability: the probability of transitioning from state i to state j is defined as follows:

aij = P(qt+1 = sj ∣ qt = si). (3)

Observation probability: the probabilistic output of the CNN (P) is fed into the HMM as the observations. Using Bayes’ theorem, the probability of observing the data in state jjj is computed as follows:

P(qt=sj∣Ot)×P(Ot)

bj(Ot) = P(Ot ∣ qt = sj) = P(q =s ) . (4)

HMM training: the parameters of the HMM are trained using the Baum–Welch algorithm (a type of Expectation–Maximization (EM) algorithm). The objective is to maximize the likelihood of the observation sequences:

λ∗ = arg⁡ max λ ⁡P(O ∣ λ). (5)

Inference: for classification, the likelihood of the observation sequence is computed for each of the four HMM models, and the class with the highest likelihood is selected as the final diagnosis.

C∗ = arg⁡ max⁡c∈{1,2,3,4} P(O ∣ λc). (6)

The training of the model is carried out in two stages:

Stage 1. CNN training I. Loss function: categorical cross-entropy II. Optimizer: adam with a learning rate of α=0.0001 III. Batch size: 32 IV. Number of epochs: 100 (with early stopping and patience = 10) V. Data split: 70% training, 15% validation, 15% testing Stage 2. HMM training I. After training the CNN, the probabilistic outputs of the Softmax layer are extracted for the training data II. For each class, one HMM is trained using the Baum–Welch algorithm III. Convergence criterion: change in log-likelihood less than 10-4 Fig. 2 illustrates the overall flow diagram of the proposed system.

The model was implemented using TensorFlow 2.10 and hmmlearn libraries in Python 3.9. The experiments were conducted on a system equipped with an Intel Core i5-11400H processor, 32 GB RAM, and an NVIDIA RTX 3050 Ti GPU.

The CNN training time was approximately 3.5 hours, and the inference time for a single signal segment was about 10 milliseconds, indicating the suitability of the proposed system for real-time applications.

The CNN network converged after 78 epochs. The performance results of the CNN alone, without integration with the HMM, are as follows:

I. Training accuracy: 96.8% II. Validation accuracy: 89.3% III. Test accuracy: 88.7% IV. Final training loss: 0.089 V. Final validation loss: 0.312 Table 2 shows the CNN network evaluation criteria for each class.

Fig. 3 shows the convergence process of the model during training. As can be seen, the model has reached a stable state after epoch 13 and the difference between training and validation accuracy is within acceptable

Four separate HMM models were trained for each class. Each model, with five hidden states, converged after an average of 45 iterations of the Baum–Welch algorithm. The average log-likelihood values for the models are as follows:

I. Healthy HMM model: −245.3 II. Anxiety HMM model: −263.7 III. Stress HMM model: −258.1 IV. Depression HMM model: −251.9

Combining the CNN with the HMM led to a significant improvement in performance. The results of the final hybrid model are as follows:

I. Overall accuracy: 92.5% II. Average class accuracy: 92.3% Table 3 presents the detailed performance metrics of the hybrid model.

The improvement in accuracy from 88.7% (CNN) to 92.5% (CNN–HMM) is statistically significant (p<0.01, paired t-test). This result indicates that adding the HMM layer significantly enhances the system’s classification capability. The confusion matrix of the final model, shown in Fig. 4, reveals that the highest classification error occurs between the “anxiety” and “chronic stress” classes. This finding is clinically justifiable, as these two disorders exhibit substantial similarities in physiological symptoms, particularly in the Beta (13–30 Hz) and Gamma (>30 Hz) frequency bands of EEG signals. Previous studies have also demonstrated that distinguishing between these two disorders is challenging even for clinical specialists.

Table 4 presents a comparison between the proposed model and other methods reported in the literature.

As can be seen, in addition to achieving higher accuracy, the proposed model also exhibits high interpretability, which is crucial for clinical applications. It is worth noting that despite using a larger number of classes (four classes compared to two classes in most studies), our method has still maintained superior accuracy.

One of the main advantages of the HMM model is its interpretability through transition matrices. Fig. 5 illustrates the transition matrices for the four HMM models (one for each class). For example, the transition matrix of the HMM corresponding to the depression class indicates that:

I. The probability of remaining in the early states (State 1 and State 2), which represent low disease severity, is relatively low (approximately 0.65–0.70).

II. The probability of transitioning to more advanced states (State 3, State 4, and State 5) is relatively high (approximately 0.25–0.30).

This pattern is consistent with the progressive worsening and temporal evolution of depressive symptoms. In contrast, the transition matrix for the healthy class shows that the probability of remaining in the initial states is very high (around 0.85), while the probability of transitioning to other states is very low (less than 0.10), indicating stability and lack of change in mental health status.

These transition matrices can be used as an auxiliary diagnostic tool for clinicians, helping them better understand the progression of the disease over time and make more informed and appropriate treatment decisions.

Fig. 6 shows a visual comparison of the performance of the CNN model against the CNN-HMM hybrid model for various metrics. As can be seen, the addition of the HMM uniformly improved all the evaluation