Hearbeat Acoustic Analysis

This analysis is performed over heart beat sound wave files, available for public in Kaggle: Source. The dataset consists of two sources: (A) from general public via iStethoscope Pro iPhone and (B) from clinic trial in hospitals, using digital stethoscope DigiScope.

This problem consists of a classification problem, over time-series data, for detection of healthy hearts and abnormal ones.
The study case presented here works on feature engineering over time-series data, as well as spectral engineering and ROC curve.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from glob import glob
from statistics import median
import re

import librosa as lr
from librosa.core import stft, amplitude_to_db
from librosa.display import specshow

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_curve, roc_auc_score, accuracy_score, recall_score
from sklearn.calibration import CalibratedClassifierCV

a_df = pd.read_csv('archive/set_a.csv')
a_df = a_df[['fname', 'label']]
a_df['fname'] = 'archive/'+a_df['fname']
a_df

	fname	label
0	archive/set_a/artifact__201012172012.wav	artifact
1	archive/set_a/artifact__201105040918.wav	artifact
2	archive/set_a/artifact__201105041959.wav	artifact
3	archive/set_a/artifact__201105051017.wav	artifact
4	archive/set_a/artifact__201105060108.wav	artifact
...	...	...
171	archive/set_a/__201108222241.wav	NaN
172	archive/set_a/__201108222244.wav	NaN
173	archive/set_a/__201108222247.wav	NaN
174	archive/set_a/__201108222254.wav	NaN
175	archive/set_a/__201108222257.wav	NaN

176 rows × 2 columns

a_df['label'].value_counts(dropna=False)

NaN         52
artifact    40
murmur      34
normal      31
extrahls    19
Name: label, dtype: int64

b_df = pd.read_csv('archive/set_b.csv')
b_df

	dataset	fname	label	sublabel
0	b	set_b/Btraining_extrastole_127_1306764300147_C...	extrastole	NaN
1	b	set_b/Btraining_extrastole_128_1306344005749_A...	extrastole	NaN
2	b	set_b/Btraining_extrastole_130_1306347376079_D...	extrastole	NaN
3	b	set_b/Btraining_extrastole_134_1306428161797_C...	extrastole	NaN
4	b	set_b/Btraining_extrastole_138_1306762146980_B...	extrastole	NaN
...	...	...	...	...
651	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
652	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
653	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
654	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
655	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal

656 rows × 4 columns

b_df['label'].value_counts(dropna=False)

normal        320
NaN           195
murmur         95
extrastole     46
Name: label, dtype: int64

Database (A) contains 176 audio files and database (B) contains 656 audio files. These records are labeled in four classes:

• Normal: indicates a healthy heart beat
• Artifact (only A): indicates unregular, too noisy, or not heartbeat waveforms
• Murmur: one of abnormalities in heart beat
• Extrahls (A) or Extrastole (B): one of abnormalities in heart beat

Some of the data does not have a label, indicating NaN.

Individual analysis

We will start with individual analysis, over one waveform from each class, from database A, to illustrate the feature engineering process, before moving on to the classification model.

files = ['archive/set_a/normal__201103170121.wav',
         'archive/set_a/artifact__201106050353.wav',
         'archive/set_a/extrahls__201101241423.wav',
         'archive/set_a/murmur__201101180902.wav'
         ]

audios = []
times = []

# Reading auditory data
for file in files:
    audio, sfreq = lr.load(file)
    audios.append(audio)
    
    # Creating time array
    indices = np.arange(0, len(audio))
    times.append(indices/sfreq)

sfreq # sampling frequency (same for every file, prior checked)

22050

df = pd.DataFrame({
    'normal': audios[0][:176400],
    'artifact': audios[1][:176400],
    'extrahls': audios[2][:176400],
    'murmur': audios[3][:176400]
}, index = times[0][:176400])
df

	normal	artifact	extrahls	murmur
0.000000	-0.002086	0.001620	-0.003448	-0.043658
0.000045	-0.002933	0.002505	-0.004786	-0.063429
0.000091	-0.002454	0.002712	-0.003911	-0.055169
0.000136	-0.003028	0.003359	-0.003882	-0.059671
0.000181	-0.002875	0.003370	-0.003303	-0.055605
...	...	...	...	...
7.999773	0.000226	-0.005503	0.000161	-0.040152
7.999819	0.000171	-0.004669	-0.000015	-0.039910
7.999864	0.000138	-0.003795	0.000310	-0.039516
7.999909	0.000653	-0.003394	0.000330	-0.039116
7.999955	0.001223	-0.002117	0.000449	-0.038825

