{'fraction_non_alphanumeric': 0.029678201872417844, 'fraction_numerical': 0.015502406852897698, 'mean_word_length': 4.461291658041813, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-Frequency Analysis and Parameterisation of Knee Sounds for Non-invasive Detection of Osteoarthritis

Objective: In this work the potential of non-invasive detection of knee osteoarthritis is investigated using the sounds generated by the knee joint during walking. Methods: The information contained in the time-frequency domain of these signals and its compressed representations is exploited and their discriminant properties are studied. Their efficacy for the task of normal vs abnormal signal classification is evaluated using a comprehensive experimental framework. Based on this, the impact of the feature extraction parameters on the classification performance is investigated using Classification and Regression Trees (CART), Linear Discriminant Analysis (LDA) and Support Vector Machine (SVM) classifiers. Results: It is shown that classification is successful with an area under the Receiver Operating Characteristic (ROC) curve of 0.92. Conclusion: The analysis indicates improvements in classification performance when using non-uniform frequency scaling and identifies specific frequency bands that contain discriminative features. Significance: Contrary to other studies that focus on sit-to-stand movements and knee flexion/extension, this study used knee sounds obtained during walking. The analysis of such signals leads to non-invasive detection of knee osteoarthritis with high accuracy and could potentially extend the range of available tools for the assessment of the disease as a more practical and cost effective method without requiring clinical setups.


I. INTRODUCTION
O STEOARTHRITIS is the most common disabling and financially burdensome of all musculoskeletal diseases, and prevalence is rising. It occurs most frequently in the knee, affecting 1 in 5 adults over the age of 45 [1]. It leads to pain, stiffness and swelling of the joint, greatly degrading the quality of life. Risk of Osteoarthritis (OA) is associated with increased mechanical wear, such as through older age and high body weight [2]. Currently, there is no cure and treatments aim to manage symptoms through lifestyle modification, physioand pharmacological therapy [1]. In severe cases, total knee replacement is required.
Clinical detection of knee OA relies on a combination of patient reported symptoms and medical imaging of cartilage This work has been submitted to the IEEE for possible publication. Copyright c 2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
C. Yiallourides and P. A. Naylor are with the Department of Electrical and Electronic Engineering, Imperial College London, UK. costas.yiallourides08@imperial.ac.uk, p.naylor@imperial.ac.uk and subchondral bone degradation. Current imaging methods such as X-ray, Magnetic Resonance Imaging (MRI) and ultrasound have poor sensitivity in early disease and as a result, at the time of diagnosis, OA is already at a progressed stage, and understanding of its cause and development is still limited. Additionally, current imaging techniques provide images of the static anatomical structure of knee joints at a particular posture and are therefore limited in assessing the dynamic integrity of the knee during a dynamic Open Chain Activity (OCA), for example when the foot leaves and makes contact again with the ground as happens during walking. This is important since OA patients experience pain and discomfort when their knee is functional. Although dynamic MRI produces good measurements in the assessment of knee function, it is normally not practical in terms of cost and accessibility [3]. Hence, there exists the need for a quick, noninvasive, portable and cheaper technique that would ideally be accessible in a non-clinical environment and could be used as a screening tool for the mild disease cases.
Joints generate sounds during movement. When the knee is active, the joint between the tibia and the femur bones moves. The regions and perhaps the quality of joint surfaces coming into contact are different at each angular position, generating therefore a number of different sounds during movement. In healthy knee joints, the bones have smooth surfaces due to a thin layer of cartilage and are separated by a protective space filled with synovial fluid to reduce friction [4]. They are able to move freely and the level of sound emitted is low. In OA knees this structure is degraded and the protective space and associated lubrication reduce, resulting in increased friction which accelerates the wear of cartilage [2]. This increased friction makes the knee more noisy during motion.
The potential for using knee joint sounds for diagnostic purposes has been known for many years. Blodgett, in 1902, reported on auscultation of the knee, with attention to sounds of normal joints and their change with repetitive motion, where a relation between an increase of sound activity and age was noted [5]. In 1913, Bircher reported that different types of meniscal injury generate distinctive sound signals [6], [7]. Steindler in 1937, used a system consisting of a cardiophone, an oscilloscope and a recorder to study 397 knees [8]. He found a relation between pathologies and the pitch, amplitude and the sequence of sounds and was able to classify the joints based on these features. It was observed, however, that it was difficult to separate other body sounds such as muscle activity from articular cartilage sounds. In [9] the authors claim that arXiv:2004.12745v1 [eess.AS] 27 Apr 2020 sounds could be detected in rheumatoid arthritis before any changes were observable in an X-Ray image but no further work was conducted to confirm this claim.
Auscultation based Phonoarthrography (PAG) utilises microphones in the audible frequency range to record sounds generated during movement. Important work on PAG by Chu et al. reported that the spectral activity of pathological knees (recorded during active motion) spanned the entire audible frequency range and the signals' acoustic power increased with severity of cartilage damage [10], [11], [12], [13]. Significant work was directed to the development of Vibroarthrography (VAG) as an alternative to PAG which relies on accelerometer sensors, operating at frequencies below 1 kHz, to pick up mechanical vibrations. Algorithms proposed for classifying the knee VAG signals according to pathological conditions, range from linear prediction modelling [14], [15] to timefrequency analysis [16], [17], [18] and wavelet matching pursuit decomposition [19], [20]. Several features have been used for classification, including spectrogram features, waveform variability parameters, statistical features [21], fundamental frequency, mean amplitude of pitch and their jitter and shimmer [22], [23]. Classifiers used in the literature range from early neural network architectures [18], [24] to maximal posterior probability decision criterion [25], bagging ensemble and multiple classifier system based on adaptive weighted fusion [19]. A thorough description of VAG analysis can be found in [26].
The use of Acoustic Emission (AE) at ultrasonic frequencies was explored as a potential biomarker for assessing the knee joint condition. In [27] piezoelectric contact sensors were used to capture ultrasonic AE signals (50 kHz to 200 kHz) emitted during sit-to-stand movements and was demonstrated, using Principal Component Analysis, that healthy and OA knees are separable in the feature space. It was further concluded that OA knees produce substantially more AE events with higher peak magnitude and average signal level [28], [29], [30].
AE analysis during knee flexion-extension was also explored in the context of knee injury rehabilitation [31], [32], [33]. In [32], a 64-dimensional feature representation of 200 ms frames of the knee sound signal was used, from which a k-Nearest Neighbors (kNN) graph was constructed. A graph based metric was then proposed to quantify the homogeneity of the feature matrix without modelling the underlying distribution. Based on this metric, it was concluded that injured knee joints produce more heterogeneous features than healthy knee joints [32]. Although this approach alleviates the need for prior algorithm training, it is only accurate when sound data from both knees of an injured subject is available since the study focused on the intra-subject knee sound differences. Inter-subject sound differences were not considered. This is one of the most challenging aspects of knee joint sound analysis as there is strong variability in the knee sounds amongst individuals which is likely due to their joints' structural differences [33].
The movement protocols most often reported in the literature, for OA and other arthritis related studies, are knee flexion-extension and sit-to-stand movements [26], [27], [34]. In the work presented here, acoustic signals are captured by a contact microphone attached to the patella, while patients are walking on a specialised treadmill. Various feature-based descriptions for these signals are investigated. In particular, discriminative features are sought that are relevant to the analysis and classification of normal (clinically healthy) and abnormal (clinically OA) knee joints. A preliminary version of this work was presented in [35]. Here this work is extended to consider time-frequency representations of the knee acoustic signals as features and examine their discriminatory power. In addition, a study is presented on the effect on the classification performance of the choice of parameters in the feature extraction step.
The main aim of this work is to answer four questions: (a) Does the classification performance improve when the Discrete Fourier Transform (DFT) spectrum is compressed using triangular filter-banks? (b) Does the classification performance improve when the natural logarithm and Discrete Cosine Transform (DCT) are used instead of only the DFT? (c) Is the classification performance better when using uniform or non-uniform frequency spacing in the analysis? (d) Which frequency ranges of AE signals contain more discriminative information and hence are important for OA classification? In answering these questions, an insight will be obtained into which features best characterise normal and OA knees. Several classifiers are used but their optimisation is beyond the scope of this paper.
The remainder of the paper is structured as follows. Section II describes the features considered in this work and introduces the notation that is used throughout the paper. This is followed by feature analysis and selection in Section III. Experiments along with results and discussion are presented in Section IV where detailed information about the data acquisition and assessment protocol can be found in IV-B. Finally, Section V concludes the paper with a summary of the proposed work.
II. FEATURE EXTRACTION Acoustic signals are recorded over the patella using a contact microphone, as will be described in Section IV-B. Let s i (n) denote the signal at discrete time index n captured by the patella microphone for the i th knee in the data-set, where i = 1, 2, . . . , I for I knees. Prior to extracting features, all recorded signals are normalised to have equal Root Mean Square (RMS) level. We hypothesize that the acoustic artifacts caused by walking on the treadmill are uncorrelated with the features used for the analysis. Furthermore, it is assumed that sounds related to abnormalities appear within time periods of τ s seconds. Accordingly, s i (n) is divided into non-overlapping segments of length τ s , denoted as s i,j (n) for j = 1, 2, . . . , J i segments. Each segment is then labelled for classification according to the condition of the knee from which it was obtained.
A signal segment s i,j (n) is further divided into frames of length l ms with 50% overlap. This creates an N f × l n matrix S where N f is the number of frames and l n is the frame length in samples. Considering a hanning window of length l n transformed into the diagonal matrix H, the DFT of S can be computed as where for each row n = 0, 1, . . . , l n − 1 and for each column k = 0, 1, . . . , l n − 1 where k is the frequency index. By taking the magnitude of each element in Ψ f and retaining only the first K = floor(1 + l n /2) columns, the matrix Ψ F is constructed.
A filter-bank with N B triangular band-pass filters linearly spaced in frequency is used to construct the matrix where each element is the magnitude of the bandwidth of a single filter at a single frequency bin k = 0, 1, . . . , K − 1.
A matrix U M is similarly constructed from triangular filters that are equally spaced along the mel-frequency axis which is defined as in [36]. A compact spectrum representation can then be obtained as The resultant matrices are N f × N B and in this way dimensionality reduction is achieved. The columns of Ψ E , Ψ D and Ψ F can be considered as distributions of the spectrum magnitude values for particular frequency bands for the former two and bins for the latter. From each such distribution a feature vector f is extracted using 11 statistical parameters that capture certain signal attributes that aim to highlight differences between OA and healthy signals. These parameters are chosen in this work to be the mean, kurtosis, variance, skewness, max, min and the 10 th , 25 th , 50 th , 75 th , 90 th percentiles. The vector f can be obtained by defining an operator D(·) which acts on a matrix and returns the 11 statistical parameters of each column. Therefore, Extraction procedures inspired by human auditory perception are widely used in many applications. Mel-frequency Cepstral Coefficients (MFCC) are a common choice of features that are successfully used in speech recognition and music genre classification. MFCC have been previously used for the analysis of VAG signals, [37], [38], but have not so far been used for OA detection from the analysis of acoustic signals emitted from the knee and sensed at the patella. Their extraction process involves mapping the power of the Short Time Fourier Transform (STFT) spectrum using triangular overlapping windows onto the mel scale which is designed to approximate the human auditory system's response. The aim is to exploit the property of the mel scale and apply it to the knee signals. In particular, for the sounds heard as pops, clicks, grindings etc. during knee motion. The use of the mel scale for knee signals is motivated by the fact that these sounds can be distinguished in a recording by even the untrained ear with minimal effort and that they are sounds that are likely generated by the friction between the tibia and the femur bones which in turn is caused by the effects of OA in the knee.
MFCC are computed according to [39] and in the same way, but replacing the mel-frequency filter-bank with linearly spaced filters, a set of Linear-frequency Cepstral Coefficients (LFCC) is also computed. Hence denote the two sets of Cepstral Coefficients (CC), where C(·) is the DCT operator. Applying the same thinking as with the STFT based representations, the D(·) operator is employed to obtain the statistical representation of the matrices in (5) as Given the fixed frame segmentation process employed using short time-frames, it is likely that a knee sound related to OA might extend to more than one frame. By taking the time derivatives of the coefficients, the information present in the evolution of these sounds across a multiple of frames can be extracted. Therefore, for each element in Ψ D , Ψ E , C M and C L , defined as a static coefficient and denoted as a t , the first derivatives in time are computed using where U = 4 and d t is termed a delta coefficient from frame t computed using the static coefficient of that frame. The second derivatives, termed as delta-delta coefficients, are also computed using (7) but replacing a t with d t and setting U = 1 which makes it a simple difference equation. The choice of U = 1 and U = 4 is adopted from speech recognition.
All the values of the first derivative obtained from a single feature in either C M , C L , Ψ E and Ψ D , (i.e. cepstral coefficient or frequency band), can be considered to form a distribution from which the 11-dimensional vector f is extracted. The same is performed for all the features in these 4 matrices for their first and second derivatives and the new vectors are appended to the appropriate φ feature matrix. Performing the above process for all signal segments produces 5 feature sets that will be used for the classification experiments. These are for c = I i=1 J i total segments. In the following Section, the effectiveness of the signal segment parameterisation using the 11 statistics is examined and their discriminant power for the task of normal vs abnormal segment classification is studied.

III. FEATURE ANALYSIS AND SELECTION
Let x refer to any of the symbols {D, E, F, L, M }. The aim is to evaluate the discriminant power of each The classifier employed for this purpose is Support Vector Machine (SVM) because it is efficient for small training data and avoids making any assumptions on the underlying data distribution [40]. This makes it a suitable choice since the distribution is unknown. The linear kernel SVM is used, denoted as SVM l , as it is the simplest form of kernel and is less prone to overfitting than other more complex kernels [40].
The knee database, obtained as described in Section IV-B and used in the experiments to follow, is comprised of 19 normal knees and 21 abnormal from which the signal segments s i,j (n) are obtained using a window size of τ s seconds.
A cross-validation procedure is employed using 5 groups, randomly constructed from the database, with a normal to abnormal knees ratio of 3:5, 3:5, 3:5, 5:3, 5:3 for each group which are then made up with the segments of their constituting knees. In this way the problem of having segments of a knee signal in more than one group is avoided. Some variability in the group sizes exists given that some knee recordings in the database are longer than others and therefore have more segments. Four groups are used for training the SVM model which is then tested on the group left out. This is repeated until all 5 groups are evaluated. Prior to this, the training data (4 groups) is scaled by subtracting the mean and normalising by the variance. The same scale values are then applied to the test data (the group left out). The above procedure is executed 100 times in order to reduce the variance of the estimator and the results are averaged at the end.
The performance of each feature is assessed based on several metrics. Relying only on the error rate (E r ) is often not sufficient to draw safe conclusions since potential classification errors other than the number of misclassified observations are not captured by the error rate. Hence, the F 0.5 measure (a variation of F 1 ) and Matthew's Correlation Coefficient (MCC) are also used. Both are computed from the confusion matrix [41]. From a clinical perspective, the false prediction of abnormal segments as normal is worse than the contrary. F 0.5 emphasises this error type more than F 1 and is thus preferred. MCC is a balanced measure ranging from -1 (prediction totally different from observation), to 1 (perfect prediction), with 0 stating no better than random prediction [42].
Following the pre-processing with feature extraction and analysis steps, the selection of feature subsets for subsequent classification and analysis is performed. The selection method used is a hybrid of a filter and a wrapper approach. First, the features are ranked in each of the 3 metric categories (the filter step) and with the application of thresholds, θ er , θ 0.5 , θ mcc for E r , F 0.5 , MCC respectively, the best N features are selected, where N is dictated by those that satisfy all thresholds. Secondly, by allowing θ er , θ 0.5 and θ mcc to vary in the range 0, 1 with discrete steps of size w, the entire feature space is searched and all possible subsets are constructed subject to these constraints (the wrapper step). The three metrics are bounded in a continuous range and therefore defining a discrete set of constraints is necessary in order to make the search space tractable as it is not practical to test all possible combinations of features as an exhaustive search. This feature selection method forms nested subset of features, q is an index that takes integer values in the range 1 to L where L = N B + 2K + 3 and J(.) is any of E r , F 0.5 , MCC , evaluated against the corresponding threshold. Each subset is used for training and testing the SVM classifier by employing the cross-validation procedure described earlier. Their classification performance is evaluated using the area under the curve (AUC) of the Receiver Operating Characteristic (ROC) curve [43]. The subset that gives the highest AUC is chosen.
The experimental framework is based on a systematic approach that aims to (a) find the best frame length l for extracting the 5 alternative feature sets in (8), (b) examine the effect of the number of filters N B on the classification performance when using either one of Φ D , Φ E , Φ L or Φ M and (c) obtain insights into the time-frequency information of normal and abnormal signals and their differences.

IV. EXPERIMENTS, RESULTS AND DISCUSSION
An investigation of the effect of frame length values based on a deterministic approach was initially conducted by defining a suitable range and quantifying the classifier performance in order to choose the best l. Subsequent, Monte Carlo simulations were performed to test the suitability of the choice and to identify performance trends of the feature sets in a larger range. Finally, experiments were conducted varying the number of filters N B .

A. Implementation details
In all the experiments that follow, the signal segment length τ s was set to 20 s. Other time periods that do not violate the assumption outlined at the beginning of Section II could also be used but would affect the total number of segments obtained. Additionally, the error rate threshold θ er was fixed at 0.456 which is the error rate obtained when the predicted class is always the largest. This is the error rate attributed to random guessing and hence anything worse than this would mean that the classifier performs poorly. By keeping θ er constant, the values of θ 0.5 and θ mcc are varied in the range [0, 1] with stepsize w = 0.05 and the possible feature subsets are constructed and subsequently used in the classifier.
The L 1 -norm soft margin formulation is used for SVM due to its advantages over the L 2 -norm in high dimensional feature spaces and in the presence of redundant features [44]. The sequential minimal optimization, [45], is employed for solving the SVM minimization problem which is the standard algorithm for this task. The penalty parameter, often called box constraint, is a term that trades off misclassification of training observations against simplicity of the decision surface. A low value makes the surface smooth (i.e. misclassification becomes less important), while a high value attempts to classify all training examples correctly. In the following experiments this parameter was set to 1. Finally, the formula used for the Gaussian kernel is exp(−γ x 1 − x 2 2 ) for which γ = 1 and x 1 , x 2 denotes the inner product between the training vectors x 1 and x 2 . Linear Discriminant Analysis (LDA), [46], is also employed in certain experiments and the empirical pooled covariance matrix is used for the multivariate normal distribution of each class. Finally, for the Classification and Regression Tree (CART) classifier, the split predictor (feature) is selected as the one that maximizes the split criterion gain (gini index) over all possible splits of all predictors [47]. The tree, once fully grown, is pruned using the gini index as the pruning criterion.

B. Data Acquisition and Test Protocol
Adults with clinical knee OA and reporting no knee pain in the last 2 weeks were recruited. Knees were classified by clinicians as: 1) normal (clinically healthy), or 2) abnormal (OA). Exclusion criteria were: aged <18 years, previous surgery, unable to provide consent. AE signals were acquired during walking on a treadmill instrumented with force plates. The recordings were made with a sampling frequency of 44.1 or 48 kHz and downsampled to F s = 16 kHz for subsequent processing, using a contact microphone with a sound port for detecting airborne sounds and an electret condenser microphone mounted inside a capsule (Basik Pro Schertler, 20 Hz -20 kHz), attached over the patella and supported by a digital preamplifier (RME Babyface; PreSonus DigiMax LT).
The assessment commenced with a 5 minute warm-up and acclimatisation to treadmill walking followed by data acquisition at increasing speeds on a flat level until maximum walking speed was achieved (speed increments of 0.5 km/h, maximum walking speed defined as the maximum pain-free speed where one foot was always in contact with the ground). Maximum speeds achieved ranged from 2.5 to 9 km/h. Data used in this work originates from 40 knees, of which 19 are normal (from 15 patients) and 21 are abnormal (from 18 patients). Table I displays the demographical characteristics of the study participants. Approximately 83 minutes of sound data from healthy knees and 99 minutes from OA is used. Following the segmentation process described in Section II, 249 normal and 297 abnormal segments of 20 s are obtained. C. Deterministic search in a specified range of frame lengths The frame length l is tested for the values 20, 24, 28, . . . , 100 ms. The limiting values were chosen so that 20 ms is a short enough window to allow good time resolution in the time-frequency representations for the sounds (clicks, pops, grindings) heard during walking and 100 ms is a large enough window to capture the two major events in a single stride, namely, the heel strike and the push off responses as captured by the patella microphone. This information was extracted from the ground reaction force signals obtained from the treadmill's force plates from which can be extracted the timings of each event in the gait cycle.
The experimental framework developed in the previous Sections is applied and, for each frame length l, the features are extracted, analysed independently and the possible subsets are constructed for evaluation using the SVM l classifier. In all cases N B = 20 is chosen, giving 60 coefficients (including the 0 th cepstral coefficient) for all the feature sets except Φ F . For Φ F , the size depends on l since the DFT length used is equal to the frame length, as shown in Section II.
For the analysis of the results, a top-down approach is followed. The overall results are summarised in Fig. 2. Each point on a line represents the highest AUC, averaged over 100 trials, obtained by any subset of the corresponding frame-size and of the particular feature set.
It is evident from Fig. 2 that classifying using features obtained from the Φ M feature set, including first and second derivatives, scores consistently higher than any other set, for any l, with its best performance occurring with l = 48 ms (AUC = 0.915). On the other hand, the sets Φ F , Φ L , Φ E and Φ D peak at 24 ms, 20 ms, 88 ms and 20 ms respectively. Initial observations suggest that the optimal framesize is different for each feature set and that reducing the dimensionality of the spectrum with triangular shaped filters generates comparable results with the full spectrum features and in some cases (l ≥ 72 ms) improves the classification performance.
1) Comparison of MFCC and LFCC results: It is seen in Fig. 2 that AUC values range from 0.878 to 0.915 and from 0.77 to 0.84 for Φ M and Φ L respectively. Hence, the performance for Φ M is only weakly sensitive to the choice of frame-size and in fact, above 32 ms the variance of the metric value drops approximately to 1 /3 of the variance obtained when including the AUC for l < 32 ms. The fluctuation is small because, for any l, the final selected subset that gave the highest AUC consisted of only f M 3 which contains the statistical parameters of the distribution of the 2 nd MFCC.   All other subsets considered from this set resulted in worse performances. For l ≥ 88 ms and l ≤ 32 ms, the AUC of Φ M follows a downward trend suggesting that, even to a small degree, the performance drops for frame-sizes outside of this range. This hypothesis is tested in Section IV-E. Higher variability is observed in the results of Φ L where more than one feature vector is chosen in the final subset for 10 out of 21 frame lengths. The highest AUC for this set occurs at 20 ms which is obtained with a subset of 17 feature vectors. The classification rate for this case, averaged over the 100 crossvalidation iterations, is 75.07% compared to 85.25% obtained from the Φ M feature set with l = 48 ms. This translates into classifying correctly, 55 more segments (from a total of 546 segments) when using the Φ M set. Fig. 3 shows a box plot of the differences in AUC for all frame lengths per Cepstral Coefficients (CC). It is evident that many CC obtained with a linear-frequency filter-bank generated better classification results than their mel-frequency filter-bank counterparts. More precisely, the feature vectors corresponding to high order (10 to 19) as well as low order coefficients (0, 1, 3, 5, 8) consistently give higher AUC values in any of the frame length cases tested. Exceptions occur in each of the 14 th and 18 th coefficient at 20 ms and one in the 19 th at 28 ms but can be considered as outliers because of their very small value. The feature vector representing the 2 nd MFCC, which was found to be the most important of the Φ M set, clearly outperforms the 2 nd LFCC with a difference that reaches −0.57 in the AUC result (absolute median value is 0.5). The same holds for the 4 th and 9 th CC but with a smaller absolute median value of 0.18 and 0.15 respectively. For the trajectory coefficients however, there is no apparent advantage of one feature set over the other that persists with changing frame length, except in very few cases (at the 2 nd and 7 th delta coefficients). , by comparing the AUC scores of the static coefficient and its delta and delta-delta. Each box in the figure consists of the values obtained, for every frame length, by subtracting the AUC of the coefficient indicated by the index on the x-axis from the top performing one found as before. The MFCC static coefficients 0, 2 to 7 and 9 to 12 exhibit better classification performance in over 65% of the frame lengths with the 0 th , 2 nd and 5 th to 7 th scoring consistently higher than their trajectories in any frame length. This is different for the 1 st , 8 th and 13 th to 18 th static CCs where they score lower than either the delta or delta-delta coefficients in over 50% of the frame lengths. The feature vector representing the 19 th MFCC scores similarly with the vectors corresponding to the coefficient's derivatives but its performance degrades for l > 84 ms. For the Φ L feature set, the overwhelming majority of the static CCs score consistently higher than their delta and delta-delta. More specifically, the feature vectors f L i for i = 4, . . . , 20, i = 5, 10, 11 score higher for any l.
Overall, the results show that the higher order LFCC outperform the corresponding MFCC (Fig. 3) but Fig. 2 suggests that, from a classification point of view, it is better to use mel-scaled instead of linear-scaled filters to extract the CC. However, to justify this, further tests are needed in order to investigate the effect on the performance, the number of filters in the filter-bank has. This is explored in Section IV-F.
2) Comparison of STFT and compressed STFT feature sets: All three feature sets, Φ F , Φ E and Φ D , achieve comparable maximum performances in AUC, with values 0.846 (l = 24 ms), 0.839 (l = 88 ms) and 0.850 (l = 20 ms) respectively. At l = 24 ms the subset selected from Φ F that gave the highest AUC only consists of features that fall within the range 0.29 to 2 kHz. Comparing the performances of the individual feature vectors (i.e. the 11-dimensional elements of φ E , φ D and φ F ) within this range, shows that the φ F feature vectors perform better than the corresponding from φ E and φ D . This is depicted in Fig. 5 that shows the average AUC scores per feature vector, against frequency. At approximately 3.5 kHz however, the performance starts to drop rapidly, becoming comparable to, or worse than, the corresponding features from the Φ E and Φ D sets. At around 6 kHz the performance starts to improve again, reaching a maximum near 7.2 kHz at a value of 0.75, comparable to that of the lower frequency features (≤ 3.5 kHz). Above 4.5 kHz it is observed that the classification performance actually improves when using the static coefficients or their second derivatives from the Φ E and Φ D sets. The line plots in Fig. 5 follow similar trends.
Adding to the above, the static coefficients from both Φ E and Φ D perform better than their corresponding derivatives for frequencies up to 2 kHz. The delta coefficients peak in the range 2.2 to 2.6 kHz but then begin to drop rapidly whereas for higher frequencies (≥ 5 kHz) the delta-delta coefficients perform better overall. This shows that additional information exists in the dynamics of the spectrum for mid to high frequencies. The above observations can also be derived from the results at other, higher values of l. From the plots of Fig. 6 it becomes apparent that, with increasing frame length, the delta coefficients capture more information than the corresponding static and delta-delta coefficients in the frequency range 1.6 to 3 kHz.
With increasing frame length, the frequency resolution improves and it becomes immediately apparent from Fig. 6 that two regions of high performing features exist, one in the range [220, 420] Hz, denoted as F r1 , and another in the range [1, 3.4] kHz, denoted as F r2 . These frequency regions show a collection of features that individually score higher than the rest in their respective set (e.g. ≥ 0.726 AUC in the Φ F set). It is also worth noticing that for l ≤ 36 ms there exists another frequency band (6.6 to 7.6 kHz), where some feature vectors that fall within the band, from all three sets (Φ F , Φ E , Φ D ), also scored high AUC (≥ 0.75). However, in the three examples of Fig. 6 it is seen that their performance gradually drops with increasing frame length to, or below, AUC = 0.5. This means that the classifier randomly assigns observations to classes (AUC = 0.5), given the input features, or the classifier failed to apply the information at hand correctly (AUC < 0.5).
For the second case one may reverse the classifier's decisions and obtain a ROC curve that would give AUC > 0.5, as long as the classifier consistently produces results falling in the lower right part of the ROC space as described in [48]. Feature vectors from Φ F that fall within F r1 consistently perform better than the corresponding vectors from both Φ D and Φ E in all frame-sizes tested. On the contrary, the performance of the feature vectors from the latter two sets, at frequencies that fall within F r2 , increases with increasing l and becomes similar to or even better than those of the Φ F set. As can be seen in Fig. 6 (c), for frequencies between 1.5 kHz and 3 kHz (within F r2 ), the delta coefficients from both Φ E and Φ D sets outperform the rest, achieving AUC ≥ 0.8. This effect is reflected in the performance curve (Fig. 2) of the Φ E and Φ D feature sets as a slightly upward trend compared to the degrading performance with the Φ F set.
The feature subsets from the Φ F set, selected based on the method described in Section III, contain only features that fall within F r1 or F r2 for l ≥ 24 and only from F r1 for l ≥ 32. For the Φ D and Φ E sets, the majority of the feature vectors included in the selected subsets fall within the two frequency bands for any l tested. Therefore, these observations highlight the importance of the spectrum features that fall within the two identified bands of frequencies and show that these features have a strong impact on the classification performance. In the analysis of our previous work, [35], it was observed that the top performing features obtained from the magnitude spectrum are primarily at the low frequencies. This is further supported here where we have also identified specific bands that carry significantly discriminant information.
From the above results and observations it can be concluded that the information in the time-frequency spectrum that enables an SVM classifier to discriminate between the two classes (normal/healthy vs abnormal/OA) is contained in a range of frequencies but in specific bands that depend on l. It can be deduced that the frame-size, and hence the DFT length, is a classification performance trade-off. Small l is preferred in order to achieve suitable frequency resolution for the knee sounds occurring in the range of 0.7 to 3.5 kHz and at frequencies ≥ 6 kHz but a larger l is preferred to capture the finer details of the spectrum in the 220−420 Hz band. For the sets Φ D and Φ E , the latter two hold true. However, for the range 0.7 to 3.5 kHz the choice of l does not significantly affect the maximum AUC values achieved but it has an effect in the performance of the individual coefficients.
The suitability of the choice N B = 20 is tested in Section IV-F, even though reducing the dimensionality of the STFT spectrum using 20 filters scaled either linearly or non-linearly (mel) in frequency improves the AUC in 20 out of the 21 l values tested. An adaptable filter-bank with narrower filters at the low (< 1 kHz) and highest frequencies (≥ 6 kHz) and broader at the mid to high frequencies could capture, with fewer coefficients than the full resolution spectrum, the information needed to discriminate the normal and abnormal signals with higher than 0.85 in AUC. This is further supported by the results showing that the Φ D low frequency features (up to 500 Hz or 1.6 kHz depending on l) consistently perform better than the corresponding Φ E whereas at high frequencies the difference in AUC is diminished and sometimes reversed. These observations are further investigated and validated in Section IV-F.

