How to Correct the baseline shift of the signal

I am trying to correct the baseline of the signal (picture attached) by using a single exponential curve fit to adjust the baseline.

The data file is attached as well.

Answer 1

[Star Strider]

I always use discrete filters to eliminate baseline drift, baseline offset, and high-frequency noise.

Try this:

[D,S] = xlsread('Book1GA.xls');
Tcl = regexp(S, '(\d\.\d+)', 'match');                              % Time Values (Cell Array)
Tch = cell2mat([Tcl{:}]');                                          % Time Values (Character Array)
Tv = str2num(Tch);                                                  % Time Vector (Double)
L = numel(Tv);
Ts = mean(diff(Tv));                                                % Sampling Interval
Tsd = std(diff(Tv));
Fs = 1/Ts;                                                          % Sampling Frequency
Fn = Fs/2;                                                          % Nyquist Frequency
Dmd = D - mean(D);
FT_D = fft(Dmd)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                                 % Frequency Vector
Iv = 1:numel(Fv);

figure(1)
plot(Fv, abs(FT_D(Iv))*2)
grid
axis([0  25    ylim])
set(gca, 'XMinorTick','on')
xlabel('Frequency')

Wp = [2.1  15.0]/Fn;                                                % Passband Frequency (Normalised)
Ws = [1.8  18.0]/Fn;                                                % Stopband Frequency (Normalised)
Rp =   1;                                                           % Passband Ripple (dB)
Rs =  50;                                                           % Stopband Ripple (dB)
[n,Ws]  = cheb2ord(Wp,Ws,Rp,Rs);                                    % Filter Order
[z,p,k] = cheby2(n,Rs,Ws);                                          % Filter Design, Sepcify Bandpass
[sos,g] = zp2sos(z,p,k);                                            % Convert To Second-Order-Section For Stability

figure(2)
freqz(sos, 2^16, Fs)                                                % Filter Bode Plot

D_Filtered = filtfilt(sos, g, D);                                   % Filter Signal

figure(3)
plot(Tv, D_Filtered)
grid
xlabel('Time')

Experiment witht he passband and stopband frequencies to get the result you want. Note that in a bandpass filter, the passband frequency limnits must always be within the stopband frequency limits.

EDIT — Adding the plot:

[amirhossein shahriari]

thank you for your great answer, How can I measure the beat rate of the signal like how many action potentials are occurring per minute

[Star Strider]

As always, my pleasure.

Probably the easiest way to measure the rate is to use the findpeaks (link) function. If you return the indices in the ‘locs’ output, you can calculate the rate as:

[pks, locs] = findpeaks(D_Filtered, 'MinPeakHeight',100, 'MinPeakDist',20);     % Determine Peaks & Indices
T_Pks = Tv(locs);                                                               % Peak Times
Rate = 1/mean(diff(T_Pks));                                                     % Calculate Rate (Peaks/Time Unit)

To then calculate the number of action potentials are occurring per minute, multiply the rate per time unit by the time units per minute, so if the time units are seconds:

RatePerMinute = Rate * 60;                                                      % Rate/Min = Rate/Sec * Sec/Min

Here, if the time units are seconds, the rate per minute is about 140.

[amirhossein shahriari]

I don't quite understand some of the steps especially signal processing ones. Can you please explain for example why would I have to find frequencies and how did you get those numbers for passband and stopband and a general description of the code? Thank You

[Star Strider]

You would use the fft results to determine the frequency content of your signal. It is then necessary to decide what frequencies are important in your signal, and then design the filter to do what you want. You would need to experiment with the frequencies to see the results of different choices.

I chose the passband frequencies to eliminate the low-frequency baseline drift, and some of the high-frequency noise, based on the plot of the fft of your original signal. (The lower stopband frequency must be less than the lower passband frequency, and the higher stopband frequency must be higher than the higher passband frequency in a bandpass filter.) I believe these frequencies do you want. If they do not, you can certainly change them.

The passband ripple and stopband ripple (attenuation) can also be whatever you want. I chose 1 dB and 50 dB respectively because that is usually adequate. You can certainly change them as well.

I chose a Chebyshev Type II filter because it has a flat passband, and tolerates short transition regions. The cheb2ord function calculates the optimal order for the filter. The cheby2 function calculates the zero-pole-gain filter coefficients here. (It can produce other outputs as well.) The zp2sos function converts the coefficients to second order section representation, since that produces as stable and reliable filter.

In figure(2), the freqz function calculates and displays the filter frequency response (Bode plot).

The filtfilt function creates a linear-phase version of the filtered signal, preventing phase distortion.

In figure(3), the code plots the filtered signal in the time domain.

If you want to eliminate the baseline drift but not the baseline offset, calculate the mean of your original data, then add it back to the filtered data, ‘D_Filtered’ in my code.

[amirhossein shahriari]

Thank you very much Sorry for bothering again, I didn’t get how you got those numbers for Wp and for Ws [2.1 15.0] [1.8 18.0] I couldn’t find them on the plot

[Star Strider]

My pleasure.

They’re not on the plot, because at the time there was no reason to plot the passband limits I chose. This code for figure(1) includes them:

figure(1)
plot(Fv, abs(FT_D(Iv))*2)
grid
axis([0  25    ylim])
set(gca, 'XMinorTick','on')
xlabel('Frequency')
hold on
plot(2.1*[1 1], ylim, '--r',    15*[1 1], ylim, '--r')
hold off

The rest of my code is unchanged.

[amirhossein shahriari]

thank you , I have one more question how did you choose these passband limits and stopband limits like [2.1 15] and [1.8 18]?

[Star Strider]

My pleasure.

I chose the passband frequencies from the fft plot, and then tweaked them to get what I considered to be an acceptable time-domain output. I chose the stopband frequencies so the ratios between the respective frequencies were about 1.2:

Ratio = [2.1/1.8  18/15];

So I could have calculated them as:

Ws = Wp.*[1/1.2  1.2];

I rounded up the lower stopband from 1.75 to 1.8. (I usually use this sort of ‘shorthand’ calculation in my own code.) The goal is to create stable filters that produce the results I want.

There is nothing magic about any of this, and different designs are certainly possible, depending on what you want as a result. Other filters, such as elliptical filters, are efficient and work well.

If you have R2018a, experiment with the bandpass (link) and related filter functions. (I believe that function uses an elliptic filter.) Note that if you choose to use the filter designed by the second ‘d’ output, the documentation recommends using the filter function, that unfortunately produces phase distortion in the filter output. I recommend using the filtfilt function because it is a linear-phase filter, and because the bandpass function itself uses filtfilt internally to calculate the first output.

[amirhossein shahriari]

Thank you so much . I have tried to experiment other sets of my data with savitzky filter but still seems like this band pass filter you design are way better and I can get way better average of my action potentials. So do you know why some other filters like savitzky are not helpful enough for this type of data

[Star Strider]

As always, my pleasure.

The Savitzky-Golay filter is primarily used for curve-smoothing (or at least that is my experience). It is quite good for that purpose. However for signal processing, I find it difficult to design functional filters with it, and prefer using frequency-selective filters such as those I use here. The Tips (link) section of the documentation describes the use of the Savitzky-Golay filter much better than I can.

Also, if you have R2018a or later, the bandpass (link) and related functions make filter design significantly easier. If you use the ‘d’ output from it with an IIR filter, remember to use the filtfilt function to do the actual filtering, not the filter function as suggested in the documentation (that I consider to be an error).

Answer 2

[Image Analyst]

You might use fitnlm(), if you have the Statistics and Machine Learning Toolbox

[D, S] = xlsread('Book1GA.xls');
Tcl = regexp(S, '(\d\.\d+)', 'match');                              % Time Values (Cell Array)
Tch = cell2mat([Tcl{:}]');                                          % Time Values (Character Array)
Tv = str2num(Tch);   % Thanks to Star for figuring out how to get the time!
plot(Tv, D, 'b-')
grid on;
erodedSignal = imerode(D, ones(51, 1));
hold on;
plot(Tv, erodedSignal, 'r-', 'LineWidth', 2);

X = Tv;
Y = erodedSignal;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X, Y);
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) + b(2) * exp(-b(3)*x(:, 1));  
beta0 = [10000, 300, 1]; % Guess values to start with.  Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!

% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) + coefficients(2) * exp(-coefficients(3)*X);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'k-', 'LineWidth', 3);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Original Signal', 'Baseline', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

The black line is the exponential fit of the bottoms of the signal.

[amirhossein shahriari]

thank you so much for your reply , how can I subtract the exponential fit of the data

[Image Analyst]

Simply subtract yFitted from D

yCorrected = D - yFitted;

[joshua marvin]

Sorry this is on an old post but is there a way to get the model to be e^-x^2 it tells me to use .^2 but then it givs me not a number or inf errors

[Image Analyst]

@joshua marvin have the model be this:

modelfun = @(b,x) b(1) + b(2) * exp(-b(3)*x(:, 1) .^ 2);  

Answer 3

[Image Analyst]

Well you can try imerode() or movmin(). Or you can try some more sophisticated spectroscopic methods. It would have been easy for people to try things if you had attached your data. Try this link