Esercitazioni Matlab per bioingegneria

Matlab

Caricare il file

close all

clear

clc

ecgData = load('ecgClean.txt');

Plottare il file

figure;

hold on

plot(ecgData);

xlabel('Campioni');

ylabel('ECG Ampiezza [mV]');

title('Segnale ECG in campioni');

grid on

Lunghezza del tracciato in n°campioni e in tempo

n = length(ecgData)

Fs = 200 % [Hz]

timeLength = n / Fs

Frequenza max osservabile e risoluzione in frequenza

fMax = fs / 2; % frequenza max osservabile (frequenza di Nyquist)

fr = fs / n; % risoluzione in frequenza [Hz]

FFT segnale con e senza media

MediaSegnale = mean(ecgData); % calcola la media del segnale

ecg_con_media = ecgData; % segnale con media

ecg_senza_media = ecgData - MediaSegnale; % rimuove la media dal segnale

FFT_con_media = fft(ecg_con_media); % calcola le FFT

FFT_senza_media = fft(ecg_senza_media);

Ricampionamento del segnale

fsOriginal = 200; %fs originale [Hz]

fsNew = 100; % fs nuova [Hz]

ecg100 = resample(ecg_senza_media, fsNew, fsOriginal); % ricampionamento

Plot

Asse temporale

t = (0:n-1) / fs; % asse temporale

durata = 10; % secondi da visualizzare

campioni10s = durata * fs; % numero campioni in 10s

figure;

subplot(3,1,1); % 1° riga: 10s di ECG

plot(t(1:campioni10s), ecg_con_media(1:campioni10s), 'b');

hold on;

plot(t(1:campioni10s), ecg_senza_media(1:campioni10s), 'r'); % grafico

sovrapposto

xlabel('Tempo [s]');

ylabel('ECG [mV]');

title('10 s di ECG: con media (blu) e senza media (rosso)');

legend('Con media','Senza media');

grid on;

Modulo dello spettro

f = (0:n-1) * (fs / n); % asse delle frequenze

subplot(3,1,2); % 2° riga: modulo dello spettro

plot(f, abs(FFT_con_media), 'b');

hold on;

plot(f, abs(FFT_senza_media), 'r');

xlim([0 fMax]);

xlabel('Frequenza [Hz]');

ylabel('|ECG| [mV]');

title('Modulo dello spettro');

legend('Con media','Senza media');

grid on;

Fase dello spettro

subplot(3,1,3); % 3° riga: fase delle frequenze

plot(f, angle(FFT_con_media), 'b');

hold on;

plot(f, angle(FFT_senza_media), 'r');

xlim([0 fMax]);

xlabel('Frequenza [Hz]');

ylabel('Fase [rad]');

title('Fase dello spettro');

legend('Con media','Senza media');

grid on;

Filtro butterworth

Filtro passa alto

ordine = 3; % ordine del filtro butter perché

la funzione filtfilt è bidirezionale

fc = 0.4; % frequenza di taglio scelta

Wn = fc / (fs / 2); % frequenza di taglio

normalizzata

[b, a] = butter(ordine, Wn, 'high'); % filtro passa-alto butterworth

ecg_filtrato = filtfilt(b, a, ecg_senza_media); % applicazione del filtro

Filtro passa basso

ordine = 3; % ordine del filtro butter per alta

frequenza

fc = 40; % frequenza di taglio scelta per alta

frequenza

Wn = fc / (fs / 2); % frequenza di taglio normalizzata

[b, a] = butter(ordine, Wn, 'low'); % filtro passa-basso butterworth

ecg_filtrato = filtfilt(b, a, ecg_senza_media); % applicazione del filtro

Filtro passa banda

ordine = 3;

Wn = [0.4 40] / (fs/2);

[b,a] = butter(ordine, Wn, 'bandpass');

ecg_filtrato = filtfilt(b, a, ecg_senza_media);

Esercitazione 1

• Caricare plottare il file ecgClean (fs=200 Hz). Qual è l’unità di misura del segnale?

• Quanto è lungo il tracciato in tempo e in numero di campioni?

• Plottare, nella stessa figura (uno sotto l’altro), il segnale in funzione del numero dei

campioni e del

tempo.

• Plottare, nella stessa figura (una colonna e tre righe) 3 finestre di ecg in funzione del tempo

da 15 s

ognuna, localizzate all’inizio del tracciato, dopo 25 s e alla fine del tracciato.

• Sapendo che ogni battito cardiaco è lungo 750 ms, plottare 6 battiti in funzione del numero

di campioni,

il primo nero, il secondo verde, il terzo rosso, il quarto blu, il quinto marrone e il sesto rosa.

NB- tutte le figure devono mostrare le unità di misura negli assi e il titolo.

%%%%%%%%%%%%%%%%ecgClean%%%%%%%%%%%%%%%%%%

%% caricare e plottare ecgClean %%

close all

clear

clc