D. Local search in the vicinity of the best frame length
The l values that gave the maximum AUC, denoted as l b,x for the frame length b of the x feature set, were found from the experiments in the previous Section in which the time step used was t = 4 ms. The existence of local maxima in the vicinity of these points can be tested by defining a grid of 6 values with a time step of 1 ms centred at l b,x i.e. l b,x − t for t = 3, 2, . . . , −3, t = 0. Steps shorter than 1 ms will not have an impact in the result given that, at F s = 16 kHz, the sample difference would be less than 16 samples.
The same experimental framework is followed and the training and test sets are standardized as before. In this experiment, 3 more classifiers are used to evaluate the subsets created at the feature selection step in order to assess how well the subsets generalise with different classifiers. These are the (a) LDA classifier that finds linear hyperplanes in the feature space which separate the two classes, (b) CART and (c) SVM with a Gaussian kernel (SVM g ) in order to look for more complex and non-linear boundaries in the feature space. Classification results are again evaluated using AUC.
From these experiments it was found that, for the sets Φ M and Φ L , SVM l and LDA generated comparable results at approximately AUC= 0.9 and 0.8 in the respective sets, with SVM l being slightly better by 1% for Φ M and 0.5% to 7% for Φ L . SVM g and CART achieved AUC scores that ranged from 9% to 55% lower than the maximum. For both Φ E and Φ D sets, SVM l outperformed by at least 5.2% and 10.5% respectively the other classifiers which generated comparable scores. For all l b,x − t in Φ F and 4 out of the 6 values of t in Φ L , SVM g performed significantly worse than the other classifiers by at least 70% and 31% respectively, producing results as low as AUC = 0.36. The power of the SVM classifier with a Gaussian kernel is limited in this case by the relatively large number of features compared to the small training set size, which increases the risk of overfitting when the data is transformed to a high dimensional feature space. Therefore, for the binary classification task of this work and based on the amount of available data, the performance advantage of linear classifiers suggests that the two classes can be linearly separated in the feature spaces explored.
As described in the previous paragraph, SVM l achieves higher AUC with all of the feature sets. However, the features were selected based on results obtained from training and testing this specific type of classifier. Therefore, the subsets created at the feature selection step are tailored to work better with this and similar classifiers. Nevertheless, the work in this paper is not concerned with finding the best classifier to use with the available data. It is rather focused on the efficacy of the features in question to separate the two classes and capture information that would eventually lead to finding the specific abnormality (OA) signatures in the signals. The experiments showed that the average AUC values for the frame-sizes l b,x − t, for t = 3, 2, . . . , −3, t = 0, were close to those of l b,x . The results are improved for all feature sets except for Φ L . Table II reports the frame length value that gave the highest AUC per feature set, found using SVM l . A corresponding score, S c , computed as S c = MCC + (1 − E r ) + F 0.5 /3, [35], is also used and consists of the MCC and F 0.5 measures which capture different attributes of the classification result than the AUC and would therefore be useful in the analysis. S c can vary between 0 and 1 (where S c = 1 indicates perfect prediction). It is shown in Table II that S c ranks the top performance of the feature sets in the same way as AUC.

