How can I draw 117 normal 292 abnormal heart sounds?

I have a total of 409 heart sounds, 117 of which are normal and 292 of which are abnormal. I have prepared separate codes normally and abnormally, how should I go about drawing each sound?

The audio files were not added because they were too large, this is the example for which I found the code: https://www.mathworks.com/matlabcentral/fileexchange/65286-heart-sound-classifier

Database: https://archive.physionet.org/pn3/challenge/2016/

For example;

% fs:2000
[PCG_normal, fs] = audioread('a0011.wav');
p_normal = audioplayer(PCG_normal, fs);
play(p_normal, [1 (get(p_normal, 'SampleRate') * 3)]);
% Plot the sound waveform
plot(PCG_normal(1:fs*3))


[PCG_abnormal, fs] = audioread('a0001.wav');
p_abnormal = audioplayer(PCG_abnormal, fs);
play(p_abnormal, [1 (get(p_abnormal, 'SampleRate') * 3)]);
%        Plot the sound waveform
plot(PCG_abnormal(1:fs*3))

Answer 1

[Star Strider]

Drawing them depends on the result you want. Heart sounds are best displayed as time-frequency plots, so I would plot them using the pspectrum function with the 'spectrogram' option. (Ideally, they should be plotted along with the corresponding Lead II EKG trace.)

Plotting 409 of them might be something of a challenge, since tiledlayout and subplot would create axes that are too small to easily visualise.

EDIT — (2 Apr 2023 at 17:32)

This emphasizes the time-frequency characteristics of the phonocardiogram signals —

Uz1 = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1342909/heart%20sound.zip');
Uz2 = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1342914/records.zip');

[PCG_normal, fsn] = audioread(Uz1{2});
p_normal = audioplayer(PCG_normal, fsn);
Ln = size(PCG_normal,1);
tn = linspace(0, Ln-1, Ln)/fsn;
play(p_normal, [1 (get(p_normal, 'SampleRate') * 3)]);
% Plot the sound waveform
figure
plot(PCG_normal(1:fsn*3))

[PCG_abnormal, fsa] = audioread(Uz1{1});
p_abnormal = audioplayer(PCG_abnormal, fsa);
La = size(PCG_abnormal,1);
ta = linspace(0, La-1, La)/fsa;
play(p_abnormal, [1 (get(p_abnormal, 'SampleRate') * 3)]);
%        Plot the sound waveform
figure
plot(PCG_abnormal(1:fsa*3))

figure
tiledlayout(2,2)
nexttile
pspectrum(PCG_normal, fsn,'spectrogram')
colormap(turbo)
nexttile
pspectrum(PCG_abnormal, fsa,'spectrogram')
colormap(turbo)
nexttile
plot(tn,PCG_normal)
xlabel('t')
ylabel('Amplitude')
title('PCG Normal')
xlim([0 max(tn)])
grid
nexttile
plot(ta,PCG_abnormal)
xlabel('t')
ylabel('Amplitude')
title('PCG Abnormal')
xlim([0 max(ta)])
grid

.

[Cey]

Thank you very much, let's say that I manually drew each sound one by one. If I want to do the filtering process later, how do you think I can do it in this code? Or how do I do automatic segmentation? @Star Strider

data=PCG_normal;
data1=PCG_abnormal;

L = numel(data);
Fs = 2000;              
Fn = Fs/2;
t = linspace(0, L-1, L)/Fs;
X=t(:,1:10000);
L1 = numel(data1);
Fs1 = 2000;              
Fn = Fs1/2;
t1 = linspace(0, L1-1, L1)/Fs1;
X1= t1(:,1:10000);

% % Frequency for 20 400 
Wn=[0.02 0.4];
[B,A] = butter(4,Wn,'bandpass');
Normal= filtfilt(B,A,data);

% % Frequency for 20 800 
Wn1=[0.02 0.8];
[B1,A1] = butter(4,Wn1,'bandpass');
Abnormal = filtfilt(B1,A1,data1);

figure  
subplot(2,1,1)
plot(X1,Abnormal(1:10000,:))
grid
title('butterworth bandpass abnormal')
subplot(2,1,2)
plot(X,Normal(1:10000,:))
grid
title('butterworth bandpass normal')

[Star Strider]

My pleasure!