176400 rows × 4 columns

# Ploting waveform
fig, ax = plt.subplots(2, 2, figsize=(16,8))
ax[0,0].plot(df.index, df['normal'])
ax[0,0].set_ylabel('Signal')
ax[0,0].title.set_text('Normal')
ax[0,1].plot(df.index, df['artifact'])
ax[0,1].set_ylabel('Signal')
ax[0,1].title.set_text('Artifact')
ax[1,0].plot(df.index, df['extrahls'])
ax[1,0].set_xlabel('Time')
ax[1,0].set_ylabel('Signal')
ax[1,0].title.set_text('Extrahls')
ax[1,1].plot(df.index, df['murmur'])
ax[1,1].set_xlabel('Time')
ax[1,1].set_ylabel('Signal')
ax[1,1].title.set_text('Murmur')
plt.show()

As noted before, the Extrahls and Murmur indicate two types of abnormalities in heart beats. Artifact files have too much background noise or inaudible heart beats.

Envelope features

Since we are dealing with sound waves, which is a particular kind of time series data, we may work on envelope. It consists of smoothing and averaging over rolling window.

# Feature engineering of the data: smoothing and envelope
window_size = 500
audio_rectified = df.apply(np.abs)
audio_envelope = audio_rectified.rolling(window_size).mean()

# Ploting waveform
fig, ax = plt.subplots(2, 2, figsize=(16,8))
ax[0,0].plot(audio_envelope.index, audio_envelope['normal'])
ax[0,0].set_ylabel('Signal')
ax[0,0].title.set_text('Normal')
ax[0,1].plot(audio_envelope.index, audio_envelope['artifact'])
ax[0,1].set_ylabel('Signal')
ax[0,1].title.set_text('Artifact')
ax[1,0].plot(audio_envelope.index, audio_envelope['extrahls'])
ax[1,0].set_ylabel('Signal')
ax[1,0].set_xlabel('Time')
ax[1,0].title.set_text('Extrahls')
ax[1,1].plot(audio_envelope.index, audio_envelope['murmur'])
ax[1,1].set_ylabel('Signal')
ax[1,1].set_xlabel('Time')
ax[1,1].title.set_text('Murmur')
plt.show()

# Rolling window
fig, ax = plt.subplots(2, 2, figsize=(16,8))
ax[0,0].plot(audio_envelope.index, audio_envelope['normal'].rolling(int(sfreq/2)).mean())
ax[0,0].set_ylabel('Signal')
ax[0,0].title.set_text('Normal')
ax[0,1].plot(audio_envelope.index, audio_envelope['artifact'].rolling(int(sfreq/2)).mean())
ax[0,1].set_ylabel('Signal')
ax[0,1].title.set_text('Artifact')
ax[1,0].plot(audio_envelope.index, audio_envelope['extrahls'].rolling(int(sfreq/2)).mean())
ax[1,0].set_ylabel('Signal')
ax[1,0].set_xlabel('Time')
ax[1,0].title.set_text('Extrahls')
ax[1,1].plot(audio_envelope.index, audio_envelope['murmur'].rolling(int(sfreq/2)).mean())
ax[1,1].set_ylabel('Signal')
ax[1,1].set_xlabel('Time')
ax[1,1].title.set_text('Murmur')
plt.show()

# Feature engineering on envelope
envelope_mean = np.mean(audio_envelope, axis=0)
envelope_std = np.std(audio_envelope, axis=0)
envelope_max = np.max(audio_envelope, axis=0)

rol_mean = np.mean(audio_envelope.rolling(int(sfreq/2)).mean())
rol_std = np.std(audio_envelope.rolling(int(sfreq/2)).mean())
rol_max = np.max(audio_envelope.rolling(int(sfreq/2)).mean())

display(pd.DataFrame([envelope_mean, envelope_std, envelope_max, rol_mean, rol_std, rol_max],
                    index=['mean', 'std', 'max', 'rol_mean', 'rol_std', 'rol_max']))

	normal	artifact	extrahls	murmur
mean	0.007253	0.004898	0.009546	0.055043
std	0.029361	0.003649	0.019131	0.069378
max	0.249543	0.028115	0.113429	0.519126
rol_mean	0.006898	0.004859	0.009663	0.056239
rol_std	0.004376	0.002886	0.003056	0.019750
rol_max	0.013228	0.010245	0.017068	0.094531

Tempo features

Tempo is another feature we can extract from waveform data and relates to number of beats per minute.

# Calculate the tempo of the sounds
tempos = []
for col, i_audio in df.items():
    tempos.append(lr.beat.tempo(y=i_audio.values, sr=sfreq, hop_length=2**6, aggregate=None))

# Convert the list to an array so we can manipulate it more easily
tempos = np.array(tempos)

# Calculate statistics of each tempo
tempos_mean = tempos.mean(axis=-1)
tempos_std = tempos.std(axis=-1)
tempos_max = tempos.max(axis=-1)

display(pd.DataFrame([tempos_mean, tempos_std, tempos_max],
                    index=['mean', 'std', 'max']))

	0	1	2	3
mean	76.511403	123.070574	112.034087	55.473390
std	1.331489	29.969391	19.965717	0.919956
max	79.507212	156.605114	136.899834	57.262812

Spectral features

Finally the most important one, spectral analysis over waveforms.
By applying STFT (short-time Fourier transform) over the waveform, we can extract the spectogram, with frequency content over the time.

# Calculate our STFT
HOP_LENGTH = 2**9
SIZE_WINDOW = 2**12

specs = []
specs_db = []

for audio in audios:
    audio_spec = stft(audio, hop_length=HOP_LENGTH)
    spec_db = amplitude_to_db(np.abs(audio_spec)) # Convert into decibels for visualization
    specs.append(audio_spec)
    specs_db.append(spec_db)

# Visualize
fig, ax = plt.subplots(2, 2, figsize=(16,8))
specshow(specs_db[0], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[0,0])
ax[0,0].title.set_text('Normal')

specshow(specs_db[1], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[0,1])
ax[0,1].title.set_text('Artifact')

specshow(specs_db[2], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[1,0])
ax[1,0].title.set_text('Extrahls')

specshow(specs_db[3], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[1,1])
ax[1,1].title.set_text('Murmur')
plt.show()

# Calculate the spectral centroid and bandwidth for the spectograms
bandwidths = []
centroids = []
times_spec = []

for i, spec in enumerate(specs):
    bandwidths.append(lr.feature.spectral_bandwidth(S=np.abs(spec))[0])
    centroids.append(lr.feature.spectral_centroid(S=np.abs(spec))[0])
    times_spec.append(np.linspace(0, times[i][-1], len(centroids[i])))
    
# Display these features on top of the spectogram
fig, ax = plt.subplots(2, 2, figsize=(16,8))
specshow(specs_db[0], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[0,0])
ax[0,0].plot(times_spec[0], centroids[0])
ax[0,0].fill_between(times_spec[0], centroids[0]-bandwidths[0]/2, centroids[0]+bandwidths[0]/2, alpha=0.5)
ax[0,0].set(ylim=[0, 11000])
ax[0,0].set(xlim=[0, times[0][-1]])
ax[0,0].title.set_text('Normal')

specshow(specs_db[1], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[0,1])
ax[0,1].plot(times_spec[1], centroids[1])
ax[0,1].fill_between(times_spec[1], centroids[1]-bandwidths[1]/2, centroids[1]+bandwidths[1]/2, alpha=0.5)
ax[0,1].set(ylim=[0, 11000])
ax[0,1].set(xlim=[0, times[1][-1]])
ax[0,1].title.set_text('Artifact')

specshow(specs_db[2], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[1,0])
ax[1,0].plot(times_spec[2], centroids[2])
ax[1,0].fill_between(times_spec[2], centroids[2]-bandwidths[2]/2, centroids[2]+bandwidths[2]/2, alpha=0.5)
ax[1,0].set(ylim=[0, 11000])
ax[1,0].set(xlim=[0, times[2][-1]])
ax[1,0].title.set_text('Extrahls')

specshow(specs_db[3], sr=sfreq, x_axis='time',
         y_axis='hz', hop_length=HOP_LENGTH, ax=ax[1,1])
