Bio Medical Signal Processing
(EC4061D)
ASSIGNMENT 2
Bachelor of Technology
in
Electronics and Communication Engineering
Muhammed Jasim M
B220412EC, EC03
Department of Electronics and Communication
Engineering
National Institute of Technology Calicut Calicut, Kerala, India
– 673 601
Winter Semester 2025
ASSIGNMENT 2 .............................................................................................................................. 1
Department of Electronics and Communication Engineering...................................................... 1
Questions ...................................................................................................................................... 2
Question 1: Detect QRS Complex Using Pan-Tompkins Algorithm and Find Heart Rate .................. 3
Theory Explanation .................................................................................................................... 3
Steps in the Pan-Tompkins Algorithm ........................................................................................ 3
Heart Rate Calculation ............................................................................................................... 4
MATLAB Implementation........................................................................................................... 4
Question 2: Detect α Wave Rhythm and Spikes in EEG Using eeg1-p4.dat ..................................... 6
2a. Detect α Wave Rhythm Present in EEG Waveform ................................................................... 6
Theory ....................................................................................................................................... 6
Steps ......................................................................................................................................... 6
MATLAB Code for Alpha Rhythm Detection ............................................................................... 7
2b. Detect Spikes in EEG Using Matched Filter ............................................................................... 9
Theory ....................................................................................................................................... 9
Steps ......................................................................................................................................... 9
MATLAB Code for Spike Detection ............................................................................................. 9
Question 3: Find RR Interval from PPG Signal............................................................................... 11
Theory ..................................................................................................................................... 11
Steps to Find RR Interval .......................................................................................................... 11
MATLAB Code for RR Interval Detection .................................................................................. 12
Questions
1. Detect QRS complex by using Pan-Tompkins Algorithm. Find
Heart Rate.
2. a) Detect α wave rhythm present in EEG waveform:
ACF: In PSD, pass signal.
PSD: Through PSD function, highest peak indicates α rhythm.
b) Detect spikes in EEG waveform using matched filter.
3. Find RR from PPG signal.
Question 1: Detect QRS Complex Using Pan-Tompkins
Algorithm and Find Heart Rate
Theory Explanation
The QRS complex in an ECG signal represents the electrical activity of ventricular
depolarization. Detecting the QRS complex is essential for calculating heart rate.
The Pan-Tompkins algorithm is a robust method designed to detect QRS
complexes accurately in noisy ECG signals.
Steps in the Pan-Tompkins Algorithm
1. Preprocessing:
a. Band-pass filtering: Removes baseline wander and high-frequency
noise while preserving QRS features.
b. Typical range: 5–15 Hz.
2. Differentiation:
a. Highlights rapid changes in the signal, emphasizing the steep slopes
of the QRS complex.
3. Squaring:
a. Non-linear operation that amplifies large differences (QRS peaks)
and suppresses smaller ones.
4. Integration:
a. A moving average filter smooths the signal and highlights regions of
interest corresponding to QRS complexes.
5. Thresholding:
a. Adaptive thresholds are applied to detect peaks corresponding to R-
waves (highest point in the QRS complex).
Heart Rate Calculation
Once R-peaks are detected, heart rate is computed using the formula:
HR=60RR_intervalHR=RR_interval60
Where RR_intervalRR_interval is the time difference between consecutive R-peaks,
measured in seconds.
MATLAB Implementation
Below is the MATLAB code to detect QRS complexes using the Pan-Tompkins
algorithm and calculate heart rate:
% Load ECG signal
load('ecg_fs_500.mat');
% Get variable name (assuming it's the first variable in the file)
varNames = who('-file', 'ecg_fs_500.mat');
ecgVarName = varNames{1};
ecg_signal = eval(ecgVarName);
% Ensure ecg_signal is a row vector
if size(ecg_signal, 1) > size(ecg_signal, 2)
ecg_signal = ecg_signal';
end
fs = 500; % Sampling frequency
% Band-pass filter
lowcut = 5; highcut = 15;
[b, a] = butter(2, [lowcut, highcut]/(fs/2), 'bandpass');
filtered_ecg = filtfilt(b, a, ecg_signal);
% Differentiation
diff_ecg = diff(filtered_ecg);
% Squaring
squared_ecg = diff_ecg.^2;
% Moving average integration
window_size = round(0.150 * fs);
integrated_ecg = movmean(squared_ecg, window_size);
% Thresholding and peak detection
threshold = 0.6 * max(integrated_ecg);
[R_peaks, locs] = findpeaks(integrated_ecg, 'MinPeakHeight', threshold);
% Calculate heart rate
RR_intervals = diff(locs) / fs;
HR = 60 ./ RR_intervals;
average_HR = mean(HR);
% Display results
disp(['Average Heart Rate: ', num2str(average_HR), ' bpm']);
% Plot results
figure;
subplot(3,1,1);
plot(ecg_signal);
title('Original ECG Signal');
xlabel('Samples');
ylabel('Amplitude');
subplot(3,1,2);
plot(filtered_ecg);
title('Filtered ECG Signal');
xlabel('Samples');
ylabel('Amplitude');
subplot(3,1,3);
plot(integrated_ecg);
hold on;
plot(locs, integrated_ecg(locs), 'ro');
title('Integrated Signal with Detected R-Peaks');
xlabel('Samples');
ylabel('Amplitude');
Output Figure:
Average Heart Rate: 1379.3958 bpm
Question 2: Detect α Wave Rhythm and Spikes in EEG
Using eeg1-p4.dat
2a. Detect α Wave Rhythm Present in EEG Waveform
Theory
Alpha waves are brain signals with a frequency range of 8–13 Hz, typically
observed during relaxed states or when the eyes are closed. To detect alpha
rhythms, we analyze the EEG signal's Power Spectral Density (PSD) using methods
like Welch's method or FFT. Peaks within the alpha band (8–13 Hz) indicate the
presence of alpha activity.
Steps
1. Autocorrelation Function (ACF):
a. ACF analyzes periodicity in the signal, helping identify dominant
frequencies.
b. It provides insight into the rhythmic nature of the alpha waves.
2. Power Spectral Density (PSD):
a. PSD is computed using methods like Welch's method to measure
signal power across frequencies.
b. Peaks in the alpha band (8–13 Hz) indicate alpha rhythm.
MATLAB Code for Alpha Rhythm Detection
% Load EEG signal from 'eeg1-p4.dat'
eeg_signal = load('eeg1-p4.dat');
% Load data from file
fs = 250; % Sampling frequency (modify based on your data)
% Step 1: Compute Autocorrelation Function (ACF)
acf = xcorr(eeg_signal, 'coeff');
% Normalized autocorrelation
acf_lags = (-length(eeg_signal)+1:length(eeg_signal)-1)/fs;
% Plot ACF
figure;
plot(acf_lags, acf);
title('Autocorrelation Function of EEG Signal');
xlabel('Lag (seconds)');
ylabel('Correlation');
% Step 2: Compute Power Spectral Density (PSD)
[psd_values, freq] = pwelch(eeg_signal, hamming(512), 256, 1024, fs);
% Extract alpha band (8–13 Hz)
alpha_band = freq(freq >= 8 & freq <= 13);
alpha_psd = psd_values(freq >= 8 & freq <= 13);
% Plot PSD
figure;
plot(freq, psd_values);
hold on;
plot(alpha_band, alpha_psd, 'r', 'LineWidth', 1.5);
% Highlight alpha band
title('Power Spectral Density of EEG Signal');
xlabel('Frequency (Hz)');
ylabel('Power');
legend('PSD', 'Alpha Band');
% Display Alpha Band Power
disp(['Alpha Band Power: ', num2str(sum(alpha_psd))]);
Output Figure:
Alpha Band Power: 198032.73
2b. Detect Spikes in EEG Using Matched Filter
Theory
A matched filter is a signal processing technique used to detect known patterns
(e.g., spikes) in noisy signals. The filter is designed to maximize the correlation
between the input signal and a predefined spike template.
Steps
1. Define a spike template based on prior knowledge or observed spike
patterns.
2. Apply the matched filter by convolving the EEG signal with the reverse of
the spike template.
3. Identify spikes as peaks exceeding a threshold.
MATLAB Code for Spike Detection
% Load EEG signal from 'eeg1-p4.dat'
eeg_signal = load('eeg1-p4.dat');
% Load data from file
% Define spike template manually or load it from a file
spike_template = [0, -1, -2, -3, -2, -1, 0];
% Example spike template
fs = 250; % Sampling frequency
% Step 1: Create matched filter using spike template
matched_filter = flip(spike_template); % Reverse the spike template
% Step 2: Apply matched filter to EEG signal
filtered_signal = conv(eeg_signal, matched_filter, 'same');
% Step 3: Thresholding for spike detection
threshold = 0.7 * max(filtered_signal); % Set threshold based on filtered
signal amplitude
spike_indices = find(filtered_signal > threshold);
% Fix: Convert spike_indices to a row vector for display
spike_indices_row = reshape(spike_indices, 1, []); % Ensure it's a row vector
% Display detected spikes
disp(['Detected spikes at indices: ', num2str(spike_indices_row)]);
% Plot results
figure;
subplot(2,1,1);
plot(eeg_signal);
title('Original EEG Signal');
xlabel('Samples');
ylabel('Amplitude');
subplot(2,1,2);
plot(filtered_signal);
hold on;
plot(spike_indices, filtered_signal(spike_indices), 'ro'); % Mark detected
spikes
title('Filtered Signal with Detected Spikes');
xlabel('Samples');
ylabel('Amplitude');
grid on;
Output Figure:
Detected spikes at indices: 90 101 102 114 115 116 129 130 181
182 183 184 194 195 196 222 229 230 231 232 239 240 241
242 243 252 253 443 444 622 623 703 704 705
Question 3: Find RR Interval from PPG Signal
Theory
Photoplethysmogram (PPG) signals measure blood volume changes in the
microvascular bed of tissue, providing information about heartbeats. The RR
interval represents the time between consecutive peaks in the PPG waveform,
corresponding to heartbeats. RR intervals are crucial for calculating heart rate
variability and other cardiovascular metrics.
Steps to Find RR Interval
1. Preprocessing:
a. Filter the PPG signal to remove noise and artifacts using a band-pass
filter.
b. Typical frequency range for heartbeats in PPG is 0.5–4 Hz.
2. Peak Detection:
a. Identify peaks in the PPG signal that correspond to heartbeats.
b. Use MATLAB's findpeaks function with appropriate thresholds.
3. RR Interval Calculation:
a. Compute the time difference between consecutive peaks (R-peaks).
b. Convert sample indices into seconds using the sampling frequency.
MATLAB Code for RR Interval Detection
% Load PPG signal from 'PPG_fs_250.mat'
data = load('PPG_fs_250.mat'); % Load data from file
ppg_signal = data.x; % Extract PPG signal (update variable name if needed)
fs = 250; % Sampling frequency (as indicated by filename)
% Step 1: Band-pass filtering
lowcut = 0.5; % Lower cutoff frequency (Hz)
highcut = 4; % Upper cutoff frequency (Hz)
[b, a] = butter(2, [lowcut, highcut]/(fs/2), 'bandpass'); % Design band-pass
filter
filtered_ppg = filtfilt(b, a, ppg_signal); % Apply zero-phase filtering
% Step 2: Peak detection
[peaks, locs] = findpeaks(filtered_ppg, 'MinPeakHeight', mean(filtered_ppg),
'MinPeakDistance', round(0.6 * fs));
% Step 3: Calculate RR intervals
RR_intervals = diff(locs) / fs; % Convert sample differences to seconds
% Display results
disp('RR Intervals (seconds):');
fprintf('%f\n', RR_intervals); % Print each RR interval on a new line
disp(['Mean RR Interval: ', num2str(mean(RR_intervals)), ' seconds']);
% Plot results
figure;
subplot(2,1,1);
plot(ppg_signal);
title('Original PPG Signal');
xlabel('Samples');
ylabel('Amplitude');
subplot(2,1,2);
plot(filtered_ppg);
hold on;
plot(locs, filtered_ppg(locs), 'ro'); % Mark detected peaks
title('Filtered PPG Signal with Detected Peaks');
xlabel('Samples');
ylabel('Amplitude');
Output Figure:
