Introduction

Chronic obstructive pulmonary disease (COPD) is a heterogeneous disorder with stable periods and acute exacerbations. Exacerbations are characterized by rapid worsening of dyspnea, cough, and/or sputum and drive treatment intensity, hospitalization, and follow-up.¹

Although lung sounds constitute a fundamental component of bedside respiratory assessment, conventional auscultation remains largely subjective. The increasing use of digital stethoscopes and computer-assisted analysis techniques has enabled the quantification of respiratory sounds across time–frequency–energy domains, allowing for more reproducible and standardized measurements; guidelines such as those from Computerized Respiratory Sound Analysis have further supported terminology and measurement standardization in this field.^2,3 Systematic reviews have demonstrated that computerized lung sound analysis holds clinically meaningful potential for detecting abnormal respiratory sounds.⁴ Despite rapid advances in artificial intelligence and machine learning, many studies emphasize classification over pathophysiological interpretability.^5-7

In COPD, airflow limitation, energy instability, heterogeneous small-airway narrowing, and nonlinear vibration may be encoded in the spectral and temporal structure of respiratory sounds.⁶ Therefore, rather than relying on a single metric, the combined evaluation of complementary indices may enhance both discriminatory performance and pathophysiological interpretability.

In this study, we proposed four recording-level indices derived from digital lung sound recordings—Median Respiratory Frequency (MRF), Long-Term Energy Index (LTEI), Subharmonic Index (SubH), and Harmonicity Deviation Index (HDI)—to jointly describe frequency distribution, temporal energy instability, nonlinear/subharmonic structure, and spectral organization, respectively.

Materials and Methods

Study Design and Ethical Approval

This single-center, cross-sectional retrospective study was conducted at the Chest Diseases Clinic of Kafkas University Faculty of Medicine Hospital. We reviewed electronic medical records and digitally archived lung sound recordings obtained during routine care between June 1, 2023 and November 30, 2025, and reported the study in accordance with the Strengthening the Reporting of Observational Studies in Epidemiology initiative.

Data were analyzed in anonymized form. The protocol complied with the Declaration of Helsinki and was approved by the Kafkas University Clinical Research Ethics Committee (KAÜ-80576354-050-099/813; December 2, 2025). The Ethics Committee explicitly waived the requirement for written informed consent, determining that the study met the criteria for waiver: the design was retrospective, all personal identifiers were removed prior to analysis, and the recordings had been obtained exclusively as part of routine clinical care with no study-specific procedures performed.

Adults (≥18 years) were eligible if lung sound recordings were obtained with the routine protocol and corresponding clinical and spirometric data were available. We excluded recordings with insufficient duration, extensive artifacts, or missing audio files, concurrent major pleural pathology or acute cardiovascular events at recording, and incomplete clinical or spirometric data.

Participants and Clinical Groups

We included 235 participants and categorized them as Healthy, Stable COPD, or COPD Exacerbation according to clinical status at the time of recording.

Demographics, smoking history (pack-years), comorbidities, clinical findings, spirometry, and computed tomography (CT)-reported emphysema status corresponding to the recording date were extracted from the hospital information system. COPD diagnosis was confirmed using Global initiative for chronic Obstructive Lung Disease (GOLD) 2024 criteria.⁸

Healthy participants were identified from archived recordings of individuals who had undergone routine occupational health surveillance at our institution. As these individuals attended for occupational screening rather than respiratory complaints, they represented a population free of known respiratory disease at the time of recording, providing an appropriate reference group for acoustic comparison. Eligible healthy participants were required to have no chronic respiratory symptoms, no known chronic lung disease, no evidence of acute lower respiratory tract infection or pneumonia at the time of recording, and normal spirometry. Chest auscultation was performed by a pulmonologist; the absence of pathological lung sounds was defined as the absence of wheeze, crackles, rhonchi, stridor, or focal/asymmetric reduction in breath sounds. Normal spirometry was defined^9,10 as postbronchodilator forced expiratory volume in 1 second (FEV₁) and forced vital capacity (FVC) ≥80% predicted and FEV₁/FVC ≥0.70. Smokers without clinical or spirometric COPD were eligible only if all of these criteria were satisfied.

