**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Bandpass filter Transient Response. How to get rid of it?

40 views (last 30 days)

Show older comments

I created this bandpass filter and plotted it's impulse response:

% Bandpass filter parameters

f_nyquist = sampling_frequency / 2;

low_cutoff_frequency = 6; % Low cutoff frequency in Hz

high_cutoff_frequency = 400; % High cutoff frequency in Hz

% Design the bandpass filter

[b, a_filter] = butter(2, [low_cutoff_frequency, high_cutoff_frequency] / (f_nyquist), 'bandpass');

% Create an impulse signal (Dirac delta function)

impulse_signal = zeros(1, 1000); % Change 1000 to the desired length

impulse_signal(1) = 1; % Set the first sample to 1

% Apply bandpass filter to the impulse signal

impulse_response = filtfilt(b, a_filter, impulse_signal);

% Time vector for plotting

time = (0:length(impulse_response)-1) / sampling_frequency;

% Plot the impulse response

figure('Name', 'Impulse Response of the Bandpass Filter');

plot(time, impulse_response);

title('Impulse Response of the Bandpass Filter');

xlabel('Time (s)');

ylabel('Amplitude');

grid on;

[h, f] = freqz(b, a_filter, 1024, 'whole'); % Frequency response

figure;

plot(f/pi * (sampling_frequency/2), abs(h)); % Plot magnitude response

xlabel('Frequency (Hz)');

ylabel('Magnitude');

title('Frequency Response of the Bandpass Filter');

grid on;

Why is there a spike from zero to -0.1?

That's the bandpass filtered sEMG signal. How can I remove the information that may be the transient response around 0 and the info around 4 seconds where I attenuate my signal in order to not include the switch (stop of an experiment)?

That's the impulse response of the bandpass filter.

### Accepted Answer

Star Strider
on 2 Oct 2024 at 14:22

Without your sampling frequeency, it is impossible to desingn the filtere and therefore difficult to figure out what the problem is. If yolu have R2018a or later, it would be eeasier to sue the bandpass function, specifically:

[signal_filtered,df] = bandpass(signal, [low_cutoff_frequency, high_cutoff_frequency], sampling_frequency, 'ImpulseReesponse','iir');

This produces a very efficient elliptiic filter and filters the signal. The second output is the digital filter object itself that you can then use with filtfilt in other places in your code, if necessary.

That aside, one way to reduce the initial transient response that filtfilt occasionally produces is:

passband = [low_cutoff_frequency, high_cutoff_frequency]/f_nyquist;

[n,Wn] = buttord(passband, passband.*[0.9 1.1], 1, 50);

[z,p,k] = butter(n, Wn, 'bandpass');

[sos,g] = zp2sos(z,p,k);

signal = [ones(50,1)*signal(1); signal];

signal_filtered = fiiltfiilt(sos, g, signal);

signal_filtered = signal_filtered(51:end);

This introduces a constant vector at the beginining of the signal that is then removed after filterinig. That should eliminate the transients. I also made other changes to your code that will produce a stable filter and a reliable filtered signal.

I did not test this code, however itt should work.

.

##### 12 Comments

Iro Liontou
on 4 Oct 2024 at 9:41

Star Strider
on 4 Oct 2024 at 10:40

My pleasure!

If I knew what you want to do, the result you want, and have your actual signal (rather than an image of it), I could probably be of more help.

The window approach works with FIR filters. You are using an IIR filter, so I doubt windowing will be effective.

If my answer helped you solve your problem, please Accept it!

Iro Liontou
on 4 Oct 2024 at 11:25

What i do is that I am trying to process a sEMG signal. What Inoticed was that when I bandpass filtered my signal and measure the FFT on the filtered signal, there was a spike on 0 Hz. That is because (what I'm thinking) is due to the 1.5 Volts my signal is. So I wanted to cut this, since maybe this is the transient response of the filtered that I used.

Thats the code:

%% sEMG Signal Analysis - Healthy Ingestion

clc;

clear;

close all;

%%%%%%%%%%%%%%%% Load sEMG Signal %%%%%%%%%%%%%%%%

% Load the saved raw sEMG signal x

load('raw_emg_signal_try_1.mat', 'x');

% Parameters

L = length(x); % Length of the signal

sampling_frequency = 1000; % Sampling frequency in Hz

duration = L / sampling_frequency;

time = (0:L-1) / sampling_frequency;

% Plot the raw signal

figure('Name', 'Raw EMG Signal');

plot(time, x);

grid on;

title('Raw EMG Signal');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, duration]);

