Help With Starting Point for a Low Pass Filter

I am brand new to signal processing and having a lot of trouble understanding the documentation ot implement a low pass filter, which is what I think I need. An exmaple of what I'm working with is shown in Fig1.

The signal is oxygen amperometry changes which are necessarily slow. I've been attempting to use matlab find peaks function which does a decent job at finding peaks, but I cannot get it to find an appropriate width as it is considering all the smaller peaks instead of the one large one. For example, see Fig2. I'm primarily interested in peaks 1,2, and 4. The findpeaks does a good job with the widths of peaks 1 and 3, but in 4 it's considering that a small peak. Changing factors within the findpeaks function also doesn't appear to help as even though the smaller peaks may not appear in the output, findpeaks still considers them in the anlysis. In Fig2 I've been able to have findpeaks suppress the output of peaks near peak 4, but that doesn't help its width estimation.

I'm pretty sure the answer to my issue is preprocessing, but a smoothing such as moving average doesn't quite capture what I need and so I think I need a more advanced filter like a low pass. Can anyone suggest starting parameters in some code for me to work with ? I'm having trouble getting any implementation to work. An example dataset is in ExData.xlsx where time point 0 is the event of interest.

Answer 1

[Star Strider]

I am not certain what result you want.

Here is a prototype lowpass filter design you can experiment with. It has a passband frequency of 0.08 Hz (it stops everything above that) and a reasonably short transition region. I will help you refine it to do what you want, and help with findpeaks as well. (With respect to findpeaks, you may find the code in this Comment (link) relevant with respect to selecting your peaks of interest.)

The Filter —

[D,S] = xlsread('ExData.xlsx');
t = D(:,1);
s = D(:,2);

figure
plot(t,s)                                               % Plot Original Signal
grid

Ts = mean(diff(t));                                     % Sampling Interval
L = numel(t);
Fs = 1/Ts;                                              % Sampling Frequency
Fn = Fs/2;                                              % Nyquist Frequency
smc = s-mean(s);                                        % Subtract Mean (0 Hz)
FTs = fft(smc)/L;                                       % Fourier Transform
Fv = linspace(0,1,fix(L/2)+1)*Fn;                       % Frequency Vector
Iv = 1:numel(Fv);                                       % Index Vector

figure
plot(Fv, abs(FTs(Iv))*2)                                % Plot Fourier Transform
grid
xlim([0 1])

Ws = 0.09/Fn;                                           % Lowpass Stopband Frequency
Wp = 0.08/Fn;                                           % Lowpass Passband Frequency
Rp =  1;                                                % Passband Ripple
Rs = 50;                                                % Stopband Ripple & Attenuation
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                         % Elliptical Filter Parameters
[z,p,k] = ellip(n,Rp,Rs,Wp,'low');                      % Elliptical Filter Design (Specify Lowpass)
[sos,g] = zp2sos(z,p,k);                                % Second-Order-Section For Stability

figure
freqz(sos, 2^14, Fs)                                    % Plot Filter Characteristics
set(subplot(2,1,1),'XLim',[0 2.5])
set(subplot(2,1,2),'XLim',[0 2.5])

s_filt = filtfilt(sos,g,s);                             % Filter Signal

figure
plot(t, s_filt)                                         % Plot Filtered Signal
grid

Experiment with the filter. Remember that with a lowpass filter, the stopband frequency (Ws) has to be greater than the passband frequency (Wp). The Fourier transform plot will help you decide on the frequencies.

[Aaron Smith]

Thanks a lot for your help. I'm about half way through this, but can you answer some questions?

Why are we mean centering?

What is the logic behind variables FV and IV?

The plot of just that looks great, and it is really similar to what I had been making but my Y valus for the fft were above 1 and your code isn't. Can you help me understand that? Sorry for basic questions but am very new to this.

[Star Strider]

My pleasure.

I subtracted the mean of the signal from the signal because otherwise the very high amplitude of the d-c (constant) offset obscures the other peaks in the fft result. The 0 Hz frequency is the mean of the signal.

The ‘Fv’ vector is the frequency vector for the one-sided Fourier transform. It creates a vector from 0 to 1 and is half the length of the fft result, then multiplies it by the Nyquist frequency, since that is the highest uniquely-resolvable frequency in a sampled signal. The index vector ‘Iv’ simply makes the plotting and other references to the fft result easier.

It is necessary to normalise the fft output by the length of the signal, here dividing it by ‘L’, so that will decrease the amplitude of the peaks to about what they should be. Because the fft function returns a two-sided Fourier transform, the energy of the waveforms in the time-domain signal are divided between the two halves. So it is then necessary to multiply by 2 to recover the approximate amplitude of the time-domain signal. (The amplitudes will not be exactly the same in the time-domain and frequency-domain because of ‘leakage’ due to computations on discrete signals.)

[Aaron Smith]

Ok, I think I understand. Overall, the code works really well and I can capture Peak 4 in Figure 2 with no problem (Fig4). Only one minor issue. Peak 2 in Fig2 is being blunted a little too much (Fig5) such that find peaks is not picking it up relative to what is around it (Fig4). I have tried tweaking the parameters but am getting large changes in the average y-int of the data (Fig3). Have to admit I don't really know what the parameters control very well.

[Star Strider]

I am not certain what parameters you are referring to.

You can improve the ‘detail’ of the filtered signal result by increasing the bandpass and bandstop frequencies of the filter. That may make it possible for findpeaks to isolate the peaks you want. You will likely have to experiment with it. (I have findpeaks, however I don’t understand the characteristics of the peaks you want to isolate, so I haven’t experimented with it with your data.)