Stable COPD patients had postbronchodilator FEV₁/FVC <0.70 and were clinically stable at recording, defined as no exacerbation, systemic corticosteroid/antibiotic use, or change in maintenance inhaler therapy within ≥4 weeks, with stable vital signs documented in physician notes.

COPD exacerbation was defined per GOLD 2024 and required recordings obtained during the exacerbation period. Exacerbation was an acute worsening of dyspnea and/or cough and sputum compared with stable state, requiring additional treatment and typically lasting <14 days.¹ Patients on invasive mechanical ventilation were excluded; technically adequate recordings under noninvasive ventilation or oxygen therapy were included.

Digitally archived posterior thoracic recordings were screened and 3 datasets were defined (370 recordings/62 individuals; 508/85; 525/88). Each 15-s recording was analyzed as a region-specific observation.

In COPD, emphysema status was determined qualitatively (present/absent) from routine chest CT reports: cases explicitly reporting emphysema were included, whereas reports stating absent emphysema or not mentioning emphysema were excluded. To provide an objective CT-based assessment of emphysema, available chest CT scans were analyzed using automated lung segmentation with 2 deep learning–based tools (TotalSegmentator and LungMask). After segmentation of the lung parenchyma, emphysema burden was quantified as the percentage of low-attenuation area below -950 Hounsfield units (LAA-950). LAA-950 was calculated separately for each segmentation model, and the mean of the 2 measurements was used as the final emphysema value in the analysis. Image analysis was performed blinded to clinical group assignment. This approach ensured that the study population reflected a radiologically confirmed emphysematous COPD subgroup.

To address the possibility of selection bias related to exclusions, the available baseline characteristics of excluded screened cases were also summarized and reviewed descriptively in relation to the final analytic cohort. Because several exclusions were related to incomplete clinical or spirometric data, the excluded cases were not treated as a formal comparison group in the primary analyses.

Lung Sound Recording Protocol

Lung sound recordings were obtained under standardized conditions using a Littmann CORE digital stethoscope (3M; St. Paul, Minnesota); all participants were seated and instructed to breathe deeply and rhythmically through the mouth. Recordings were performed in a quiet outpatient examination room or, for hospitalized patients, in a single-occupancy inpatient room with the door closed to minimize ambient noise; a seated position was maintained in both settings. Patients with COPD exacerbation who were receiving intermittent noninvasive ventilation (NIV) were recorded during a period when NIV was not in use; no recordings were obtained while NIV was actively applied. Auscultation was performed at 6 posterior thoracic sites (bilaterally at the interscapular, midscapular, and infrascapular levels). A 15-second recording was obtained from each site, yielding a total of 90 seconds of lung sound data per participant. All recordings were digitally stored with a sampling frequency of 44.1 kHz and 16-bit resolution.

Segments containing prominent artifacts due to speech, coughing, patient or stethoscope movement, or cable contact were identified through combined visual and auditory inspection of waveforms, spectra, and spectrograms and were excluded. If artifacts affected the majority of a 15-second segment, the entire segment was discarded. Only recordings of sufficient technical quality were included in the statistical analyses. Because each 15-second recording contained multiple respiratory cycles, inspiratory and expiratory phases were not manually separated; all acoustic analyses were conducted at the recording level, and the calculated acoustic indices (MRF, LTEI, SubH, and HDI) were expressed as single representative values for each segment.

Signal Preprocessing and Time-Frequency Analysis

Raw lung sound recordings were first processed using a 100–1600Hz band-pass filter to attenuate low-frequency chest wall noise and high-frequency interference. Subsequently, amplitude normalization was applied to reduce inter-recording amplitude variability, and potential DC offsets were corrected to center the signal around zero.

Each 15-second recording was segmented into overlapping windows of 200ms duration with 50% overlap for short-time analysis. For each window, the short-time Fourier transform (STFT) was computed, and the power spectrum was obtained as:

These power spectra served as the primary data source for the calculation of all acoustic markers.

Calculation of Acoustic Parameters

