How to band pass filter the EMG signal in real time?

%#######################################
% In signal processing, the power spectrum of a continuous time signal describes 
% the distribution of power into frequency components. composing that signal. According to Fourier analysis, 
% any physical signal can be decomposed into a number of discrete frequencies, 
% or a spectrum of frequencies over a continuous range.
% Get the full wave rectification of the EMG signal by measuring it's
% absolute value.
% ######################################
clc;
clear;
close all;
% Parameters
sampling_frequency = 1000;              % Sampling frequency in Hz
duration = 5;                           % Duration of signal acquisition in seconds
num_samples = sampling_frequency * duration; % Sampling_frequency = (No. of Samples) / Time Interval
% Initialize variables
time = (0:num_samples-1)/sampling_frequency;
frequency = linspace(0, sampling_frequency/2, num_samples/2);
x = zeros(1, num_samples);
a = arduino();                          % Check port from device manager and pin is A0
% Create figure for real-time EMG signal
figure('Name', 'Real-Time EMG Signal Visualization');
h_real_time = plot(time, x);
grid on;
title('Real-Time EMG Signal');
xlabel('Time (s)');
ylabel('Voltage (V)');
xlim([0, duration]);        % Adjust xlim as needed
ylim([0, 3]);               % Adjust ylim as needed
% Create figure for spectrogram
figure('Name', 'Real-Time Spectrogram EMG Signal Visualization');
h_spectrogram = imagesc('XData', [], 'YData', [], 'CData', []);
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Spectrogram of EMG Signal');
colorbar;
% Pre-allocate spectrogram variables
Nspec = 256;  % Spectrogram window size
overlap = round(Nspec * 0.85);  % 85% overlap
% Filter parameters
filter_order = 4;                    % Filter order
bandpass_low_cutoff = 10;            % Bandpass low cutoff frequency in Hz
bandpass_high_cutoff = 500;          % Bandpass high cutoff frequency in Hz
% Design Butterworth bandpass filter
assert(filter_order > 0 && round(filter_order) == filter_order, 'Filter order must be a positive integer');
Wn = [bandpass_low_cutoff, bandpass_high_cutoff] / (sampling_frequency / 2);
[b, a] = butter(filter_order, Wn, 'bandpass');
% Start loop to acquire signal
for i = 1:num_samples
    % Read voltage
    voltage = abs(readVoltage(a, 'A0'));
    
    % Store voltage
    x(i) = voltage;
    
    % Update real-time plot
    set(h_real_time,'YData',abs(x));
    set(h_real_time,'XData',time);
    drawnow;
    
    % Apply bandpass filtering
    if i > filter_order
        x_filtered = filtfilt(b, a, abs(x(1:i)));
        
        % Calculate spectrogram for the current signal
        if i > Nspec
            [S, F, T] = spectrogram(x_filtered, hamming(Nspec), overlap, Nspec, sampling_frequency, 'yaxis');
            % Update the spectrogram
            set(h_spectrogram, 'XData', T, 'YData', F, 'CData', 10*log10(abs(S)));
        end
    end
    
    drawnow;
end

How do I band pass this signal? And is it right to rectify it after I bandpass filter it?

[Iro Liontou]

That's really helpful. Thank you very much!!!

Answer 1

[Balavignesh]

Hi Iro,

In signal processing, especially when dealing with Electromyography (EMG) signals, the order of operations can significantly affect both the outcome and the interpretability of the processed signal. Bandpass filtering is designed to remove noise and unwanted frequencies outside the band of interest. This typically involves eliminating both very low-frequency motion artifacts and very high-frequency electrical noise. Performing this step before rectification helps to ensure that the signal's frequency content is limited to the physiologically relevant range, thereby improving the Signal-to-Noise Ratio (SNR) and making subsequent analysis more meaningful.

Here is a rough approach for implementation:

• Read the raw voltage without immediately applying rectification.
• Apply the bandpass filter to the raw signal or the most recent set of readings.
• Rectify the filtered signal by taking the absolute value after filtering.

The following code snippet may help you achieve this:

for i = 1:num_samples
    % Read raw voltage
    raw_voltage = readVoltage(a, 'A0');  % No immediate rectification
    
    % Store raw voltage
    x(i) = raw_voltage;
    
    % Apply bandpass filtering first
    if i > filter_order
        x_filtered = filtfilt(b, a, x(1:i));
        
        % Rectify after filtering
        x_rectified = abs(x_filtered);
        
        % Update real-time plot with rectified signal
        set(h_real_time,'YData',x_rectified);
        set(h_real_time,'XData',time);
        drawnow;
        
        % Calculate spectrogram for the rectified signal
        if i > Nspec
            [S, F, T] = spectrogram(x_rectified, hamming(Nspec), overlap, Nspec, sampling_frequency, 'yaxis');
            % Update the spectrogram
            set(h_spectrogram, 'XData', T, 'YData', F, 'CData', 10*log10(abs(S)));
        end
    end
    
    drawnow;
end

Kindly have a look at the following documentation links to have more information on:

• bandpass filters: https://www.mathworks.com/help/signal/ref/bandpass.html

Hope that helps!

Balavignesh