Vibroarthrographic Signal Spectral Features in 5-Class Knee Joint Classification

Vibroarthrography (VAG) is a non-invasive and potentially widely available method supporting the joint diagnosis process. This research was conducted using VAG signals classified to five different condition classes: three stages of chondromalacia patellae, osteoarthritis, and control group (healthy knee joint). Ten new spectral features were proposed, distinguishing not only neighboring classes, but every class combination. Additionally, Frequency Range Maps were proposed as the frequency feature extraction visualization method. The results were compared to state-of-the-art frequency features using the Bhattacharyya coefficient and the set of ten different classification algorithms. All methods evaluating proposed features indicated the superiority of the new features compared to the state-of-the-art. In terms of Bhattacharyya coefficient, newly proposed features proved to be over 25% better, and the classification accuracy was on average 9% better.


Introduction
The knee is one of the most loaded joints within the human body, highly susceptible to injuries and an increased risk of early degeneration of the articular surfaces. Classical radiography is a basic diagnostic tool for imaging knee injuries. In advanced degeneration of this joint, X-ray examination correlates with arthroscopy evaluation which is used as a "gold standard". However, the lower sensitivity and specificity of X-rays are limitations for diagnosis of early stages of chondral disorders, e.g., chondromalacia. On the other hand, the availability of modern imaging methods such as magnetic resonance imaging is limited due to high expense. One of the experimental methods developed for sensitive assessment of articular function is vibroarthrography (VAG) [1][2][3]. This relatively cheap and potentially widely available non-invasive method is based on the analysis of high frequency vibroacoustic emission, which is a natural phenomenon acquired from the relative motion of articular surfaces of the synovial joint. Although the VAG method is still under development, it reveals high accuracy, sensitivity, and specificity. Previously, it has been used to differentiate disorders of the patellofemoral joint, due to the specific, disorder-related character of the VAG signal pattern [4][5][6].
Exemplary VAG signals, generated by knee joints with different conditions were presented in Figure 1. Conditions included in the Figure and in the rest of this paper are as follows: control group (healthy knee joint, abbreviated as ctrl), three stages of chondromalacia patellae (cmp1, cmp2, and cmp3, consecutively), and osteoarthritis (oa). As it can be seen, a lot of information about the specific condition could potentially be embedded in the frequency spectrum of the signal. Although frequency analysis is present in literature devoted to vibroarthrography, the state-ofthe-art spectral features are designed to distinguish all classes at once [1][2][3][4][5][6], i.e., to utilize a single value to differentiate between specific conditions. The objective of this study was to find better frequency ranges, focusing on the distinction between each particular class pairs (i.e., ctrl-cmp1, ctrl-cmp2, ctrl-cmp3, …, cmp3-oa).
This results in ten spectral features, instead of one. There are ten spectral features, as there are two classes in a pair and five distinct classes, i.e., 5 choose 2 combinations. Specific combinations were presented in Tables 1-3. That way, a single numeric value was changed to a 10-element feature vector, distinguishing conditions more precisely, ensuring that each class is described by the feature set as unambiguously as possible. That approach prevents the potential problem of some values of the feature being typical for more than one condition, making it impossible to state a precise diagnosis.
To amplify the higher frequencies in a way and enable better distinction between classes with potentially useful information hidden in wider frequency spectrum, analysis of the derivative of the signal was conducted.

