-
Notifications
You must be signed in to change notification settings - Fork 2
/
BodySignals.m
150 lines (132 loc) · 4.97 KB
/
BodySignals.m
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
146
147
148
149
150
%{
The ECG, EMG, and EEG was taken from a patient.
ECG was recorded in a .csv file
EMG and EEG was recorded in a .txt file with different delimiters
Our task:
Load the data from the given files into a MATlab script and plot the data.
Then take the discrete fourier transforms of each signal then plot the
frequency Spectrums. In addition, create a digital notch filter to remove
the noise coming from mains power in the ECG signal, determine if the patient was tired
during the recording of the EMG signal based on the frequency spectrum of
the EMG signal, and determine the mental state of the patient during
recording of the EEG signal based on the frequency spectrum of the EEG
signal. Finally, break down how much of the EEG signal is comprised of the
four EEG wave components (Delta, Theta, Alpha and Beta), and display it in a bar graph.
Created by Tyler Adam Martinez
Date: 12/12/2019
%}
%% ECG Portion
%Gathering ECG data from file
ECGfile = csvread("ECGLab.csv");
ECGtime = ECGfile(:,1) + 2;
ECGsignal = ECGfile(:, end);
%frequency spectrum analysis of the RAW ECG data
ECGN = length(ECGtime);
ECGfs = 2500;
ECGf = ECGfs*(0:(ECGN/2))/ECGN;
ECGy = fft(ECGsignal);
ECGP2 = abs(ECGy/ECGN);
ECGP1 = ECGP2(1:ECGN/2+1);
ECGP1(2:end-1) = 2*ECGP1(2:end-1);
%filering out the 60Hz noise from the signal
Num = designfilt('bandstopiir','FilterOrder',2, ...
'HalfPowerFrequency1',59,'HalfPowerFrequency2',61, ...
'DesignMethod','butter','SampleRate',ECGfs);
ECGfiltered = filtfilt(Num, ECGsignal);
%frequency spectrum analysis of the Filtered ECG data
ECGN = length(ECGtime);
ECGfs = 2500;
ECGf = ECGfs*(0:(ECGN/2))/ECGN;
ECGyfiltered = fft(ECGfiltered);
ECGP2fil = abs(ECGyfiltered/ECGN);
ECGP1fil = ECGP2fil(1:ECGN/2+1);
ECGP1fil(2:end-1) = 2*ECGP1fil(2:end-1);
%% EMG Portion
%Gathering EMG data from file
EMGfile = fopen('EMG.txt', 'r');
EMGdata = fscanf(EMGfile, '%f');
%frequency spectrum analysis of the EMG data
EMGN = length(EMGdata);
EMGfs = 1000;
EMGtime = 1/(EMGfs):1/(EMGfs):30;
EMGtime = transpose(EMGtime);
EMGf = EMGfs*(0:(EMGN/2))/EMGN;
EMGy = fft(EMGdata);
EMGP2 = abs(EMGy/EMGN);
EMGP1 = EMGP2(1:EMGN/2+1);
EMGP1(2:end-1) = 2*EMGP1(2:end-1);
%Calculating Area underneath curve using Reimann Sums
IntFatigue = sum(EMGP1(1:1501));
IntNFatigue = sum(EMGP1(1:3001, 1));
Fatigue = IntFatigue/IntNFatigue;
%% EEG Portion
%Gathering EEG data from file
EEGdata = dlmread('EEGSignal.txt', '\t');
EEGsignal = EEGdata(1, :);
EEGtime = EEGdata(2, :);
%frequency spectrum analysis of the EEG data
EEGN = length(EEGtime);
EEGfs = 2500;
EEGf = EEGfs*(0:(EEGN/2))/EEGN;
EEGy = fft(EEGsignal);
EEGP2 = abs(EEGy/EEGN);
EEGP1 = EEGP2(1:EEGN/2+1);
EEGP1(2:end-1) = 2*EEGP1(2:end-1);
%% plotting the data for ECG Unfiltered
figure('Name','Electrocardiography Unfiltered','NumberTitle','off');
subplot(2, 1, 1);
plot(ECGtime, ECGsignal);
title("RAW ECG Signal"); xlabel("Seconds (s)"); ylabel("Amplitude (100 mV)");
subplot(2,1,2);
plot(ECGf, ECGP1);
xlabel('frequency (hz)'); ylabel('Amplitude'); title('Frequency Spectrum');
xlim([0 100]);
%% plotting the data for ECG Filtered
figure('Name','Electrocardiography Filtered','NumberTitle','off');
subplot(2, 1, 1);
plot(ECGtime, ECGfiltered);
title("Filtered ECG Signal"); xlabel("Seconds (s)"); ylabel("Amplitude (100 mV)");
subplot(2,1,2);
plot(ECGf, ECGP1fil);
xlabel('frequency (hz)'); ylabel('Amplitude'); title('Frequency Spectrum');
xlim([0 100]);
%% plotting the data for ECG Unfiltered vs. Filtered
figure('Name','Electrocardiographs Unfiltered Vs Filtered','NumberTitle','off');
subplot(2, 1, 1);
plot(ECGtime, ECGsignal);
title("Unfiltered ECG Signal"); xlabel("Seconds (s)"); ylabel("Amplitude (100 mV)");
subplot(2, 1, 2);
plot(ECGtime, ECGfiltered);
title("Filtered ECG Signal"); xlabel("Seconds (s)"); ylabel("Amplitude (100 mV)");
%% plotting the data for EMG
figure('Name','Electromyographs ','NumberTitle','off');
subplot(2, 1, 1);
plot(EMGtime, EMGdata);
title("EMG Signal"); xlabel("Seconds (s)"); ylabel("Amplitude (100 mV)");
subplot(2,1,2);
plot(EMGf, EMGP1);
xlabel('frequency (hz)'); ylabel('Amplitude'); title('Frequency Spectrum');
xlim([0 100]);
%% plotting the data for EEG
figure('Name','Electroencephalograhy','NumberTitle','off');
subplot(3, 1, 1);
plot(EEGtime, EEGsignal);
title("EEG Signal"); xlabel("Seconds(s)"); ylabel("Amplitude (uV)");
subplot(3,1,2);
plot(EEGf, EEGP1);
xlabel('frequency (hz)'); ylabel('Amplitude'); title('Frequency Spectrum');
xlim([0 30]);
%% Sorting the EEG signal into its four components and displaying it in a histogram
Delta = EEGP1(1:17);
Delta = sum(Delta);
Theta = EEGP1(18:33);
Theta = sum(Theta);
Alpha = EEGP1(34:49);
Alpha = sum(Alpha);
Beta = EEGP1(50:162);
Beta = sum(Beta);
EEGbarlabels = categorical({'Delta','Theta','Alpha','Beta'});
EEGbar = [Delta Theta Alpha Beta];
subplot(3,1,3);
EEGb = bar(EEGbarlabels, EEGbar);
title('EEG Signals'); ylabel('Amplitude');