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

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.

Answer 1

[Star Strider]

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.

.

[Iro Liontou]

It does work. It hasn't that much of difference. Maybe i should try a windowing method, like Hann Window. Thank you very much for the help @Star Strider!

[Star Strider]

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]

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]

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]

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);

% Two-sided spectrum P2 for normalization of the outpower with respect to the length of the input 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);

% 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 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]

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]

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

[Star Strider]

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]

I think now should be fine. If you would like to recommend anything different, it would be a huge help. Thank you!

[Star Strider]

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
% Two-sided spectrum P2 for normalization of the outpower with respect to the length of the input signal
Psg = abs(Psg/L);
% Compute single-sided spectrum by taking the positive part of double-sided spectrum and multiply by 2
% % 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); 
% 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 = 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); 
% Two-sided spectrum P2 for normalization of the outpower with respect to the length of the input 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); 
% 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 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]

The analysis that I do is from a sEMG signal of a healthy ingestion. I'M still working on what the 450 Hz might be. I'm thinking taht it is some kind of noise or artifact since it is high frequency and also because the frequency range of a muscle contraction is around 50-150 Hz. Thank you once more very much for your time and help!!!

[Star Strider]

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.