Standard Descriptors
Wu in [7] specified three main methods of the knee joint VAG signal analysis: the spatiotemporal, time-frequency, and statistical analysis methods. Using the spatiotemporal analysis method Wu recommends to use adaptive segmentation [8,9] to avoid redundancy in given segments, and emphasizes its usefulness for calculating sets of features and signal classification. The statistical analysis is often used to determine distribution measures and statistical parameters [10,11], giving the possibility to display data in the form of a tabular summary [12] and graphic forms such as the histogram [6] or the box plot [13][14][15], and to analyze variations using various types of tests (t-test [16], Kruskal-Wallis test [13], one-way ANOVA test [4], Wilcoxon rank-sum test [11]). The time-frequency analysis and the frequency analysis supported by using statistical analysis are commonly used in the field of signal processing. The frequently used VAG signal analysis is the short-time Fourier transform (STFT) extended by the visual representation of a spectrogram. Łysiak et al. [14] used STFT and spectrograms, on whose specific analysis the new three descriptors have been proposed. Dołęgowski et al. [15] proposed the incremental decomposition of voltage in time and the spectrogram as the methods to identify the knee joint disease stage and compared the statistical Although frequency analysis is present in literature devoted to vibroarthrography, the state-of-the-art spectral features are designed to distinguish all classes at once [1][2][3][4][5][6], i.e., to utilize a single value to differentiate between specific conditions. The objective of this study was to find better frequency ranges, focusing on the distinction between each particular class pairs (i.e., ctrl-cmp1, ctrl-cmp2, ctrl-cmp3, . . . , cmp3-oa).
This results in ten spectral features, instead of one. There are ten spectral features, as there are two classes in a pair and five distinct classes, i.e., 5 choose 2 combinations. Specific combinations were presented in Tables 1-3. That way, a single numeric value was changed to a 10-element feature vector, distinguishing conditions more precisely, ensuring that each class is described by the feature set as unambiguously as possible. That approach prevents the potential problem of some values of the feature being typical for more than one condition, making it impossible to state a precise diagnosis.
To amplify the higher frequencies in a way and enable better distinction between classes with potentially useful information hidden in wider frequency spectrum, analysis of the derivative of the signal was conducted.

Standard Descriptors
Wu in [7] specified three main methods of the knee joint VAG signal analysis: the spatiotemporal, time-frequency, and statistical analysis methods. Using the spatiotemporal analysis method Wu recommends to use adaptive segmentation [8,9] to avoid redundancy in given segments, and emphasizes its usefulness for calculating sets of features and signal classification. The statistical analysis is often used to determine distribution measures and statistical parameters [10,11], giving the possibility to display data in the form of a tabular summary [12] and graphic forms such as the histogram [6] or the box plot [13][14][15], and to analyze variations using various types of tests (t-test [16], Kruskal-Wallis test [13], one-way ANOVA test [4], Wilcoxon rank-sum test [11]). The time-frequency analysis and the frequency analysis supported by using statistical analysis are commonly used in the field of signal processing. The frequently used VAG signal analysis is the short-time Fourier transform (STFT) extended by the visual representation of a spectrogram. Łysiak et al. [14] used STFT and spectrograms, on whose specific analysis the new three descriptors have been proposed. Dołęgowski et al. [15] proposed the incremental decomposition of voltage in time and the spectrogram as the methods to identify the knee joint disease stage and compared the statistical parameters of normalized energy values of the band 50-250 Hz (P1) and 250-450 Hz (P2). Befrui et al. [2] in their measurement system used the two accelerometers, the piezoelectric disk and the potentiometer (four channels) in Sensors 2020, 20, 5015 3 of 15 their measurement system. The signals were acquired at 16 kHz sampling frequency and extension and flexion cycles were extracted by using semi-automatic segmentation. The power spectra were calculated, frequency features were normalized and averaged, and the classification by a linear support vector machine (SVM) using the knee-specific feature vectors was performed.
Wavelet transformation is another tool which gives the possibility to analyze the change of signal frequency in the function of time. Mascarenhas et al. [17], for the analysis of VAG signals, proposed in their paper the tunable Q wavelet transform (TQWT) in comparison to the traditional wavelet packet decomposition (WPD). To overcome the imbalance set problem they used the synthetic minority oversampling technique (SMOTE) and afterwards they compared the performance of the random forest classifier (RF) they used to the soft margin support vector machine (SVM).
Nalband et al. [18] proposed for the VAG signals analysis a CAD system using time-frequency analysis and nonstationary signal processing techniques. Their methods were the Wavelet packet decomposition algorithm (WPD), the smoothed pseudo Wigner-Ville distribution (SPWVD) as a nonstationary time-frequency analysis, and a modified version of Hilbert-Huang transform (HHT) where instead of empirical mode decomposition (EMD) for computing intrinsic mode functions (IMF) they proposed complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN). The Least square support vector machine (LS-SVM) was used.
Bączkowicz et al. [1,4] in their research presented the frequency characteristics of VAG signals by STFT. The spectra were obtained by computing Discrete Fourier Transform (DFT) of the segments (150 samples/segment), the Hanning window, and 100 samples overlap/segment. They analyzed spectral activity by summing spectral power in two different bands. The first parameter P1 concerned the range of 50-250 Hz, and the second P2 the range of 250-450 Hz. In [1] Bączkowicz et al. obtained two additional parameters by computing power of spectral density using Fast Fourier Transform (FFT): F470 at 470 Hz, and F780 at 780 Hz. They performed 2-class classification (normal/abnormal) and 5-class classification (healthy knee, the first to the third stages of chondromalacia, osteoarthritis), used the genetic search algorithm to select the best features of VAG signals, and applied four different algorithms to classify the selected features, wherein one signal feature was distinguishing all classes at once.