There is nothing to be gained by filtering them, especially with those filter passbands (that appear to be normalised, so equate to 20-400 and 20-800 Hz, respectively). As is apparent in the spectrogram plots, there is little appreciable noise, and much of the most relevant information would be filtered out with those filters (for example the high-frequency component at about 850 Hz in the ‘abnormal’ record). I wouild also do spectrogram plots on several of the ’abnormal’ records to be certain that the spectral characteristic are relatively constant. Systolic ejection murumurs have a high-frequency component, however the low-frequency diastolic rumble present in atrial septal defect phonocardiograms, in addition to the widely-split second heart sound, would be worth additional consideration if those records are present in the ‘abnormal’ set.

So for signal processing purposes, I would just leave them as they are.

If you want to classify them or do some other sort of feature segmentation, I would choose two (or more) filters with slightly overlapping passbands, perhaps 20-500 Hz and 500-950 Hz, and then compute the relative power (the sum of the squared amplitudes) in those filter outputs. That should allow you to differentiate normal from abnormal phonocardiogram records (although phonocardiography is largely considered to be obsolete currently, having been replaced by ultrasound, and particularly point-of-care ultrasound, or POCUS).

The most computationally efficient filters are elliptic filters, so I would also replace the Butterworth filters with them. Start with ellipord, then ellip then zp2sos.

.

[Cey]

Thank you very much for your help. I'm going to segment the heart sounds in my hand and then classify them. I want to do the work automatically because there is a lot of sound. I'll also try the elliptic filtering you said for the filter. @Star Strider

[Star Strider]

Probably the best way to do this would be to create four or five filters that evenly divide the spectrum. Then calculate the total power in each band for each signal, and use those for features. Heart sounds are functions of both time and frequency, however that would require some sort of image pattern recognition on the spectrogram plots.

I suspect that this is likely an area of research, and a PubMed search will discover previous efforts.

[Cey]

Thank you again for your help. I divided the sounds into abnormal and normal, but as you said, it will be difficult for me to filter or segment them all together. Obviously, it is impossible to deal with 409 total heart sounds one by one, how do you think I can follow? @Star Strider

[Star Strider]

I’m not sure what to suggest. As I mentioned, pohoncardiography is essentially obsolete, having been supplanted by POCUS. A PubMed search is likely your best option.

I am apparently not able to solve your problem, so I will delete my Answer in a few hours.

[Cey]

First of all, thank you very much for your help. I will have another question, with your permission, how can I make a normal filtered heart sound in 6 or other different segments for each s1 and s2 sound so that i can do feature extraction on individual segments. @Star Strider

[Star Strider]

My pleasure.

Perhaps something like this —

LD = load('normal.mat')

LD = struct with fields:
    normal: [10000×1 double]

normal = LD.normal;
Fs = 2E3;
L = numel(normal);
t = linspace(0, L-1, L)/Fs;

figure
tiledlayout(2,1)
nexttile
plot(t,normal)
grid
nexttile
pspectrum(normal, Fs,'spectrogram')
colormap(turbo)

Lf = 7;                                                                 % Number Of Filters+1
PBv = linspace(0.01, 0.975, Lf);                                        % Passband Vector For 6 Filters
for k = 1:Lf-1
    Wp = [PBv(k) PBv(k+1)];                                             % Passband Frequency (Normalised)
    Ws = [0.975 1.025].*Wp;                                             % Stopband Frequency (Normalised)
    Rp =  1;                                                            % Passband Ripple
    Rs = 60;                                                            % Passband Ripple (Attenuation)
    [n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
    [z,p,kf] = ellip(n,Rp,Rs,Wp);                                       % Elliptic Filter Design: Zero-Pole-Gain 
    [sos{k},g(k)] = zp2sos(z,p,kf);                                     % Second-Order Section For Stability
    
    figure
    freqz(sos{k}, 2^16, Fs)                                             % Filter Bode Plot
end

figure
tiledlayout(fix((Lf-1)/2),2)
for k = 1:Lf-1
    y(:,k) = filtfilt(sos{k},g(k),normal);                               % Filter Signal
    nexttile
    plot(t,y(:,k))
    grid
    xlim([min(t) max(t)])
    title(sprintf('%.2f - %.2f Hz',PBv(k),PBv(k+1)))
end

Displaying the freqz results is not absolutely necessary. I did it here to demonstrate the different segmental filters.

(It would be helpful to have the sampling frequency or time vector as one of the variables included in the .mat files. Guessing that it‘ 2k here.)

.

[Cey]

thank you very much just as what i wanted

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.