ecgData = load('ecgClean.txt');

figure;

hold on

plot(ecgData);

xlabel('Campioni');

ylabel('ECG Ampiezza [mV]');

title('Segnale ECG in campioni');

grid on

%% lunghezza del tracciato in tempo e in numero di campioni %%

numSamples = length(ecgData)

Fs = 200 % [Hz]

timeLength = numSamples / Fs

%% plottare il segnale in funzione del numero dei campioni e del tempo %%

t = (0:numSamples-1) / Fs; % vettore del tempo [s] per creazione

dell'asse temporale

figure;

subplot(2,1,1); % prima riga della figura

plot(1:numSamples, ecgData);

xlabel('Campione n°');

ylabel('Ampiezza [mV]');

title('Segnale ECG vs numero di campioni');

grid on;

subplot(2,1,2); % seconda riga della figura

plot(t, ecgData);

xlabel('Tempo [s]');

ylabel('Ampiezza [mV]');

title('Segnale ECG vs tempo');

grid on;

%% plottare 3 finestre di ECG da 15s l'una %%

windowLength = 15; % lunghezza di ogni finestra in [s]

numWindows = windowLength * Fs; % numero finestre da 15s

%Prima finestra: inizio del tracciato

start1 = 1;

stop1 = start1 + numWindows - 1;

%Seconda finestra: dopo 25 s

start2 = 25 * Fs+ 1;

stop2 = start2 + numWindows - 1;

%Terza finestra: alla fine del tracciato

stop3 = numSamples;

start3 = stop3 - numWindows + 1;

%creazione figura 1x3

figure;

subplot(1, 3, 1);

plot(t(start1:stop1), ecgData(start1:stop1));

xlabel('Tempo [s]');

ylabel('Ampiezza [mV]');

title('Inizio tracciato');

grid on;

subplot(1, 3, 2);

plot(t(start2:stop2), ecgData(start2:stop2));

xlabel('Tempo [s]');

ylabel('Ampiezza [mv]');

title('Dopo 25 s');

grid on;

subplot(1, 3, 3);

plot(t(start3:stop3), ecgData(start3:stop3));

xlabel('Tempo [s]');

ylabel('Ampiezza [mv]');

title('Fine tracciato');

grid on;

%% plottare 6 battiti in funzione del numero di campioni %%

beatSamples = round(0.75 * Fs); % campioni per battito (150 se Fc=200)

numBeats = 6;

totalSamples = beatSamples * numBeats;

% estrai i 6 battiti consecutivi

b1 = ecgData(1:beatSamples);

b2 = ecgData(beatSamples+1 : 2*beatSamples);

b3 = ecgData(2*beatSamples+1 : 3*beatSamples);

b4 = ecgData(3*beatSamples+1 : 4*beatSamples);

b5 = ecgData(4*beatSamples+1 : 5*beatSamples);

b6 = ecgData(5*beatSamples+1 : 6*beatSamples);

% figure con 6 segmenti consecutivi, colori diversi

figure;

hold on;

plot(1:beatSamples, b1, 'k', 'LineWidth', 1.2); %

nero

plot(beatSamples+1:2*beatSamples, b2, 'g', 'LineWidth', 1.2); %

verde

plot(2*beatSamples+1:3*beatSamples, b3, 'r', 'LineWidth', 1.2); %

rosso

plot(3*beatSamples+1:4*beatSamples, b4, 'b', 'LineWidth', 1.2); %

blu

plot(4*beatSamples+1:5*beatSamples, b5, 'Color', [0.6 0.3 0], 'LineWidth',

% marrone

1.2);

plot(5*beatSamples+1:6*beatSamples, b6, 'Color', [1 0.4 0.7], 'LineWidth',

1.2); % rosa

xlabel('Campione n°');

ylabel('Ampiezza [mV]');

title('6 battiti cardiaci consecutivi (0.75 s ciascuno)');

legend('Battito 1','Battito 2','Battito 3','Battito 4','Battito 5','Battito

6');

grid on;

Esercitazione 2

• Caricare ecgClean (fs=200 Hz).

• Indicare frequenza massima osservabile e risoluzione in frequenza.

• Calcolarne la FFT del segnale con e senza media.

• Considerare una figura 3x1 e plottare :

• 1 riga della figura: 10 sec di ECG con la media in blu, e sovrapposto gli stessi 10 s del

segnale a

media nulla

• 2 e 3 riga: plottare il modulo e la fase dello spettro del segnale con (in blue) e senza (in

rosso)

media rispettivamente.

%%%%% ecgClean %%%%%

%% caricare il file

close all

clear

clc

ecgData = load('ecgClean.txt');

%% frequenza max osservabile e risoluzione in frequenza

fs = 200; % frequenza di campionamento [Hz]

n = length(ecgData); % numero di campioni

fMax = fs / 2; % frequenza max osservabile (frequenza di Nyquist)

fr = fs / n; % risoluzio