Enhanced Descriptors
As can be seen in Figure 1, VAG signals of different conditions consist of different frequency distortions, which could be potentially used as the features to classify them. The question arises of which frequency ranges would differentiate those classes in the best possible manner. Another question concerns the method of measuring the quality of obtained features.

Quality Measure of the Feature
Boxplot, as a visualization tool, gives very intuitive concept of the quality of a feature. It shows the variation of some numerical data, enabling the visual comparison. Particularly, it shows the median of a sample (indicated by the red line), interquartile range (indicated by the blue box), variability outside the interquartile range (indicated by the whiskers) and outliers (indicated by the red crosses). The question arises, though, how to quantitatively compare two features visualized on the boxplot? Lots of different measures are available in literature [19]. One of the most straightforward ones is the Bhattacharyya coefficient.
In a brief preliminary research, the Bhattacharyya coefficient was compared to several different coefficients (some existing in literature, like DBM/OVS [20][21][22][23] or Jaccard index [24], some newly defined ones). This comparison was conducted in the following way:

•
The best frequency ranges were generated by every coefficient.

•
Obtained frequency ranges were used to train 10 different classifiers (two decision trees, LDA, naïve Bayes, SVM, two knn classifiers, two random forests and a neural network).

•
The largest mean classification accuracies were compared.

•
The Bhattacharyya coefficient proved to be the best coefficient in this application.
It is defined as [25]: where p and q are probability distributions of two current classes and x is a domain of current feature values. The probability densities were obtained using kernel density estimation with a window: as defined in [26]. The Bhattacharyya coefficient indicates the overlap between two statistical samples. It ranges from 0 when distributions are completely separated to 1, when distributions overlap entirely. Therefore, in our study, the coefficient is minimized. The coefficient is equal to 0.401 for the g example from Figure 7 and 0.876 for the j.

Optimal Frequency Range
Finding frequency ranges in a strict brute-force manner would be too computationally expensive, since the DFT of the analyzed signal is composed of too many samples. Consequently, the search was done in three iterations, starting sparsely in the whole frequency domain, gradually narrowing the searched range and increasing precision. To visualize the results, a kind of map was generated, called in the rest of this paper the frequency range map (FRM). An exemplary map is shown in Figure 2.
where p and q are probability distributions of two current classes and x is a domain of current feature values. The probability densities were obtained using kernel density estimation with a window: as defined in [26]. The Bhattacharyya coefficient indicates the overlap between two statistical samples. It ranges from 0 when distributions are completely separated to 1, when distributions overlap entirely. Therefore, in our study, the coefficient is minimized. The coefficient is equal to 0.401 for the g example from Figure 7 and 0.876 for the j.