ax[1,1].plot(times_spec[3], centroids[3])
ax[1,1].fill_between(times_spec[3], centroids[3]-bandwidths[3]/2, centroids[3]+bandwidths[3]/2, alpha=0.5)
ax[1,1].set(ylim=[0, 11000])
ax[1,1].set(xlim=[0, times[3][-1]])
ax[1,1].title.set_text('Murmur')
plt.show()

w = 48

# Rolling window of spectral content
fig, ax = plt.subplots(2, 2, figsize=(16,8))
ax[0,0].plot(times_spec[0], pd.DataFrame(centroids[0]).rolling(w).mean())
ax[0,0].set_ylabel('Hz')
ax[0,0].title.set_text('Normal')
ax[0,1].plot(times_spec[1], pd.DataFrame(centroids[1]).rolling(w).mean())
ax[0,1].set_ylabel('Hz')
ax[0,1].title.set_text('Artifact')
ax[1,0].plot(times_spec[2], pd.DataFrame(centroids[2]).rolling(w).mean())
ax[1,0].set_ylabel('Hz')
ax[1,0].set_xlabel('Time')
ax[1,0].title.set_text('Extrahls')
ax[1,1].plot(times_spec[3], pd.DataFrame(centroids[3]).rolling(w).mean())
ax[1,1].set_ylabel('Hz')
ax[1,1].set_xlabel('Time')
ax[1,1].title.set_text('Murmur')
plt.show()

Classification model

Feature extraction of training set

For our classification model, we will develop a simple binary classifier. We will use files of normal and murmur labels only. We will apply database (B) as training set and database (A) as test set, to evaluate our classifier. As classification features, we will apply all those presented previously, over the envelope, tempo and frequency domains.

files = glob('archive/set_b/*.wav')

df = b_df[b_df['label'].isin(['normal','murmur'])]
df

	dataset	fname	label	sublabel
46	b	set_b/Btraining_murmur_112_1306243000964_A.wav	murmur	NaN
47	b	set_b/Btraining_murmur_112_1306243000964_B.wav	murmur	NaN
48	b	set_b/Btraining_murmur_112_1306243000964_D.wav	murmur	NaN
49	b	set_b/Btraining_murmur_116_1306258689913_A.wav	murmur	NaN
50	b	set_b/Btraining_murmur_116_1306258689913_C.wav	murmur	NaN
...	...	...	...	...
651	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
652	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
653	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
654	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
655	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal

415 rows × 4 columns

for file in files:
    try:
        file_number = re.search('([0-9]+.*)', file)[0]
        ind = df[df['fname'].str.contains(file_number)].index[0]
        df.loc[ind, 'file'] = file
    except:
        pass
    
df = df[~df['file'].isna()]
df.set_index('file', drop=True, inplace=True)

/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1599: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.obj[key] = infer_fill_value(value)
/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1720: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._setitem_single_column(loc, value, pi)
/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1637: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._setitem_single_block(indexer, value, name)
/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:692: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  iloc._setitem_with_indexer(indexer, value, self.name)

df

	dataset	fname	label	sublabel
file
archive/set_b/murmur__112_1306243000964_A.wav	b	set_b/Btraining_murmur_112_1306243000964_A.wav	murmur	NaN
archive/set_b/murmur__112_1306243000964_B.wav	b	set_b/Btraining_murmur_112_1306243000964_B.wav	murmur	NaN
archive/set_b/murmur__112_1306243000964_D.wav	b	set_b/Btraining_murmur_112_1306243000964_D.wav	murmur	NaN
archive/set_b/murmur__116_1306258689913_A.wav	b	set_b/Btraining_murmur_116_1306258689913_A.wav	murmur	NaN
archive/set_b/murmur__116_1306258689913_C.wav	b	set_b/Btraining_murmur_116_1306258689913_C.wav	murmur	NaN
...	...	...	...	...
archive/set_b/normal_noisynormal_284_1311168471850_A.wav	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
archive/set_b/normal_noisynormal_284_1311168471850_B.wav	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
archive/set_b/normal_noisynormal_285_1311169246969_C.wav	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
archive/set_b/normal_noisynormal_296_1311682952647_C.wav	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal
archive/set_b/normal_noisynormal_296_1311682952647_D.wav	b	set_b/Btraining_normal_Btraining_noisynormal_2...	normal	noisynormal

412 rows × 4 columns

window_size = 500
HOP_LENGTH = 2**9
SIZE_WINDOW = 2**12
i = 0