To quantitatively reflect different aspects of airway physiology, 4 acoustic markers were defined: MRF, LTEI, SubH, and HDI. All markers were calculated as median-based summaries of spectral and energy measures derived from short-time frames.

MRF represents the dominant frequency of lung sounds and reflects turbulence related to airway narrowing. For each frame, the power spectrum within the 100–1600Hz range was examined, and the frequency at which the maximum spectral power occurred was defined as the dominant frequency:

For a given 15-second recording, the median of all frame-level dominant frequencies fpeak(i) was recorded as the MRF.

LTEI is an energy-based index reflecting the long-term stability of the respiratory sound signal. Instantaneous energy was calculated for each 200ms frame. The mean energy μE and standard deviation σE across all frames were then computed, and LTEI was defined as the ratio σE/μE.

Across all frames, the mean energy μE and standard deviation σE were calculated, and LTEI was defined as the ratio σE/μE.

This ratio represents the relative degree of temporal energy fluctuation and yields a single LTEI value for each recording.

SubH quantifies the relative contribution of subharmonic components in lung sounds. For each frame, after identifying the dominant frequency f₀, 3 narrow frequency bands (±20Hz) centered at f₀​​/2, f₀​/3, and f₀​/4 were defined. The sum of spectral power within these bands was calculated as subharmonic energy. The total spectral power within the 100–1600Hz band was defined as EtotalE_{\text{total}}Etotal. The frame-level subharmonic ratio was computed as the ratio of subharmonic energy to total energy. For each 15-second recording, the median of all frame-level SubH values was used as the recording-level SubH.

Within the same frame, the total power of all spectral components in the 100–1600Hz band was defined as E_total​:

The frame-level subharmonic ratio was calculated as:

For each 15-second recording, the median of all frame-level SubHi_ii values was used as the SubH value of the corresponding recording:

HDI was defined to quantify spectral irregularity in the harmonic structure. For each frame, the power spectrum Pi(k) was evaluated based on the differences between adjacent frequency bins:

where M denotes the number of frequency bins. In signals with a regular harmonic structure, HDI is expected to be lower, whereas higher values are anticipated in patterns with increased noise and structural distortion. For each recording, the median of all frame-level values HDIi was defined as the HDI value of that recording:

Thus, MRF, LTEI, SubH, and HDI were used as 4 complementary acoustic markers summarizing each 15-second lung sound recording in the time–frequency domain.

Statistical Analysis

Continuous variables are presented as mean ± standard deviation or as median (interquartile range [IQR]), as appropriate. Categorical variables are presented as counts and percentages. Between-group comparisons (Healthy, Stable COPD, COPD Exacerbation) used the chi-square test for categorical variables. For continuous variables, one-way analysis of variance was used when applicable; otherwise, the Kruskal–Wallis test was used. Acoustic analyses were performed on 15-second recording segments. The primary acoustic indices were MRF, LTEI, SubH, and HDI. Overall group differences for these indices were tested using the Kruskal–Wallis test. Pairwise posthoc comparisons were performed using Dunn’s test with Holm adjustment. For multivariate analyses, the 4 indices were z-standardized and projected using principal component analysis (PCA). Group separation in the standardized 4-index space was tested using permutational multivariate analysis of variance (PERMANOVA) with Euclidean distance and permutation-based P-value values. All tests were 2-sided, and P-value<0.05 was considered statistically significant.

Results

A total of 322 archived cases were screened, of which 87 were excluded according to the predefined criteria, leaving 235 participants in the final analytic cohort: 62 healthy individuals, 85 patients with stable COPD, and 88 patients with COPD exacerbation (Supplementary Figure S1 in the online supplement). After recording-level quality control, the acoustic dataset comprised 1403 posterior thoracic recordings of 15 seconds’ duration, including 370 recordings from healthy individuals, 508 from patients with stable COPD, and 525 from patients with COPD exacerbation. Baseline characteristics of the analytic cohort are presented in Table 1.

Quantitative emphysema burden (LAA–950) was significantly higher in the COPD exacerbation group than in the stable COPD group (median 27.61% [IQR 14.18–35.39] versus 15.27% [IQR 4.81–32.99], P-value= 0.002). Detailed quantitative emphysema data are presented in Supplementary Table 1 in the online supplement.