Optimal Frequency Range
Finding frequency ranges in a strict brute-force manner would be too computationally expensive, since the DFT of the analyzed signal is composed of too many samples. Consequently, the search was done in three iterations, starting sparsely in the whole frequency domain, gradually narrowing the searched range and increasing precision. To visualize the results, a kind of map was generated, called in the rest of this paper the frequency range map (FRM). An exemplary map is shown in Figure 2. In the FRM, the abscissa x-axis is the lower frequency range and the ordinate is the upper frequency range. As a result, in the first iteration, only the upper-left half of the map consists of possible results. Applicate, or color, indicates quality of feature, defined as the Bhattacharyya coefficient. Similar maps were defined and used in [2], but with lower resolution and for slightly different features. The z-axis, or feature evaluation coefficient used in [2] was the area under the ROC curve, not the Bhattacharyya coefficient. Consequently, the coefficient in [2] was maximized, not minimized.
On each map, the best frequency range was emphasized by a red circle. The exemplary FRM in Figure 2 shows that the best range for this (first) iteration is 114-456 Hz, in which neighboring consecutive iteration will be conducted.

Definition of the Features
Four types or "families" of features were defined. The first family is the sum of some frequency range of the discrete Fourier transform of the VAG signal: In the FRM, the abscissa x-axis is the lower frequency range and the ordinate is the upper frequency range. As a result, in the first iteration, only the upper-left half of the map consists of possible results. Applicate, or color, indicates quality of feature, defined as the Bhattacharyya coefficient. Similar maps were defined and used in [2], but with lower resolution and for slightly different features. The z-axis, or feature evaluation coefficient used in [2] was the area under the ROC curve, not the Bhattacharyya coefficient. Consequently, the coefficient in [2] was maximized, not minimized.
On each map, the best frequency range was emphasized by a red circle. The exemplary FRM in Figure 2 shows that the best range for this (first) iteration is 114-456 Hz, in which neighboring consecutive iteration will be conducted.

Definition of the Features
Four types or "families" of features were defined. The first family is the sum of some frequency range of the discrete Fourier transform of the VAG signal: d 1 is the first feature family, δf is the normalization factor, equal to the 1 It ensures that the value of the feature is affected only by the shape of the spectrum and not by its size, f 1 is the lower frequency range, f 2 is the upper frequency range, f 2 is the upper frequency range, s VAG is the VAG signal, DFT is the Discrete Fourier Transform operator, f i is the i-th frequency amplitude.
The visualization of the first feature family is provided in the (b) plot of Figure 3.   The second family of features was defined similarly, but instead of the signal, the derivative of the signal is used to obtain DFT. One of the fundamental traits of the Fourier transform is that the transform of the derivative is the transform of the original signal multiplied by the frequency: Sensors 2020, 20, 5015 where F is the Fourier transform. The derivative can be then interpreted as high-pass filter. The motivation for using the derivative was the certainty that potentially useful information, which could be embedded in wider frequency spectrum, will not be omitted. The second feature family was then defined as The derivative and this family of features were plotted on a and b of Figure 4 respectively.  The calculation of the feature values was followed by the evaluation by Bhattacharyya coefficient. The point with coordinates: x = lower frequency range, y = upper frequency range, z = the The third (Equation (6)) and the fourth (Equation (7)) families of features were defined like the first and the second, but instead of summing the DFT amplitudes, the energy was summed (so the square of the values of the previous feature families): The illustrations of those families were shown in Figures 3c and 4c respectively. Interpretation of the derivative of the signal as a high pass filter can clearly be seen in Figures 3 and 4. Higher frequencies are much more potent in Figure 4, while lower frequencies seem relatively unaltered. Energies, or the c plots of both Figures appear "sharper". It is to be expected, as the energy is the square of the spectrum.
The calculation of the feature values was followed by the evaluation by Bhattacharyya coefficient. The point with coordinates: x = lower frequency range, y = upper frequency range, z = the Bhattacharyya coefficient of the feature generated by x-y frequency range, is then plotted on FRM. Then, the next feature, composed of slightly different frequency range is evaluated, plotted, etc.