ylim([min(x)-0.5, max(x)+0.5]); % Adjust ylim according to the raw data range

%%%%%%%%%%%%%%%% Remove DC Offset %%%%%%%%%%%%%%%%

x_dc_removed = x - mean(x); % Remove the DC component

% Plot EMG after DC Offset Removal

figure('Name','Raw EMG Signal with no DC Offset');

plot(time, x_dc_removed);

grid on;

title('EMG Signal after DC Offset Removal');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, duration]);

ylim([min(x_dc_removed)-0.5, max(x_dc_removed)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis %%%%%%%%%%%%%%%%

% Compute FFT of the raw signal with no DC offset to determine the cutoff frequencies

% Compute two sided spectrum (-Fmax:Fmax)

P1 = fft(x_dc_removed);

% Two-sided spectrum P2 for normalization of the outpower with respect to the length of the input signal

P1 = abs(P1/L);

% Compute single sided spectrum by taking the positive part of double sided spectrum and multiply by 2

P1 = P1(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

% FFT frequency range

f = sampling_frequency*(0:(L/2))/L; % Frequency vector

% Plot FFT of raw signal (Magnitude Spectrum)

figure('Name', 'Magnitude Spectrum of Raw EMG Signal');

plot(f, P1)

title('Single-Sided Amplitude Spectrum of Raw EMG Signal')

xlabel('Frequency (Hz)')

ylabel('|P1(f)|')

grid on;

%{

The plot in Figure 3 shows:

1. The majority of the signal's energy is concentrated in the lower

frequency range, below 100 Hz, which is typical for sEMG (20-150Hz).

2. No significant peak at 0 Hz means that the DC offset has been

successfully removed.

3. Small peak in the higher frequency, around 450 Hz, means there is

noise or artifacts - electrical interference from devices, electrode

movement artifacts, or high frequency muscle activity. It needs

filetring

%}

%%%%%%%%%%%%%%%% Bandpass Filter %%%%%%%%%%%%%%%%

% Bandpass filter parameters

f_nyquist = sampling_frequency / 2;

low_cutoff_frequency = 6; % Low cutoff frequency in Hz

high_cutoff_frequency = 200; % High cutoff frequency in Hz

% Design the bandpass filter

[b, a_filter] = butter(2, [low_cutoff_frequency, high_cutoff_frequency] / (f_nyquist), 'bandpass');

% Apply bandpass filter to the signal

filtered_signal_1 = filtfilt(b, a_filter, x_dc_removed);

% Plot FFT of filtered signal

figure('Name', 'Filtered EMG Signal');

plot(time, filtered_signal_1);

title('Filtered EMG Signal');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, duration]);

ylim([min(filtered_signal_1)-0.5, max(filtered_signal_1)+0.5]);

grid on;

% %%%%%%%%%%%%%%%% FFT Analysis on Filtered signal %%%%%%%%%%%%%%%%

P2_filtered_1 = fft(filtered_signal_1);

P2_filtered_1 = abs(P2_filtered_1/L); % Normalize

P2_filtered_1 = P2_filtered_1(1:L/2+1); % Single-sided spectrum

P2_filtered_1(2:end-1) = 2*P2_filtered_1(2:end-1); % Double the amplitude except for DC and Nyquist

f1 = sampling_frequency*(0:(L/2))/L;

% Plot FFT of filtered signal

figure('Name', 'Magnitude Spectrum of Filtered EMG Signal');

plot(f1, P2_filtered_1);

title('Single-Sided Amplitude Spectrum of Filtered EMG Signal')

xlabel('Frequency (Hz)')

ylabel('|P2(f)|')

grid on;

Star Strider
on 4 Oct 2024 at 12:17

That looks as though it should work. Subtracting the mean of the signal from the entire signal will eliminate tthe D-C (constant) offset, as will highpass-filttering it with the passband set at about 1 Hz (and with a sampling frequency of 1 kHz, that should be enough), or bandpass-filtering with the lower stopband set to about 1 Hz..

A frequency-selective filter will not eliminate broadband noise, so if you want to minimise the noise, use the sgolayfilt function. I generally choose a 3-order polynomial and then set ‘framelen’ to give the result I want (least noise with best signal fidelity, so that I do not lose any significant signal features). I generally start with ‘framelen’ set to about 1% of the signal length, and go from there. (The sgolayfilt function essentially acts as a FIR comb filter, as can be seen by dividing the FFT of the filtered signal by the FFT of the original signal.)

My guess is that you are acquiring the EMG yourself. Instrumentation that uses shielded coaxial cables for tthe leads, and a neutral reference electrode (erroneously sometimes called a ‘ground’ lead) can eliminate a significant amount of ambient signal noise, and 50-60 Hz mains frequency noise. (The reference signal is subtracted in analog hardware from the ‘active’ signal leads, prior to the ADC stage. Its use can also eliminate the D-C offset.)

I will do what I can to help you with this.

It would probably help to have your signal.

I wrote a functon for the one-sided Fourier transform:

function [FTs1,Fv] = FFT1(s,t)

% Arguments:

% s: Signal Vector Or Matrix

% t: Associated Time Vector

t = t(:);

L = numel(t);

if size(s,2) == L

s = s.';

end

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));

