AutomatedLungSoundSegmentation/labelInspExp.m at master · titusaj/AutomatedLungSoundSegmentation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
clear all
close all
clc

files = dir(fullfile(pwd, '*.txt'));
cycleCount = 1; % This is numfiles * the cycle count (based on each file thats read into memory)


for i = 8:length(files)
i
%Intialize and clear variables each time
rawWholeSignal = [];
rawTime = [];
cycleStart = [];
cycleEnd = [];
cracklePresent = [];
wheezePresent = [];
numCycles = 0 ;


eventFilename = files(i).name;

temp = strsplit(eventFilename,'.');
recordingLabel = temp{1};


temp2 = strsplit(temp{1},'_');
hardware = temp2{5}

patientNum = eventFilename(1:3);


wavFilename = strcat(temp2{1},'_',temp2{2},'_',temp2{3},'_',temp2{4},'_',temp2{5},'.wav');

%% This is for reading the segmenetation events
audioLabelData =  textread(eventFilename); %Read the text file into workspace
cycleStart = audioLabelData(:,1);
cycleEnd = audioLabelData(:,2);
cracklePresent = audioLabelData(:,3);
wheezePresent = audioLabelData(:,4);
numCycles = length(cycleStart);


if isfile(wavFilename) && sum(wheezePresent)  > 0 || sum(cracklePresent) > 0   % right now just looking the normal fiels
       %% This is to calculate information about cycles
         for j = 1:numCycles
             allCycleLengths(j) = cycleEnd(j) - cycleStart(j);
             cycleCount = cycleCount +1;
         end

         halfMedianCycleLength = median(allCycleLengths)/2;

        %% Read the raw signal
        [rawWholeSignal,fs] = audioread(wavFilename ); %Read the signal in if applicable
        dt = 1/fs;
        Norig = length(rawWholeSignal);
        rawTime = 0:dt:(Norig*dt)-dt;

        %% Extract the PCG features
%      [PCG_Features, featuresFs] = getSpringerPCGFeatures(rawWholeSignal, fs);

        %% Filter the signal
        filter_out = bandpass(rawWholeSignal,[200 1000],fs);

        % norm the signal magnitude
        normfilterOut = filter_out - min(filter_out(:));
        normfilterOut = normfilterOut ./ max(normfilterOut(:)) ;% *


        %% Segment out the cycles
        for j = 1:numCycles
            [d, indexStarts(j)] = min(abs( rawTime-round(cycleStart(j),3)));
            [d, indexEnds(j) ] = min(abs( rawTime-round(cycleEnd(j),3)));
%             groundTruthSegmentedSignal = normfilterOut(indexStart:indexEnd);
%             THRESHOLD = 2.5
%             runLTSD(fs,groundTruthSegmentedSignal,THRESHOLD )
        end
%         close all
%
%         figure(1);
%         t1 = (1:length(groundTruthSegmentedSignal))./fs;
%         plot(t1, groundTruthSegmentedSignal)

        %% Plot the signal with the ground truth
            figure('Name', 'PCG features');
            t1 = (1:length(rawWholeSignal))./fs;
            plot(t1,rawWholeSignal);
            hold on;
            %t2 = (1:length(PCG_Features))./featuresFs;
            %plot(t2,PCG_Features);


            for k=1:length(cycleStart)
                vline([cycleStart(k)],['g'])
                %%this draw a line midway adding the half median cycle
                %%count to the start of the cycle
                 %if k ~=length(cycleStart) % check not at the end
                 %    vline([cycleStart(k)+halfMedianCycleLength],['k'])
                 %end
           end

            hold on
            for k=1:length(cycleEnd)
                vline([cycleEnd(k)],['r'])
            end
            hold off


            %% Create Ground Truth Envelope bsed on the segmented cycle annotation
            groundTruthEnvolope = zeros(1, length(t1));
            for j = 1:numCycles

                if wheezePresent(j) == 1
                    groundTruthEnvolope(indexStarts(j):indexEnds(j)) = 1; %Set wheeze to 1

                elseif cracklePresent(j) == 1
                    groundTruthEnvolope(indexStarts(j):indexEnds(j)) = 2; %Set crackle to two

                elseif cracklePresent(j) == 1 && wheezePresent(j) == 1
                    groundTruthEnvolope(indexStarts(j):indexEnds(j)) = 3; %Set crackle to three

                end
            end

            hold on
            plot(t1,groundTruthEnvolope,'k');

%% Write the raw vector to a file
%         cd('rawHilberts/')
%         %csvwrite(strcat(temp{1},'_data.csv'),hilbertEnv)
%         length(PCG_Features(:,3))
%         dlmwrite(strcat(temp{1},'_data.csv'),PCG_Features(:,3)); % 30 times faster!
%         cd ..
%


end

end