Research Methodology
For the purposes of this research, VAG signals were acquired from groups of patients qualified by specialists as having first stage chondromalacia (cmp1), second stage chondromalacia (cmp2), third stage chondromalacia (cmp3), and osteoarthritis (oa) and from the control group, healthy knee (ctrl). All subjects underwent routine medical interviews, physical examination, and imaging via MRI. Patients with stages I-III of chondromalacia patellae were classified in accordance with the Outerbridge classification [27]. In turn osteoarthritis (OA) patients were diagnosed with mild to moderate knee OA (grade II and III according to the Kellgren-Lawrence grading system [28]) with a disease duration of more than 2 years. All imaging examinations were analyzed by single radiologist, who was blinded to patients' symptoms. To prevent any signal artifacts from deteriorations other than chondral lesions, individuals with a history of knee fracture, knee surgery, history of meniscal tears, significant knee instability, or patellar maltracking were not enrolled in the study. Moreover, due to the methodology of the VAG assessment, individuals with a restricted knee joint range of motion (required 0 • to 100 • ), significantly weakened muscles, and substantially swollen knees in the affected lower limb were excluded from the study. The VAG examination was performed in a sitting position and consisted of exactly four full cycles of alternating extension and flexion of the knee joint (90 • -0 • -90 • ) lasting 6 s. Mounted 1 cm above the apex of patella, the acceleration sensor was recording vibration and acoustic processes occurring during the knee movements. The acceleration sensor attachment was shown in Figure 5. Bhattacharyya coefficient of the feature generated by x-y frequency range, is then plotted on FRM. Then, the next feature, composed of slightly different frequency range is evaluated, plotted, etc.

Research Methodology
For the purposes of this research, VAG signals were acquired from groups of patients qualified by specialists as having first stage chondromalacia (cmp1), second stage chondromalacia (cmp2), third stage chondromalacia (cmp3), and osteoarthritis (oa) and from the control group, healthy knee (ctrl). All subjects underwent routine medical interviews, physical examination, and imaging via MRI. Patients with stages I-III of chondromalacia patellae were classified in accordance with the Outerbridge classification [27]. In turn osteoarthritis (OA) patients were diagnosed with mild to moderate knee OA (grade II and III according to the Kellgren-Lawrence grading system [28]) with a disease duration of more than 2 years. All imaging examinations were analyzed by single radiologist, who was blinded to patients' symptoms. To prevent any signal artifacts from deteriorations other than chondral lesions, individuals with a history of knee fracture, knee surgery, history of meniscal tears, significant knee instability, or patellar maltracking were not enrolled in the study. Moreover, due to the methodology of the VAG assessment, individuals with a restricted knee joint range of motion (required 0° to 100°), significantly weakened muscles, and substantially swollen knees in the affected lower limb were excluded from the study. The VAG examination was performed in a sitting position and consisted of exactly four full cycles of alternating extension and flexion of the knee joint (90°-0°-90°) lasting 6 s. Mounted 1 cm above the apex of patella, the acceleration sensor was recording vibration and acoustic processes occurring during the knee movements. The acceleration sensor attachment was shown in Figure 5. After acquiring the signals, their spectra were obtained, and four feature families were proposed. For each family, a set of features containing different frequency ranges was generated, and every feature from this set was evaluated by the Bhattacharyya coefficient. The best features constructed final 10-element vector. This vector was compared to the state-of-the-art feature vector [1], by using both of those vectors as an input to some classification algorithms. The accuracies of trained algorithms determined which feature vector is superior.

Acquisition of the VAG Signal
The signals analyzed in this paper were obtained from knee joints of patients diagnosed by physicians into five groups: control group (containing 66 signals), I stage chondromalacia patellae (26 signals), II stage chondromalacia patellae 30 signals), III stage chondromalacia patellae (36 signals), and osteoarthritis (26 signals). The diagnoses were based on the X-ray. After acquiring the signals, their spectra were obtained, and four feature families were proposed. For each family, a set of features containing different frequency ranges was generated, and every feature from this set was evaluated by the Bhattacharyya coefficient. The best features constructed final 10-element vector. This vector was compared to the state-of-the-art feature vector [1], by using both of those vectors as an input to some classification algorithms. The accuracies of trained algorithms determined which feature vector is superior.

Acquisition of the VAG Signal
The signals analyzed in this paper were obtained from knee joints of patients diagnosed by physicians into five groups: control group (containing 66 signals), I stage chondromalacia patellae (26 signals), II stage chondromalacia patellae 30 signals), III stage chondromalacia patellae (36 signals), and osteoarthritis (26 signals). The diagnoses were based on the X-ray.