features = pd.DataFrame(columns=['envelope_mean','envelope_std','envelope_max',
                                 'rol_env_mean','rol_env_std','rol_env_max',
                                 'tempo_mean','tempo_std','tempo_max',
                                 'band_mean','cent_mean',
                                 'rol_spec_mean','rol_spec_std','rol_spec_max'])

for file, row in df.iterrows():
    i += 1
    if i%50 == 0:
        print(f'{file} : {i}/412')
    
    # Reading audio file
    audio, sfreq = lr.load(file)
        
    # Creating time array
    indices = np.arange(0, len(audio))
    time = indices/sfreq
    
    # Envelope features
    df_audio = pd.DataFrame(audio)
    audio_envelope = df_audio.apply(np.abs).rolling(window_size).mean()
    envelope_mean = np.mean(audio_envelope, axis=0)[0]
    envelope_std = np.std(audio_envelope, axis=0)[0]
    envelope_max = np.max(audio_envelope, axis=0)[0]
    features.loc[file, 'envelope_mean'] = envelope_mean
    features.loc[file, 'envelope_std'] = envelope_std
    features.loc[file, 'envelope_max'] = envelope_max
    
    # Rolling envelope features
    rol_env_mean = np.mean(audio_envelope.rolling(int(sfreq/2)).mean())[0]
    rol_env_std = np.std(audio_envelope.rolling(int(sfreq/2)).mean())[0]
    rol_env_max = np.max(audio_envelope.rolling(int(sfreq/2)).mean())[0]
    features.loc[file, 'rol_env_mean'] = rol_env_mean
    features.loc[file, 'rol_env_std'] = rol_env_std
    features.loc[file, 'rol_env_max'] = rol_env_max
    
    # Tempo features
    tempo = lr.beat.tempo(y=audio, sr=sfreq, hop_length=2**6, aggregate=None)
    tempo_mean = tempo.mean(axis=-1)
    tempo_std = tempo.std(axis=-1)
    tempo_max = tempo.max(axis=-1)
    features.loc[file, 'tempo_mean'] = tempo_mean
    features.loc[file, 'tempo_std'] = tempo_std
    features.loc[file, 'tempo_max'] = tempo_max
    
    # Spectral features
    spec = stft(audio, hop_length=HOP_LENGTH)
    centroid = lr.feature.spectral_centroid(S=np.abs(spec))
    band_mean = np.mean(lr.feature.spectral_bandwidth(S=np.abs(spec)))
    cent_mean = np.mean(centroid)
    features.loc[file, 'band_mean'] = band_mean
    features.loc[file, 'cent_mean'] = cent_mean
    
    # Rolling spectral features
    rol_spec_mean = np.mean(pd.DataFrame(centroid[0]).rolling(48).mean())[0]
    rol_spec_std = np.std(pd.DataFrame(centroid[0]).rolling(48).mean())[0]
    rol_spec_max = np.max(pd.DataFrame(centroid[0]).rolling(48).mean())[0]
    features.loc[file, 'rol_spec_mean'] = rol_spec_mean
    features.loc[file, 'rol_spec_std'] = rol_spec_std
    features.loc[file, 'rol_spec_max'] = rol_spec_max
    
print(f'{file} : {i}/412')

archive/set_b/murmur__244_1309198148498_B.wav : 50/412
archive/set_b/normal__137_1306764999211_C.wav : 100/412
archive/set_b/normal__170_1307970562729_A.wav : 150/412
archive/set_b/normal__209_1308162216750_A.wav : 200/412
archive/set_b/normal__282_1311166081161_C.wav : 250/412
archive/set_b/normal_noisynormal_107_1305654946865_A.wav : 300/412
archive/set_b/normal_noisynormal_137_1306764999211_D1.wav : 350/412
archive/set_b/normal_noisynormal_264_1309356143724_C.wav : 400/412
archive/set_b/normal_noisynormal_296_1311682952647_D.wav : 412/412

features

	envelope_mean	envelope_std	envelope_max	rol_env_mean	rol_env_std	rol_env_max	tempo_mean	tempo_std	tempo_max	band_mean	cent_mean	rol_spec_mean	rol_spec_std	rol_spec_max