E. Monte Carlo experiments
In the previous experiments, the l values to be tested were defined in a deterministic approach. In this Section, a stochastic approach is followed in which the values are randomly chosen from a range. In this way the optimality of the result in the previous Section is assessed and the effect of increasing or decreasing the frame length even further is examined.
Firstly, plausible limits for l need to be set. Using frames larger than a single stride will cause overlap of the sounds from two strides resulting in poor modelling of those sounds. In addition, applying the DFT operation will become inappropriate because the signal in a single frame will be nonstationary. The maximum recorded speed is 9 km/h giving an average stride duration of 0.7 s. Given that a single frame length value is applied at the feature extraction stage for all s i,j (n), the upper bound is set at 0.7 s. For the lower bound 2 ms is chosen which gives a time resolution that allows fine localisation of the sounds and approximately 470 Hz frequency resolution at F s = 16 kHz for the DFT.
Monte Carlo simulations are performed by randomly assigning 20 values to l within the range 2 to 700 ms, excluding the range 20 to 100 ms, since it was examined in the previous experiments. The same experimental framework is executed as before using only SVM l and the results are shown in Fig. 7. In comparison to Table II using fixed frame segmentation, only the Φ M feature set generates a classification performance that behaves smoothly over the entire range and has a global maximum. The Φ L feature set on the other hand generates a clear local maximum in the range of 20 to 100 ms. However, for l ≥ 205 ms (except at 272 and 400 ms) the AUC scores are higher than the previously best one obtained at l = 20 ms. The maximum is achieved at l = 567 ms giving AUC = 0.875 which is even higher than that achieved with Φ M . The Φ D and Φ E feature sets achieve an AUC of 0.86, both at 132 ms, whereas for l ≥ 174 ms, AUC values less than 0.81 are achieved with both sets. The Φ F set on the other hand exhibits a more variable behaviour, with the AUC ranging between 0.76 and 0.85. However, it shows a clear global maximum at l = 23 ms.