Between-group differences for the 4 prespecified indices (MRF, LTEI, SubH, and HDI) were assessed using the Kruskal–Wallis test, followed by Dunn’s posthoc tests with Holm adjustment for multiple comparisons (Table 2 andf Table 3).

These median patterns are visualized in Figure 1.

MRF

Overall group differences were significant (Kruskal–Wallis p=1.36e-05). MRF was lower in Stable COPD than in Healthy (P-value=3.29e-05) and COPD exacerbation (P-value=0.0007), while Healthy and COPD exacerbation did not differ (P-value=0.2561).

SubH

A pronounced overall difference was observed (Kruskal–Wallis p=3.29e-06). SubH was higher in Stable COPD than in Healthy (P-value=0.0001) and COPD exacerbation (P-value=1.74e-05), whereas Healthy and COPD exacerbation were similar (P-value=0.9121).

HDI

HDI showed the strongest overall separation (Kruskal–Wallis p=4.74e-08). HDI was higher in COPD exacerbation than in Healthy (P-value=1.15e-06) and Stable COPD (P-value=3.09e-06), while Healthy and Stable COPD did not differ (P-value=0.5048).

LTEI

LTEI did not show a statistically significant overall difference across groups (Kruskal–Wallis p=0.0858), and no pairwise comparison remained significant after Holm correction (Table 3).

Multidimensional Separation Using the 4-Index Acoustic Map

To test whether the 4 indices jointly capture a multivariate shift across clinical states, we applied a permutation-based multivariate analysis of variance (PERMANOVA; Euclidean distance on z-scored indices). Overall group separation was significant (pseudo-F=5.19, R²=0.0074, permutation p=0.0001; 9,999 permutations).

A 2-dimensional PCA projection of the z-scored 4-index space showed partially overlapping clusters with modest centroid shifts across clinical groups (Figure S1 in the online supplement), consistent with the small but significant multivariate separation observed by PERMANOVA (Table 4).

A schematic summary of the multidimensional acoustic physiology map is shown in Figure 2.

Discussion

Although lung sounds remain one of the fundamental tools in the clinical assessment of patients with COPD, quantitative analysis of respiratory sounds is a relatively underexplored field, and most existing studies have focused on isolated markers such as the presence of wheezes or crackles, dominant frequency, or a limited number of spectral features.^4,8 In contrast, the present study aimed to evaluate airway physiology within a multidimensional acoustic framework by using 4 composite acoustic indices (MRF, LTEI, SubH, and HDI) derived from digital lung sound recordings obtained from healthy individuals, patients with stable COPD, and patients during COPD exacerbation.

It is well established that respiratory sounds originate from turbulent airflow within the airways and that their frequency content is determined by airway caliber, airflow velocity, and the transmission properties of the lung parenchyma and chest wall.^3,11 Normal lung tissue behaves as a low-pass filter, attenuating higher-frequency components while transmitting lower frequencies more efficiently. Alterations in parenchymal structure, fluid content, or consolidation may modify this filtering effect.¹² Flow-standardized studies have shown that median frequency tends to increase in asthma compared with healthy controls, whereas in COPD it may be similar to or lower than in controls, potentially reflecting changes in parenchymal structure and sound transmission.¹³

In our results, MRF differed among healthy individuals, patients with stable COPD, and those with COPD exacerbation, with the lowest median values observed in the stable COPD group, while the healthy and COPD exacerbation groups exhibited relatively similar median levels (Table 2 and Table 3). This pattern is consistent with a pathophysiological framework in which chronic airway narrowing accompanied by hyperinflation and parenchymal destruction attenuates higher-frequency components generated in central airways, thereby shifting the chest wall-recorded spectrum toward lower frequencies.^3,12,13 During exacerbations, increased turbulence and wheeze episodes in narrowed segments may tend to elevate dominant frequencies, whereas coarse rhonchi, prolonged expiration, secretion retention, and reflex limitation of inspiratory flow may counterbalance this effect. Moreover, the use of median-based recording-level MRF likely limits the influence of brief high-frequency events.^3,4 Together, these opposing forces may account for the absence of a significant difference between Healthy and COPD exacerbation groups: rather than reflecting physiological equivalence, their convergent MRF values likely arise from distinct mechanisms—chronic low-frequency shift in one direction and acute turbulence-driven elevation in the other—that net to a similar recording-level median.