Fv = Fs*(0:(NFFT/2))/NFFT;

% Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

Fv = Fv(:);

FTs1 = FTs(Iv,:);

end

There is nothing magic about it, other than saving me the time and trouble of typing all that in whenever II want to calculate the Fourier transform of a signal. You are free to use it as you wish. It automatically eliminates any D-C offseet, and windows the signal to increase the reliability of the result.

.

Iro Liontou
on 6 Oct 2024 at 14:48

Sorry for the delay. I was experimenting. Well my signal is that. Is the x variable.

%% sEMG Signal Analysis - Healthy Ingestion

clc;

clear;

close all;

%%%%%%%%%%%%%%%% Load sEMG Signal %%%%%%%%%%%%%%%%

% Load the saved raw sEMG signal x

load('raw_emg_signal_try_1.mat', 'x');

% Parameters

L = length(x); % Length of the raw sEMG signal

sampling_frequency = 1000; % Sampling frequency in Hz

duration = L / sampling_frequency;

time = (0:L-1) / sampling_frequency;

% Plot the raw signal

figure('Name', 'Raw sEMG Signal - Healthy Ingestion');

plot(time, x);

grid on;

title('Raw sEMG Signal - Healthy Ingestion');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, duration]);

ylim([min(x)-0.5, max(x)+0.5]); % Adjust ylim according to the raw data range

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Pre-Processing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% Truncate raw sEMG signal %%%%%%%%%%%%%%%%

% Truncate signal to 4 seconds

max_time = 4.3;

num_samples = max_time * sampling_frequency; % Number of samples for 4.3 seconds

% x is truncated to 4.3 seconds of data

x_truncated = x(1:num_samples);

time_truncated = (0:num_samples-1) / sampling_frequency;

% Plot the truncated signal

figure('Name', 'sEMG Signal (Truncated to 4 seconds)');

plot(time_truncated, x_truncated);

grid on;

title('sEMG Signal Truncated to 4 Seconds');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([1.5 - 0.5, 1.5 + 0.5]);

%%%%%%%%%%%%%%%% Remove DC Offset %%%%%%%%%%%%%%%%

% Truncated signal with no DC offset

x_truncated_dc_removed = x_truncated - mean(x_truncated);

% Plot sEMG after DC Offset Removal

figure('Name', 'sEMG Signal with no DC Offset (Truncated to 4 seconds)');

plot(time_truncated, x_truncated_dc_removed);

grid on;

title('sEMG Signal with no DC Offset Truncated to 4 Seconds');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_truncated_dc_removed)-0.5, max(x_truncated_dc_removed)+0.5]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 1st Filter - Apply Savitzky-Golay Filter on raw sEMG signal %%%%%%%%%%%%%%%%

% Parameters

frame_len = round(0.01 * length(x_truncated_dc_removed)); % Set framelen to 1% of signal length