Unnamed: 0
archive/set_b/murmur__112_1306243000964_A.wav	0.029616	0.039656	0.256565	0.029443	0.009835	0.052083	122.563272	3.171611	142.564655	388.663388	307.757889	305.204159	18.299142	340.806170
archive/set_b/murmur__112_1306243000964_B.wav	0.080671	0.117259	0.751429	0.081459	0.023949	0.148503	142.623971	38.528759	196.875000	307.096908	246.360413	246.321898	39.731223	334.974950
archive/set_b/murmur__112_1306243000964_D.wav	0.039154	0.065133	0.495011	0.039050	0.014602	0.073650	127.124808	26.138088	186.233108	351.627674	303.114160	303.192225	36.031353	385.726461
archive/set_b/murmur__116_1306258689913_A.wav	0.040288	0.035095	0.269973	0.040634	0.008101	0.058151	112.359995	28.192894	150.889599	352.644700	283.358025	277.767161	12.194930	307.199390
archive/set_b/murmur__116_1306258689913_C.wav	0.016960	0.014725	0.125012	0.016975	0.004363	0.026639	129.135469	14.767350	154.267724	427.153801	377.463748	377.961077	29.499340	456.322541
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
archive/set_b/normal_noisynormal_284_1311168471850_A.wav	0.036412	0.042268	0.228297	0.036477	0.006010	0.051108	121.177779	23.834819	165.375000	354.665334	322.932842	322.194223	36.420620	437.845739
archive/set_b/normal_noisynormal_284_1311168471850_B.wav	0.039781	0.051964	0.364323	0.039622	0.006845	0.054320	99.646657	4.939613	126.821319	384.058825	389.328554	388.813428	22.038180	445.606226
archive/set_b/normal_noisynormal_285_1311169246969_C.wav	0.018664	0.025828	0.146681	0.018568	0.004709	0.030776	115.262789	1.411450	118.803879	491.346528	424.323041	429.251745	67.338075	571.937269
archive/set_b/normal_noisynormal_296_1311682952647_C.wav	0.048895	0.052445	0.436355	0.048987	0.012326	0.096403	127.684934	9.144677	136.899834	395.547843	377.467977	377.268749	27.595868	463.881173
archive/set_b/normal_noisynormal_296_1311682952647_D.wav	0.051436	0.060907	0.288201	0.051921	0.014091	0.113271	132.010060	1.342449	133.366935	355.760134	333.697390	330.987210	20.294923	386.362722

412 rows × 14 columns

files_na = ['archive/set_b/murmur__171_1307971016233_E.wav',
            'archive/set_b/normal__286_1311170606028_A1.wav',
            'archive/set_b/normal__296_1311682952647_A1.wav']

features = features[~features.index.isin(files_na)]
df = df[~df.index.isin(files_na)]

X_train = features
y_train = df['label']
y_train[y_train == 'murmur'] = 1
y_train[y_train == 'normal'] = 0
y_train = y_train.astype('int')

/home/silas/.local/lib/python3.8/site-packages/pandas/core/series.py:992: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._where(~key, value, inplace=True)

Feature extraction of test set

df2 = a_df[a_df['label'].isin(['normal','murmur'])]
df2.set_index('fname', inplace=True)
df2

	label
fname
archive/set_a/murmur__201101051104.wav	murmur
archive/set_a/murmur__201101051108.wav	murmur
archive/set_a/murmur__201101051114.wav	murmur
archive/set_a/murmur__201101180902.wav	murmur
archive/set_a/murmur__201102051443.wav	murmur
...	...
archive/set_a/normal__201106221450.wav	normal
archive/set_a/normal__201108011112.wav	normal
archive/set_a/normal__201108011114.wav	normal
archive/set_a/normal__201108011115.wav	normal
archive/set_a/normal__201108011118.wav	normal

65 rows × 1 columns

window_size = 500
HOP_LENGTH = 2**9
SIZE_WINDOW = 2**12
i = 0

features2 = pd.DataFrame(columns=['envelope_mean','envelope_std','envelope_max',
                                 'rol_env_mean','rol_env_std','rol_env_max',
                                 'tempo_mean','tempo_std','tempo_max',
                                 'band_mean','cent_mean',
                                 'rol_spec_mean','rol_spec_std','rol_spec_max'])