LTEI reflects the temporal stability of respiratory sound energy—a property mechanistically linked to the near-quadratic relationship between instantaneous airflow and recorded sound power,^14,15 such that LTEI can be interpreted as a surrogate for the temporal variability of the airflow–acoustic power transfer function.

In our dataset, LTEI values were numerically higher in the healthy group than in both COPD groups, while stable COPD and COPD exacerbation groups exhibited similar and relatively lower LTEI values (Table 2 and Table 3). However, the overall group effect did not reach statistical significance (P-value= 0.086), and no pairwise comparison remained significant after multiplicity correction; therefore, this pattern should be interpreted as exploratory. This directionality is compatible with the well-established presence of tidal expiratory flow limitation, chronic hyperinflation, and dynamic airway compression in COPD, which may constrain airflow around a plateau despite increased respiratory effort.^8,16,17 Furthermore, physiological studies of lung sound generation and transmission have shown that lung volume, regional ventilation distribution, and parenchymal properties strongly modulate the temporal and spatial characteristics of chest wall-recorded sounds.^18,19 In heterogeneous small airway disease and emphysema, respiratory sound generation may arise from multiple asynchronous microturbulence sources rather than a few dominant bronchial jets; the statistical averaging of these sources, together with the viscoelastic damping properties of emphysematous parenchyma, could result in a flatter, less fluctuating temporal energy profile. Consequently, the observation of lower and similar LTEI values in both stable COPD and COPD exacerbation groups may reflect a “locking” of acoustic output within a narrowed energy band due to chronic airway and parenchymal damage, further reinforced by effort limitation. As this constraining mechanism operates similarly across both COPD states, the resulting narrowing of the acoustic dynamic range may reduce rather than amplify between-group differences, which may account for the absence of statistical significance despite the numerically observed trend.

Unlike binary wheeze detection, SubH quantifies coherent subharmonic energy arising from airway wall vibration and nonlinear dynamics in narrowed small airways. Prior studies using wavelet-based and bicoherence analyses have demonstrated that wheezes represent harmonically and subharmonically rich structures rather than simple single-frequency tones.⁶ More recent deep learning–based lung sound studies have also shown that nonlinear and harmonic features of adventitious sounds, particularly wheezes, contribute meaningfully to classification performance.^20,21

In the present study, median SubH values were highest in the stable COPD group, while healthy individuals and patients during COPD exacerbation exhibited similar and lower values (Table 2 and Table 3). This pattern suggests that SubH does not increase monotonically with clinical severity. Instead, it appears to reflect persistent nonlinear subharmonic patterns associated with chronic small airway remodeling. In stable COPD, heterogeneous luminal narrowing and segmental differences in elasticity may provide favorable conditions for repetitive and relatively stable airway wall flutter, leading to energy accumulation in subharmonic bands. During exacerbations, although the underlying structural substrate persists, edema, increased secretions, and abrupt airway closure–reopening events may disrupt the coherence of these oscillations, rendering subharmonic energy more episodic and temporally dispersed. As a result, the median SubH calculated over 15-second recordings may not increase relative to the stable state. This mechanism may also clarify the apparent convergence between Healthy and COPD Exacerbation groups: while healthy individuals likely exhibit low SubH due to the absence of the structural substrate required for coherent subharmonic generation, exacerbation patients may show similarly low values despite persistent small airway remodeling, as acute disruption of oscillatory coherence may limit the consistent accumulation of subharmonic energy at the recording level. The similar SubH values in these 2 groups may therefore reflect mechanistically distinct states rather than comparable airway physiology.

In normal lungs, particularly during inspiration, respiratory sounds tend to exhibit smoother spectral profiles, reflecting relatively laminar or preturbulent airflow in larger airways and homogeneous parenchymal.^18,22 In contrast, in COPD, small airway narrowing, luminal irregularity, wall thickening, secretion accumulation, and heterogeneous time constants generate multiscale turbulence and flow separation throughout the bronchial tree, resulting in fragmented, “sawtooth-like” spectral patterns.^17,18 HDI thus represents a distinct acoustic dimension capturing the degree of spectral organization versus chaotic fragmentation.²³