if mod(frame_len, 2) == 0 % Ensure framelen is odd

frame_len = frame_len + 1;

end

% Apply Savitzky-Golay filter for noise reduction

x_sgolay_filtered = sgolayfilt(x_truncated_dc_removed, 4, frame_len);

% Plot the smoothed signal after sgolayfilt

figure('Name', 'Savitzky-Golay Filtered sEMG signal');

plot(time_truncated, x_sgolay_filtered);

grid on;

title('Savitzky-Golay Filtered sEMG signal');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_sgolay_filtered)-0.5, max(x_sgolay_filtered)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis after Savitzky Filter %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_sgolay = length(x_sgolay_filtered);

% Compute two-sided spectrum (-Fmax:Fmax)

Psg = fft(x_sgolay_filtered);

% Two-sided spectrum P2 for normalization of the outpower with respect to the length of the input signal

Psg = abs(Psg/L_sgolay);

% Compute single-sided spectrum by taking the positive part of double-sided spectrum and multiply by 2

Psg = Psg(1:L_sgolay/2+1);

Psg(2:end-1) = 2*Psg(2:end-1);

% FFT frequency range

fsg = sampling_frequency*(0:(L_sgolay/2))/L_sgolay; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Savitzky-Golay Filtering');

plot(fsg, Psg);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Savitzky-Golay Filtering')

xlabel('Frequency (Hz)');

ylabel('|Psg(f)|');

grid on;

%%%%%%%%%%%%%%%% Create and Apply the Hann Window %%%%%%%%%%%%%%%%

% Create a Hann window

hann_window = hann(num_samples); % Create a Hann window of the same length as the truncated signal

% Apply Hann window to the sgolayfilt signal

x_windowed_sgolay_bd = x_sgolay_filtered .* hann_window;

% Plot the effect of the Hann Window

figure('Name', 'Effect of Hann Window on sEMG Signal');

plot(time_truncated, x_windowed_sgolay_bd);

grid on;

title('sEMG Signal (After Hann Window Applied)');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_windowed_sgolay_bd)-0.5, max(x_windowed_sgolay_bd)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis after Savitzky Filter and Hanning %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_w_sgolay = length(x_windowed_sgolay_bd);

% Compute two-sided spectrum (-Fmax:Fmax)

P_w_sg = fft(x_windowed_sgolay_bd);

% Two-sided spectrum P2 for normalization of the outpower with respect to the length of the input signal

P_w_sg = abs(P_w_sg/L_w_sgolay);

% Compute single-sided spectrum by taking the positive part of double-sided spectrum and multiply by 2

P_w_sg = P_w_sg(1:L_w_sgolay/2+1);

P_w_sg(2:end-1) = 2*P_w_sg(2:end-1);

% FFT frequency range

f_w_sg = sampling_frequency*(0:(L_w_sgolay/2))/L_w_sgolay; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Savitzky-Golay Filtering and Hanning');

plot(f_w_sg, P_w_sg);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Savitzky-Golay Filtering and Hanning')

xlabel('Frequency (Hz)');

ylabel('|P_w_sg(f)|');

grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 2nd Filter - Bandpass Filter on raw sEMG signal %%%%%%%%%%%%%%%%

% Design the bandpass filter

f_nyquist = sampling_frequency / 2;

low_cutoff_frequency = 6; % Low cutoff frequency in Hz

high_cutoff_frequency = 100; % High cutoff frequency in Hz

[b, a_filter] = butter(4, [low_cutoff_frequency, high_cutoff_frequency] / (f_nyquist), 'bandpass');

% Apply bandpass filter to the windowed signal

bandpass_filtered_signal = filtfilt(b, a_filter, x_truncated_dc_removed);

% Plot the filtered signal after hanning window

figure('Name', 'Bandpass Filtered sEMG Signal');

plot(time_truncated, bandpass_filtered_signal);

title('Bandpass Filtered sEMG Signal ');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(bandpass_filtered_signal)-0.5, max(bandpass_filtered_signal)+0.5]);

grid on;

%%%%%%%%%%%%%%%% FFT Analysis after Bandpass Filter %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_bandpass = length(bandpass_filtered_signal);