for file, row in df2.iterrows():
    i += 1
    if i%10 == 0:
        print(f'{file} : {i}/65')
    
    # Reading audio file
    audio, sfreq = lr.load(file)
        
    # Creating time array
    indices = np.arange(0, len(audio))
    time = indices/sfreq
    
    # Envelope features
    df_audio = pd.DataFrame(audio)
    audio_envelope = df_audio.apply(np.abs).rolling(window_size).mean()
    envelope_mean = np.mean(audio_envelope, axis=0)[0]
    envelope_std = np.std(audio_envelope, axis=0)[0]
    envelope_max = np.max(audio_envelope, axis=0)[0]
    features2.loc[file, 'envelope_mean'] = envelope_mean
    features2.loc[file, 'envelope_std'] = envelope_std
    features2.loc[file, 'envelope_max'] = envelope_max
    
    # Rolling envelope features
    rol_env_mean = np.mean(audio_envelope.rolling(int(sfreq/2)).mean())[0]
    rol_env_std = np.std(audio_envelope.rolling(int(sfreq/2)).mean())[0]
    rol_env_max = np.max(audio_envelope.rolling(int(sfreq/2)).mean())[0]
    features2.loc[file, 'rol_env_mean'] = rol_env_mean
    features2.loc[file, 'rol_env_std'] = rol_env_std
    features2.loc[file, 'rol_env_max'] = rol_env_max
    
    # Tempo features
    tempo = lr.beat.tempo(y=audio, sr=sfreq, hop_length=2**6, aggregate=None)
    tempo_mean = tempo.mean(axis=-1)
    tempo_std = tempo.std(axis=-1)
    tempo_max = tempo.max(axis=-1)
    features2.loc[file, 'tempo_mean'] = tempo_mean
    features2.loc[file, 'tempo_std'] = tempo_std
    features2.loc[file, 'tempo_max'] = tempo_max
    
    # Spectral features
    spec = stft(audio, hop_length=HOP_LENGTH)
    centroid = lr.feature.spectral_centroid(S=np.abs(spec))
    band_mean = np.mean(lr.feature.spectral_bandwidth(S=np.abs(spec)))
    cent_mean = np.mean(centroid)
    features2.loc[file, 'band_mean'] = band_mean
    features2.loc[file, 'cent_mean'] = cent_mean
    
    # Rolling spectral features
    rol_spec_mean = np.mean(pd.DataFrame(centroid[0]).rolling(48).mean())[0]
    rol_spec_std = np.std(pd.DataFrame(centroid[0]).rolling(48).mean())[0]
    rol_spec_max = np.max(pd.DataFrame(centroid[0]).rolling(48).mean())[0]
    features2.loc[file, 'rol_spec_mean'] = rol_spec_mean
    features2.loc[file, 'rol_spec_std'] = rol_spec_std
    features2.loc[file, 'rol_spec_max'] = rol_spec_max
    
print(f'{file} : {i}/65')

archive/set_a/murmur__201104291843.wav : 10/65
archive/set_a/murmur__201108222235.wav : 20/65
archive/set_a/murmur__201108222252.wav : 30/65
archive/set_a/normal__201102260502.wav : 40/65
archive/set_a/normal__201104122156.wav : 50/65
archive/set_a/normal__201106221418.wav : 60/65
archive/set_a/normal__201108011118.wav : 65/65

features2

	envelope_mean	envelope_std	envelope_max	rol_env_mean	rol_env_std	rol_env_max	tempo_mean	tempo_std	tempo_max	band_mean	cent_mean	rol_spec_mean	rol_spec_std	rol_spec_max