Consistent with this framework, median HDI values were lowest in the healthy group and modestly higher in both stable COPD and COPD exacerbation groups (Table 2 and Table 3). In healthy individuals, smoother airway geometry and homogeneous parenchyma preserve spectral continuity at the chest wall, yielding lower HDI values.^18,22 In COPD, patchy airway narrowing, mucosal irregularity, mucus plugs, and emphysema-related elastic heterogeneity generate numerous small turbulence foci that disrupt spectral continuity and increase HDI.^24,25 The lack of a marked difference in HDI between stable COPD and exacerbation groups suggests that HDI is primarily sensitive to chronic structural remodeling and persistent airway geometric irregularity, with acute exacerbation-related processes exerting a comparatively limited additional effect.

The observed separation in MRF, SubH, and HDI alongside the nonsignificant overall behavior of LTEI (Tables 2-4) supports the view that COPD-related acoustic change is distributed across distinct, complementary domains of sound generation and transmission. Together, these 4 domains—frequency shift, energy stability, subharmonic activity, and spectral disorganization—constitute complementary acoustic dimensions of COPD-related airway physiology. The triangular heatmaps (Figure 1) provide a compact, region-resolved view of how these domains vary across clinical states, whereas the 3-axis acoustic physiology model (Figure 2) integrates the indices into an interpretable structure that links each metric to its dominant mechanistic axis. Together, these elements define an acoustic physiology map that organizes complementary information across mechanisms, rather than compressing COPD lung sounds into a single unidimensional summary.

Beyond univariate comparisons, PERMANOVA in the z-scored 4-index space indicated a statistically significant but small multivariate separation among clinical states (Table 4; Figure S1 in the online supplement). Pairwise results suggested that Stable COPD was the most consistently separated state, whereas separation between Healthy and COPD Exacerbation was limited. The low R² values are consistent with the primary aim of this work: rather than classifying clinical states, the 4 indices were designed to characterize the acoustic-physiological mechanisms that underlie each clinical state—a dimension not captured by standard assessment. In this sense, the proposed map complements rather than replaces existing diagnostic tools, with potential relevance for disease monitoring, phenotyping, and the development of mechanistically interpretable digital biomarkers.

Limitations of the Study

Lung sound recordings were obtained without simultaneous airflow or volume measurements (e.g., pneumotachography); consequently, interindividual differences in respiratory rate, tidal volume, and inspiratory flow—all known to influence lung sound frequency content and amplitude—could not be controlled for. This constraint is particularly relevant for MRF, whose dominant frequency is sensitive to airflow velocity, and for LTEI, which captures energy fluctuations in respiratory sounds that are inherently linked to instantaneous airflow variability. Some exacerbation recordings were acquired under oxygen therapy or noninvasive ventilation, which may also have affected acoustic characteristics. Because only COPD cases with radiologically reported emphysema were included, findings may not generalize to nonemphysematous phenotypes. Analyses were conducted primarily at the recording level; multiple recordings per participant may have affected independence assumptions, and future studies should incorporate patient-level validation or models accounting for repeated measures and potential confounders (e.g., age and body mass index). As a single-center study, external validity remains to be established; multicenter validation across diverse COPD phenotypes and clinical settings is warranted. In addition, although a substantial number of screened cases were excluded, the reasons for exclusion were predefined and are now reported transparently, which reduces but does not eliminate the possibility of selection bias.

Acknowledgements

Author contributions: IT and CA were responsible for the study design. IT, CA, and AB were in charge of data collection and data cleaning. GO and MOG were responsible for the clinical review and manuscript review. YB and AMY were in charge of signal processing and software. DSU, SK, and FH were responsible for analysis support. All authors contributed to the writing. All authors approved the final manuscript.

Data availability: The data supporting the findings are available within the article; de-identified data may be provided by the corresponding author upon reasonable request.

Declaration of Interest

The authors declare no competing interests.