% Compute two-sided spectrum (-Fmax:Fmax)

Pbd = fft(bandpass_filtered_signal);

Pbd = abs(Pbd/L_bandpass);

% Compute single-sided spectrum by taking the positive part of double-sided spectrum and multiply by 2

Pbd = Pbd(1:L_bandpass/2+1);

Pbd(2:end-1) = 2*Pbd(2:end-1);

% FFT frequency range

fbd = sampling_frequency*(0:(L_bandpass/2))/L_bandpass; % Frequency vector

% Call FFT1 on the processed signal (after bandpass filtering)

% [FTT_bandpass, freq_bandpass] = FFT1(bandpass_filtered_signal, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Bandpass Filter');

plot(fbd, Pbd);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Bandpass Filter')

xlabel('Frequency (Hz)');

ylabel('|Pbd(f)|');

grid on;

%%%%%%%%%%%%%%%% Create and Apply the Hann Window %%%%%%%%%%%%%%%%

% Create a Hann window

hann_window = hann(num_samples); % Create a Hann window of the same length as the truncated signal

% Apply Hann window to the sgolayfilt signal

x_windowed_bandpass = bandpass_filtered_signal .* hann_window;

% Plot the effect of the Hann Window

figure('Name', 'Effect of Hann Window on sEMG Signal');

plot(time_truncated, x_windowed_bandpass);

grid on;

title('sEMG Signal (After Hann Window Applied)');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_windowed_bandpass)-0.5, max(x_windowed_bandpass)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis after Bandpass Filter and Hanning %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_w_bandpass = length(x_windowed_bandpass);

% Compute two-sided spectrum (-Fmax:Fmax)

P_w_sg = fft(x_windowed_bandpass);

P_w_sg = abs(P_w_sg/L_w_sgolay);

P_w_sg = P_w_sg(1:L_w_sgolay/2+1);

P_w_sg(2:end-1) = 2*P_w_sg(2:end-1);

% FFT frequency range

f_w_sg = sampling_frequency*(0:(L_w_sgolay/2))/L_w_sgolay; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Bandpass Filtering and Hanning');

plot(f_w_sg, P_w_sg);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Bandpass Filtering and Hanning')

xlabel('Frequency (Hz)');

ylabel('|P_w_sg(f)|');

grid on;

Well maybe that's a lot of code, but i was experimenting to see the differences. What I did is to apply the Savitzky-Golay Filter as you said and then the Hanning window. I did the same for the bandpass filter which I use the IIR filter (butterworth). I hope this filter does not create other problems to the processing of the signal. Last but not least, i tried to apply first the Savitzly-Golay filter then the Bandpass filter the that and then the Hanning Window.

What i noticed is that probably for the bandpass filter the low cutoff frequency should be 1 Hz and not 6Hz and also the Hanning window that I created tapers some data of my signal like the one peak that there is in the 4th second. How can I fix taht? And is it a good approach to sgolayfilt first, then bandpass and apply teh hanning window? Thank you!

Star Strider
on 6 Oct 2024 at 15:28

Your posted code would appear to work. I cannot determine that however, or reply to specific questions about how your code works without having‘raw_emg_signal_try_1.mat’ to run with it.

I generally do any frequency-selective filtering first, and use sgolayfilt after that, although I am not certain that the order is truly important..

Iro Liontou
on 6 Oct 2024 at 15:51

I uploaded the 'raw_emg_signal_try_1.mat', but it;s ok Thank you very much for your help!!!

Star Strider
on 6 Oct 2024 at 17:08

As always, my pleasure!

I sttill do not see the .mat file, only the .m file. I looked through the .m file to see if there was any attached or included data, however it was the same as you posted in your earlier Comment, without the data. If the .mat file is too large, use the zip function to zip it to a .zip file to see if it is then small enough to upload. (I believe the size limit is 5 MB, regardless of the file type.)

Iro Liontou
on 6 Oct 2024 at 18:26

Star Strider
on 6 Oct 2024 at 20:21

I made a feww changes, specifically:

xc = find(ischange(x, 'Threshold',1))-1; % Find Index Just Before Step Change

txc = time(xc)

time = time(1:xc); % Truncate Vector

x = x(1:xc); % Truncate Vector