The Process of Defining New Frequency Features
For every sample of 184 VAG signals, the frequency analysis was conducted, resulting in four spectral vectors: • the DFT of the VAG signal, from which the first feature family is obtained, using Equation (3), • the DFT of the derivative of the VAG signal, from which the second feature family is obtained, using Equation (5), • squared DFT of the VAG signal, from which the third feature family is obtained, using Equation (6), • squared DFT of the derivative of the VAG signal, from which the fourth feature family is obtained, using Equation (7).
To generate FRMs, Bhattacharyya coefficient was calculated for the sums of desired frequency ranges. As already mentioned, evaluating all possible frequency ranges would be too computationally expensive, therefore the search was done in three iterations. The search started from the full frequency range (0-5 kHz) with quite big frequency step and gradually narrowing the range and reducing the step. The last iteration was done with highest precision (the smallest frequency step) in quite narrow frequency range. Specifically, three iterations were conducted as follows: 1.
the first one in range 0-5 kHz (since signal was sampled with the frequency of 10 kHz), with about 10 Hz step. Features were then defined as sums, defined by Equations (3) and (5) the second iteration was conducted in range ±800 Hz from the best range obtained in previous iteration with about 2.5 Hz step, 3.
the last iteration was conducted in range ±80 Hz from the best range obtained in previous iteration with about 0.16 Hz step.
Three exemplary maps were shown on Figure 6, illustrating three consecutive iterations.
Sensors 2020, 20, 5015 8 of 14 For every sample of 184 VAG signals, the frequency analysis was conducted, resulting in four spectral vectors:  the DFT of the VAG signal, from which the first feature family is obtained, using Equation (3),  the DFT of the derivative of the VAG signal, from which the second feature family is obtained, using Equation (5),  squared DFT of the VAG signal, from which the third feature family is obtained, using Equation (6),  squared DFT of the derivative of the VAG signal, from which the fourth feature family is obtained, using Equation (7).
To generate FRMs, Bhattacharyya coefficient was calculated for the sums of desired frequency ranges. As already mentioned, evaluating all possible frequency ranges would be too computationally expensive, therefore the search was done in three iterations. The search started from the full frequency range (0-5 kHz) with quite big frequency step and gradually narrowing the range and reducing the step. The last iteration was done with highest precision (the smallest frequency step) in quite narrow frequency range. Specifically, three iterations were conducted as follows: 1. the first one in range 0-5 kHz (since signal was sampled with the frequency of 10 kHz), with about 10 Hz step. Features were then defined as sums, defined by Equations (3) and (5)-(7), in subsequent ranges: 0-10 Hz, 0-20 Hz, …, 0-5000 Hz, 10-20 Hz, 10-30 Hz, …, 4980-4990 Hz, 4980-5000 Hz, 4990-5000 Hz. Every range was evaluated by Bhattacharyya coefficient and plotted as a point on the FRM. 2. the second iteration was conducted in range ±800 Hz from the best range obtained in previous iteration with about 2.5 Hz step, 3. the last iteration was conducted in range ±80 Hz from the best range obtained in previous iteration with about 0.16 Hz step.
Three exemplary maps were shown on Figure 6, illustrating three consecutive iterations. The FRMs were separately generated for every class combination (details in Tables 1-3), for every feature family, so 120 maps in total, 40 maps for final iteration. Generating separate maps for each class combination may seem redundant. Only distinction between neighboring classes appears to be sufficient. Different approach, i.e., temporary ensemble all classes but one, generating maps only for distinguishing between, for example, the second stage chondromalacia condition and the rest, generating characteristic features for each class also seems enough. As it certainly would be less computationally expensive, it probably would be less robust. Neighboring classes are not completely separate and independent, as the first stage chondromalacia precedes the second stage etc. Consequently, it would be quite difficult to find the frequency range typical for particular class only. Minimizing Bhattacharyya coefficient for every class pair ensures two things: The FRMs were separately generated for every class combination (details in Tables 1-3), for every feature family, so 120 maps in total, 40 maps for final iteration. Generating separate maps for each class combination may seem redundant. Only distinction between neighboring classes appears to be sufficient. Different approach, i.e., temporary ensemble all classes but one, generating maps only for distinguishing between, for example, the second stage chondromalacia condition and the rest, generating characteristic features for each class also seems enough. As it certainly would be less computationally expensive, it probably would be less robust. Neighboring classes are not completely separate and independent, as the first stage chondromalacia precedes the second stage etc. Consequently, it would be quite difficult to find the frequency range typical for particular class only. Minimizing Bhattacharyya coefficient for every class pair ensures two things: • that the found frequency ranges for each classes are as optimal as possible, without sacrificing the distinction between more different conditions (maybe one frequency range would be appropriate for the distinction between the first and the second stage chondromalacia, but not as effective with distinguishing between first stage chondromalacia and the osteoarthritis; the question would arise, which distinction should be dominant, and which should be sacrificed), • that the obtained frequency ranges provide as unambiguous distinction as possible; the situation can be imagined in which two classes neighboring one another (for example cmp1 and cmp3 neighboring cmp2) can be distinguished from the one with very similar frequency ranges. Then the statement that the particular frequency range is typical for, e.g., cmp2 would be true, but the contrary would not unambiguously point out cmp1 or cmp3.
Analysis of the final 40 maps led to definition of 10 frequency ranges in different feature families, indicated by bolded font in the Table 2. Obtained features were compared to standard frequency descriptors defined in [1].