F. Experiments on the number of filters
The effect of the N B parameter on the classification performance of all the feature sets except Φ F is examined in this section. N B is varied from 10 to 75 while l is kept fixed at the values of Table II. The outcome is evaluated using AUC and S c . The experimental outcomes are summarised in Fig. 8 for each set, showing the highest values obtained by any feature subset in each N B case based on the plots' metric. As expected, the specific values of the two performance metrics are classifier depended. However, the general outcomes and the observations derived are similar across the different classifiers. Fig. 8 (a) and (b) show an overall negative trend with a small variance in the final value in both metrics as N B increases. It can be deduced that there is a strong indication that a small number of filters (≤ 20) is more suitable for the extraction of Ψ D . The two metrics are shown to improve slightly for N B > 50 only with CART but later drop and never exceed the highest score obtained with N B = 14. With SVM l and CART, the subsets that generated the highest results for any N B always included f M 3 . In fact, ∀N B = 10, 11, 13, this particular feature vector was the only one selected in the final subset, one of which generated the overall highest result (AUC = 0.921 for N B = 33 and SVM l ). For LDA, the selected feature subsets included f M 3 together with at most 2 more features (both static and delta coefficients) for N B = 11, 12, . . . , 18 whereas for all other cases f M 3 was the only feature vector in the subset.
The maximum classification performance measured in both AUC and S c is achieved when using only f M 3 . However, if an exhaustive search is performed through the entire Φ M feature set and all possible combinations are used with a classifier, it is likely that a subset containing more vectors than only f M 3 would generate better performance. For this experiment however, such a method is computationally very costly as it would generate possible combinations per classifier. Given that there is only a limited amount of data, the classifier results would become meaningless when n becomes larger than a certain value because the feature space will eventually become very sparse. When this happens, the classifier's decision boundaries will be formed due to the sparseness of the feature space and not due to the information captured by the features which can also lead to overfitting, among other problems, [49]. For these reasons, a suboptimal search method, like the one employed in this work, is favoured and was found to generate good performance as shown in Fig. 8 (a) and (b).
Compared to Φ M , the results of the Φ L feature set are more variable as shown by plots (a) to (d) in Fig. 8. High classification scores are obtained for 13 ≤ N B ≤ 50 (depending on the classifier) and the best result for this set is obtained with SVM l and N B = 17 giving an AUC of 0.853. A slightly upward trend for N B ≥ 51 is observed, which is more noticeable with SVM l . From Fig. 8 (e) and (f), it is clear that, when using LDA with the Φ E set, the classification performance in terms of AUC and S c is higher for N B ≤ 20. With CART, an increase in the values of both metrics is observed for N B ≥ 46 which exceeds that achieved with LDA for those values. This jump is attributed to the inclusion of features falling in the frequency band F r1 compared to only using features that fall within F r2 for N B < 32. This effect is less obvious with SVM l where the scores achieved in the range of the N B values tested are comparable. For the final set, Φ D , the plots are similar amongst the classifiers and suggest that, for best results, a small number of filters (N B ≤ 20) is also preferable. It is observed that the selected subsets with which the classifiers scored the highest AUC and S c only contain feature vectors that fall in either F r1 or F r2 or both. This further supports the observations and results discussed in Section IV-C, stressing the importance of these frequency bands, as well as the importance of the information in the dynamics of the spectrum for classifying normal and abnormal knee signals.

V. CONCLUSIONS
This paper investigates the discriminant power of timefrequency features for the task of classifying knee condition, and explores effective parameterisations of the knee sounds collected from healthy and OA knees during walking. The efficacy of knee condition classification was evaluated by qualitative and quantitative analysis using the AUC of the ROC curve and metrics derived from the confusion matrix (E r , F 0.5 , MCC). Additionally, a study of the effect of varying the values of the feature extraction parameters of frame length and number of filters in the filter-bank was presented in order to examine their impact on the classification performance. The results of our work enable the extraction, from AE signals, of spectrum features and CC focused on specific frequency bands that were shown to carry significantly discriminant information. The answers to the four research questions that were defined in the introduction, as the main aim of this work, are summarised in the following two paragraphs.
The results show that reducing the dimensionality of the STFT spectrum using a mel-spaced triangular filter-bank improves the classification performance compared to using the full resolution spectrum. Furthermore, taking the natural logarithm of the STFT spectrum and subsequently computing the DCT can also improve the performance. The analysis signifies that the results are also improved when using melfrequency scaling instead of linear frequency. In fact, using CART, LDA and SVM as the tools for classification, the findings demonstrate that low order coefficients from the Φ M feature set (especially the f M 3 feature vector) can distinguish between healthy and OA knees with the highest AUC amongst the 5 feature sets examined.
The experiments conducted to investigate the effect of the frame length and the number of filters, revealed two frequency regions, namely 220 to 420 Hz and 1 to 3.4 kHz. These regions contain a collection of features (both static and derivative coefficients) that individually score higher classification results than the rest in their respective feature set. The analysis performed highlighted the importance for classification performance of the spectrum features within these two bands. Finally, the very good classification performance in the experiments (Table II, Fig. 5 to Fig. 8) validates the hypothesis outlined at the beginning of Section II which stated that the acoustic artifacts caused by walking on the treadmill are uncorrelated with the analysed features.
Several fruitful directions for future research can be identified from the conclusions described above. For example, the combination of features from different domains was not explored for classification as the research interest was focused at investigating the feature domains independently in order to identify which ones contained more descriptive information of OA. Therefore, feature domain combination can be an interesting topic for future research. In addition, given the results and the insights derived from the experiments it follows that it will be interesting to examine whether an adaptive filterbank can be beneficial for the analysis.
Contrary to other studies that focus on sit-to-stand movements and similar variants (e.g. knee flexion and extension), this study analysed signals obtained from knees performing a dynamic action. The outcomes presented in this paper suggest that the analysis of such signals can lead to non-invasive detection of knee OA with high accuracy and could potentially extend the range of available tools for the assessment of the disease as a more practical and cost effective method without requiring clinical setups.