mean_x = mean(x)

x = x - mean_x; % Subtract ‘mean’

to elilminate the ‘step’ discontinuity at the end. This made the rest of the signal processiing a bit easier. It also eliminated the filtter transients. (I saved ‘mean_x’ separately, in the event you want to add it back later.)

I also added a Fourier transform of the original signal (after these changes), and note that a significant amount of the signal energy is above 400 Hz, and half of it is above 250 Hz. Your bandpass filter has an upper cutoff of 100 Hz.

This is not a criticism of your analysis, simiply an observation. I am not certain what sort of analysis you are doing.

Your code with my changes —

%% sEMG Signal Analysis - Healthy Ingestion

clc;

clear;

close all;

%%%%%%%%%%%%%%%% Load sEMG Signal %%%%%%%%%%%%%%%%

% Load the saved raw sEMG signal x

load('raw_emg_signal_try_1.mat', 'x');

% Parameters

L = length(x); % Length of the raw sEMG signal

sampling_frequency = 1000; % Sampling frequency in Hz

nyquist_frequency = sampling_frequency/2;

duration = L / sampling_frequency;

time = (0:L-1) / sampling_frequency;

xc = find(ischange(x, 'Threshold',1))-1; % Find Index Just Before Step Change

txc = time(xc)

txc = 4.5730

time = time(1:xc); % Truncate Vector

x = x(1:xc); % Truncate Vector

mean_x = mean(x)

mean_x = 1.4814

x = x - mean_x; % Subtract ‘mean’

% Plot the raw signal

figure('Name', 'Raw sEMG Signal - Healthy Ingestion');

plot(time, x)

% hold on

% plot(time(xc), x(xc), 'r+')

% hold off

grid on;

title('Raw sEMG Signal - Healthy Ingestion');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, duration]);

ylim([min(x)-0.5, max(x)+0.5]); % Adjust ylim according to the raw data range

%%%%%%%%%%%%%%%% FFT Analysis — Original Signal %%%%%%%%%%%%%%%%

% Compute FFT of the original signal to determine the cutoff frequencies

L = length(x);

NFFT = 2^nextpow2(L); % For Efficiency

% Compute two-sided spectrum (-Fmax:Fmax)

Psg = fft((x(:) - mean(x)).*hann(L), NFFT); % Subtract ‘mean’, Add ‘hann’ Window

Psg = abs(Psg/L);

% % Psg = Psg(1:L/2+1);

Psg = Psg(1:NFFT/2+1);

Psg(2:end-1) = 2*Psg(2:end-1);

% FFT frequency range

% % fsg = sampling_frequency*(0:(L/2))/L; % Frequency vector

fsg = sampling_frequency*(0:(NFFT/2))/NFFT; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of Original sEMG Signal');

plot(fsg, Psg);

title('Single-Sided Amplitude Spectrum of Original sEMG Signal')

xlabel('Frequency (Hz)');

ylabel('|Psg(f)|');

grid on;

ylim([0 0.006])

Psg_Int = cumtrapz(fsg, Psg);

Psg_Int_mdn = median(Psg_Int)

Psg_Int_mdn = 0.0688

freq_Psg_Int_mdn = interp1(Psg_Int, fsg, Psg_Int_mdn);

figure

plot(fsg, Psg_Int)

grid

xline(freq_Psg_Int_mdn,'--',"Associated Frequency ("+freq_Psg_Int_mdn+" Hz)")

yline(Psg_Int_mdn, '--', 'Median Signal Energy')

xlabel('Frequency (Hz)')

ylabel('Cumulative Signal Energy')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Pre-Processing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% Truncate raw sEMG signal %%%%%%%%%%%%%%%%

% Truncate signal to 4 seconds

max_time = 4.3;

num_samples = max_time * sampling_frequency; % Number of samples for 4.3 seconds

% x is truncated to 4.3 seconds of data

x_truncated = x(1:num_samples);

time_truncated = (0:num_samples-1) / sampling_frequency;

% Plot the truncated signal

figure('Name', 'sEMG Signal (Truncated to 4 seconds)');

plot(time_truncated, x_truncated);

grid on;

title('sEMG Signal Truncated to 4 Seconds');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