The Verification of Defined Features as a Classifier Input
The comparison between newly defined features and the state-of-the-art was done in two ways. First, the Bhattacharyya coefficients of features were compared. However, as mentioned before, a 4-element (P1, P2, F470, F780) spectral feature vector was proposed in [1]. In this study, 10 features were proposed, making it hard to directly compare the quality of obtained features.
To compare obtained features, both newly defined and state-of-the-art were used to train a couple of classifier algorithms. Utilized algorithms were listed in Table 4. To find which feature vector is superior, the final classification accuracies were compared.
The classification accuracy of every classifier strongly depends on the division of the data into the teaching, testing, and validation sets. To surpass that potential bias, 1024 classifiers of each type were constructed with random division of the data into mentioned sets.

Results and Discussion
The results obtained for every combination of classes, for each feature family, were presented in Tables 1 and 2. The use of the derivative was quite fruitful for some class pairs, even though intuition behind it, i.e., the informativeness of higher frequency spectrum, was not fully correct. For example, the distinctions between the first and the second stage chondromalacia or the third stage chondromalacia and osteoarthritis were better with the use of the derivative, but the range of the frequencies was very narrow. Generally, the most optimal results obtained for the primitive VAG signal and its derivative are comparable.
Boxplots of the best results for each combination of classes were shown in Figure 7.
The new features were defined with the emphasis on distinction between each particular class pair. Then, it can be seen, that some of the pairs are well distinguished, especially when the knee joint conditions strongly differ (for example, the control group and the third stage chondromalacia or the first stage chondromalacia and the osteoarthritis). The problem persists, though, with the distinction between neighboring pairs, especially the third stage chondromalacia and the osteoarthritis. This is to be expected to some extent, as the chondromalacia condition commonly occurs with the osteoarthritis.
To better evaluate the quality of the obtained features, they were compared to the features defined in [1] (used signal database was the same in this research and in [1]). Firstly, the Bhattacharyya coefficient was generated for the frequency features from [1], then some classification algorithms were used to quickly check potential classification accuracies, using both the new features and the features defined in [1].
The comparison between features defined in [1] and features defined in this paper was presented in Table 3. The boxplots of the features from [1] were given in Figure 8.
None of the newly obtained frequency ranges was close to the range of P1 or P2 from [1]. The frequencies around 250 Hz, the interface of P1 and P2, occur in the Table 2 relatively often, however the ranges are significantly narrower. The new feature closest to F470 from [1] is feature e from Table 2, so the sum of the DFT of the derivative of the signal for the frequencies around 400 Hz. Nonetheless, the difference in frequencies is quite big. The F780 feature from [1] does not overlap with any features obtained in conducted research.
The new features are superior, in terms of Bhattacharyya coefficient, for every class combination. Improvement of distinction is quite significant for some pairs (like the pair cmp1-cmp3, for which the coefficient dropped by more than 40%), and for others is smaller (line the cmp3-oa pair, for which the coefficient dropped by around 5%).
To test the features as a diagnostic tool, ten different types of classifiers were constructed for the newly defined features and another ten classifiers for the features defined in [1]. The comparison between their accuracies is presented in Table 4.   Boxplots of the best results for each combination of classes were shown in Figure 7. The new features were defined with the emphasis on distinction between each particular class pair. Then, it can be seen, that some of the pairs are well distinguished, especially when the knee joint conditions strongly differ (for example, the control group and the third stage chondromalacia or the first stage chondromalacia and the osteoarthritis). The problem persists, though, with the distinction between neighboring pairs, especially the third stage chondromalacia and the osteoarthritis. This is to be expected to some extent, as the chondromalacia condition commonly occurs with the osteoarthritis.
To better evaluate the quality of the obtained features, they were compared to the features defined in [1] (used signal database was the same in this research and in [1]). Firstly, the Bhattacharyya coefficient was generated for the frequency features from [1], then some classification algorithms were used to quickly check potential classification accuracies, using both the new features and the features defined in [1].   None of the newly obtained frequency ranges was close to the range of P1 or P2 from [1]. The frequencies around 250 Hz, the interface of P1 and P2, occur in the Table 2 relatively often, however the ranges are significantly narrower. The new feature closest to F470 from [1] is feature e from Table  2, so the sum of the DFT of the derivative of the signal for the frequencies around 400 Hz. Nonetheless, the difference in frequencies is quite big. The F780 feature from [1] does not overlap with any features obtained in conducted research.
The new features are superior, in terms of Bhattacharyya coefficient, for every class combination. Improvement of distinction is quite significant for some pairs (like the pair cmp1-cmp3, for which the coefficient dropped by more than 40%), and for others is smaller (line the cmp3-oa pair, for which the coefficient dropped by around 5%).
To test the features as a diagnostic tool, ten different types of classifiers were constructed for the newly defined features and another ten classifiers for the features defined in [1]. The comparison between their accuracies is presented in Table 4.  All accuracies presented in the Table 4 are actually the mean values of 1024 classifiers of a given type.
As expected, classifiers constructed on new feature set proved to be more accurate. Additional visual comparison was given with the use of boxplot on the Figure 9. All the data came from the Table 4.
Although the new features are more useful for classifying different knee joint conditions, obtaining them is more computationally expensive. Additionally, they are probably somewhat correlated to each other, which could be undesirable while using them as an input for some classification algorithms [29]. A principal component analysis could be a useful step before constructing classifier, but it would increase computational cost even further.
All accuracies presented in the Table 4 are actually the mean values of 1024 classifiers of a given type.
As expected, classifiers constructed on new feature set proved to be more accurate. Additional visual comparison was given with the use of boxplot on the Figure 9. All the data came from the Table 4.  Table 4.
Although the new features are more useful for classifying different knee joint conditions, obtaining them is more computationally expensive. Additionally, they are probably somewhat correlated to each other, which could be undesirable while using them as an input for some classification algorithms [29]. A principal component analysis could be a useful step before constructing classifier, but it would increase computational cost even further.  Table 4.

Conclusions
Conducted research resulted in the definition of ten new spectral features, emphasizing differences between every possible knee condition class pair. This ensured that the found frequency ranges were as optimal as possible, without the risk of sacrificing potentially valuable information about the differences between some class pairs. Proposed features were used to train ten different classification algorithms proving their superiority as signal descriptors. In comparison to [1], in terms of Bhattacharyya coefficient, newly proposed features proved to be over 25% better. The classification accuracy has increased on average by 9%.
Frequency Range Maps were proposed as a spectral feature visualization tool, utilizing Bhattacharyya coefficient as a feature quality measure. Visual analysis of those maps may help better understand the nature of the signal. The maps are not limited to VAG signals and can be used in different fields of research.
Proposed features may constitute the reference point in future studies, utilizing proposed frequency ranges as a method of measuring effect of some factor on quality of knee joint itself.