Unnamed: 0
archive/set_a/murmur__201101051104.wav	0.006933	0.010810	0.058625	0.006836	0.002482	0.011689	118.779842	33.241787	161.499023	2073.421953	1299.849966	1299.922538	116.808762	1608.506848
archive/set_a/murmur__201101051108.wav	0.004211	0.010139	0.053660	0.004147	0.002112	0.008425	130.200238	0.659588	131.667994	2339.223848	2601.622592	2613.825375	196.012075	3001.610146
archive/set_a/murmur__201101051114.wav	0.035074	0.037657	0.312929	0.034892	0.010351	0.057709	121.143638	11.161158	147.656250	1925.549710	966.560701	956.785975	43.127071	1045.849637
archive/set_a/murmur__201101180902.wav	0.054968	0.069239	0.519126	0.056127	0.019774	0.094531	55.481259	0.925666	57.262812	997.225121	242.264160	235.430952	28.805209	306.117101
archive/set_a/murmur__201102051443.wav	0.017667	0.019476	0.163166	0.018079	0.003972	0.028357	120.750035	10.057531	145.576585	1901.542487	1102.638283	1098.912886	73.235945	1246.264629
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
archive/set_a/normal__201106221450.wav	0.016585	0.030484	0.195280	0.016409	0.004014	0.026297	115.956375	13.170953	179.755435	1469.167352	668.009789	667.748721	50.712077	812.783618
archive/set_a/normal__201108011112.wav	0.028129	0.055734	0.328286	0.028913	0.010081	0.044716	123.866022	0.233707	124.529367	556.484492	164.342986	157.089317	9.853297	192.038193
archive/set_a/normal__201108011114.wav	0.020082	0.039920	0.261232	0.020800	0.010132	0.035562	128.255823	3.934739	148.718525	602.256006	179.089155	171.955421	7.062927	199.609915
archive/set_a/normal__201108011115.wav	0.025655	0.046168	0.264670	0.025057	0.008818	0.037298	121.357722	0.336778	121.599265	553.202637	156.667061	158.183435	8.488763	179.160240
archive/set_a/normal__201108011118.wav	0.038313	0.059735	0.270689	0.038239	0.008150	0.049307	90.329693	0.141479	90.666118	488.965049	137.959979	133.220252	5.805155	152.869878

65 rows × 14 columns

file_na = ['archive/set_a/murmur__201104021355.wav']

features2 = features2[~features2.index.isin(file_na)]
df2 = df2[~df2.index.isin(file_na)]

X_test = features2
y_test = df2['label']
y_test[y_test == 'murmur'] = 1
y_test[y_test == 'normal'] = 0
y_test = y_test.astype('int')

/home/silas/.local/lib/python3.8/site-packages/pandas/core/series.py:992: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._where(~key, value, inplace=True)

It can be seen that the training set contains 412 samples and our test set contains 65 samples. Moreover, on our test set, we have 34 labels murmur and 31 labels normal, so our data is evenly distributed.

Classifier ROC (Receiver Operating Characteristic Curve)

For our problem we have a total of 14 features extracted from the waveform files. We will apply a Logistic Regression, since this is a binary classification problem and analyze the model with a dummy classifier, by applying the ROC curve.

# Fitting the model
model = LogisticRegression(max_iter=1000)
model.fit(X_train, y_train)

LogisticRegression(max_iter=1000)

# Probabilities
# generate a no skill prediction (majority class)
ns_probs = [0 for _ in range(len(y_test))]

# predict probabilities
lr_probs = model.predict_proba(X_test)
# keep probabilities for the positive outcome only
lr_probs = lr_probs[:, 1]

# auc
# calculate scores
ns_auc = roc_auc_score(y_test, ns_probs)
lr_auc = roc_auc_score(y_test, lr_probs)

# summarize scores
print('No Skill: ROC AUC=%.3f' % (ns_auc))
print('Logistic: ROC AUC=%.3f' % (lr_auc))

No Skill: ROC AUC=0.500
Logistic: ROC AUC=0.855

# calculate roc curves
ns_fpr, ns_tpr, _ = roc_curve(y_test, ns_probs)
lr_fpr, lr_tpr, _ = roc_curve(y_test, lr_probs)

# plot the roc curve for the model
plt.plot(ns_fpr, ns_tpr, linestyle='--', label='No Skill')
plt.plot(lr_fpr, lr_tpr, marker='.', label='Logistic')

# axis labels
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
# show the legend
plt.legend()
# show the plot
plt.show()

# Predictions
model = LogisticRegression(max_iter=1000)
model.fit(X_train, y_train)

probs = model.predict_proba(X_test)
probs = probs[:,1]
predictions = probs >= median(np.sort(probs))
percent_score = accuracy_score(y_test, predictions)
recall = recall_score(y_test, predictions)

print(percent_score)
print(recall)

0.796875
0.7878787878787878

As we can see, our model has an accuracy of 0.797. It also has a recall of 0.788, which in this case is more important than accuracy, since we need to detect abnormalities, in case of their real occurrence.

This was just a small simple study case of time series waveforms and we could still improve more the feature engineering and model performance, by removing noise, normalizing data and extracting its variation, instead of using the absolute value.