% ylim([1.5 - 0.5, 1.5 + 0.5]);

ylim('padded')

%%%%%%%%%%%%%%%% Remove DC Offset %%%%%%%%%%%%%%%%

% Truncated signal with no DC offset

x_truncated_dc_removed = x_truncated - mean(x_truncated);

% Plot sEMG after DC Offset Removal

figure('Name', 'sEMG Signal with no DC Offset (Truncated to 4 seconds)');

plot(time_truncated, x_truncated_dc_removed);

grid on;

title('sEMG Signal with no DC Offset Truncated to 4 Seconds');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_truncated_dc_removed)-0.5, max(x_truncated_dc_removed)+0.5]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 1st Filter - Apply Savitzky-Golay Filter on raw sEMG signal %%%%%%%%%%%%%%%%

% Parameters

frame_len = round(0.01 * length(x_truncated_dc_removed)); % Set framelen to 1% of signal length

if mod(frame_len, 2) == 0 % Ensure framelen is odd

frame_len = frame_len + 1;

end

% Apply Savitzky-Golay filter for noise reduction

x_sgolay_filtered = sgolayfilt(x_truncated_dc_removed, 4, frame_len);

% Plot the smoothed signal after sgolayfilt

figure('Name', 'Savitzky-Golay Filtered sEMG signal');

plot(time_truncated, x_sgolay_filtered);

grid on;

title('Savitzky-Golay Filtered sEMG signal');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_sgolay_filtered)-0.5, max(x_sgolay_filtered)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis after Savitzky Filter %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_sgolay = length(x_sgolay_filtered);

% Compute two-sided spectrum (-Fmax:Fmax)

Psg = fft(x_sgolay_filtered);

Psg = abs(Psg/L_sgolay);

Psg = Psg(1:L_sgolay/2+1);

Psg(2:end-1) = 2*Psg(2:end-1);

% FFT frequency range

fsg = sampling_frequency*(0:(L_sgolay/2))/L_sgolay; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Savitzky-Golay Filtering');

plot(fsg, Psg);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Savitzky-Golay Filtering')

xlabel('Frequency (Hz)');

ylabel('|Psg(f)|');

grid on;

%%%%%%%%%%%%%%%% Create and Apply the Hann Window %%%%%%%%%%%%%%%%

% Create a Hann window

hann_window = hann(num_samples); % Create a Hann window of the same length as the truncated signal

% Apply Hann window to the sgolayfilt signal

x_windowed_sgolay_bd = x_sgolay_filtered .* hann_window;

% Plot the effect of the Hann Window

figure('Name', 'Effect of Hann Window on sEMG Signal');

plot(time_truncated, x_windowed_sgolay_bd);

grid on;

title('sEMG Signal (After Hann Window Applied)');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_windowed_sgolay_bd)-0.5, max(x_windowed_sgolay_bd)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis after Savitzky Filter and Hanning %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_w_sgolay = length(x_windowed_sgolay_bd);

% Compute two-sided spectrum (-Fmax:Fmax)

P_w_sg = fft(x_windowed_sgolay_bd);

P_w_sg = abs(P_w_sg/L_w_sgolay);

P_w_sg = P_w_sg(1:L_w_sgolay/2+1);

P_w_sg(2:end-1) = 2*P_w_sg(2:end-1);

% FFT frequency range

f_w_sg = sampling_frequency*(0:(L_w_sgolay/2))/L_w_sgolay; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Savitzky-Golay Filtering and Hanning');

plot(f_w_sg, P_w_sg);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Savitzky-Golay Filtering and Hanning')

xlabel('Frequency (Hz)');

ylabel('|P_w_sg(f)|');

grid on;

%%%%%%%%%%%%%%%% 2nd Filter - Bandpass Filter on raw sEMG signal %%%%%%%%%%%%%%%%

% Design the bandpass filter

f_nyquist = sampling_frequency / 2;

low_cutoff_frequency = 2; % Low cutoff frequency in Hz

high_cutoff_frequency = 100; % High cutoff frequency in Hz

bandwidth = [low_cutoff_frequency, high_cutoff_frequency] / (f_nyquist);

[b, a_filter] = butter(4, [low_cutoff_frequency, high_cutoff_frequency] / (f_nyquist), 'bandpass');

% Apply bandpass filter to the windowed signal

% x_truncated_dc_removed = [ones(100,1)*x_truncated_dc_removed(1); x_truncated_dc_removed; ones(100,1)*x_truncated_dc_removed(end)];

bandpass_filtered_signal = filtfilt(b, a_filter, x_truncated_dc_removed);

% bandpass_filtered_signal = bandpass_filtered_signal(101:end-100);

% Plot the filtered signal after hanning window

figure('Name', 'Bandpass Filtered sEMG Signal');

plot(time_truncated, bandpass_filtered_signal);

title('Bandpass Filtered sEMG Signal ');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(bandpass_filtered_signal)-0.5, max(bandpass_filtered_signal)+0.5]);

grid on;

%%%%%%%%%%%%%%%% FFT Analysis after Bandpass Filter %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_bandpass = length(bandpass_filtered_signal);

% Compute two-sided spectrum (-Fmax:Fmax)

Pbd = fft(bandpass_filtered_signal);

Pbd = abs(Pbd/L_bandpass);

Pbd = Pbd(1:L_bandpass/2+1);

Pbd(2:end-1) = 2*Pbd(2:end-1);

% FFT frequency range

fbd = sampling_frequency*(0:(L_bandpass/2))/L_bandpass; % Frequency vector

% Call FFT1 on the processed signal (after bandpass filtering)

% [FTT_bandpass, freq_bandpass] = FFT1(bandpass_filtered_signal, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Bandpass Filter');

plot(fbd, Pbd);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Bandpass Filter')

xlabel('Frequency (Hz)');

ylabel('|Pbd(f)|');

grid on;

%%%%%%%%%%%%%%%% Create and Apply the Hann Window %%%%%%%%%%%%%%%%

% Create a Hann window

hann_window = hann(num_samples); % Create a Hann window of the same length as the truncated signal

% Apply Hann window to the sgolayfilt signal

x_windowed_bandpass = bandpass_filtered_signal .* hann_window;

% Plot the effect of the Hann Window

figure('Name', 'Effect of Hann Window on sEMG Signal');

plot(time_truncated, x_windowed_bandpass);

grid on;

title('sEMG Signal (After Hann Window Applied)');

xlabel('Time (s)');

ylabel('Voltage (V)');

xlim([0, max_time]);

ylim([min(x_windowed_bandpass)-0.5, max(x_windowed_bandpass)+0.5]);

%%%%%%%%%%%%%%%% FFT Analysis after Bandpass Filter and Hanning %%%%%%%%%%%%%%%%

% Compute FFT of the smoothed signal to determine the cutoff frequencies

L_w_bandpass = length(x_windowed_bandpass);

% Compute two-sided spectrum (-Fmax:Fmax)

P_w_sg = fft(x_windowed_bandpass);

P_w_sg = abs(P_w_sg/L_w_sgolay);

P_w_sg = P_w_sg(1:L_w_sgolay/2+1);

P_w_sg(2:end-1) = 2*P_w_sg(2:end-1);

% FFT frequency range

f_w_sg = sampling_frequency*(0:(L_w_sgolay/2))/L_w_sgolay; % Frequency vector

% Call FFT1 on the processed signal

% [FTT_sgolay, freq_sgolay] = FFT1(x_sgolay_filtered, time_truncated);

% Plot FFT of the smoothed signal with no DC Offset after Savitzky Filtering

figure('Name', 'Magnitude Spectrum of sEMG Signal after Bandpass Filtering and Hanning');

plot(f_w_sg, P_w_sg);

title('Single-Sided Amplitude Spectrum of sEMG Signal after Bandpass Filtering and Hanning')

xlabel('Frequency (Hz)');

ylabel('|P_w_sg(f)|');

grid on;

I have no further suggestions.

.

Iro Liontou
on 7 Oct 2024 at 5:05

Star Strider
on 7 Oct 2024 at 10:13

As always, my pleasure!

I am not certain what your experimental setup is, however consider adding a reference electrode, if you are not already using one. Subtracting that signal from your signal-of-interest might eliminate the high-frequency noise, if it is (as I suspect) external interference.

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)