Convert stepped signal to a continuous signal

Hi, I'm recording a voltage output at a sampling frequency of 1000Hz from a device that measures liquid flow. The refresh rate of the output is 125Hz. The result is a stepped curve as shown below.

In the end, I'm interested in the steepest slope of the curve. What should I do to construct the original curve from the recorded data? I've had a look at decimate(data, 8), and y=movmean(data, 8) followed by downsample(y, 8). Both ways result in better looking graphs, but which one would be the most accurate? I'm new to the field, so maybe there's a better way of doing things.

Answer 1

[Star Strider]

There are two ways I would consider.

The easiest is to use the Signal Processing Toolbox Savitzky-Golay filter function sgolayfilt. It creates what is essentially a low-pass filter from your data.

A somewhat more involved but more general approach is to take the Fourier transform of your signal (with the fft (link) function), visually determine signal frequency and the noise frequencies, and create a lowpass filter (or bandpass filter if you also want to remove baseline offset and low-frequency baseline variations). It is easiest to use the Signal Processing Toolbox designfilt function for this. Then filter your signal using the filtfilt function.

[RJ]

Thank you for your response.

I've been trying different things and reading about filters, but I'm not sure how to continue.

I generated the following spectrum from my data:

At 1.5 Hz, 3 Hz, 4.5 Hz ... I get some peaks. The 1.5 Hz makes sense to me because the flow I'm measuring pulses at about that frequency. Not sure about the other frequencies...

I created the following filter:

   Specifications:
      FrequencyResponse: 'lowpass'
        ImpulseResponse: 'iir'
             SampleRate: 1000
    StopbandAttenuation: 60
      StopbandFrequency: 12
      PassbandFrequency: 10
         PassbandRipple: 1
           DesignMethod: 'cheby2'

It actually produces a nice looking graph, but I have to admit I almost randomly picked these settings. What would be a more scientifically accurate way to create this filter?

[Star Strider]

My pleasure.

The Chebyshev Type II design is an excellent choice. These filters are best for low-frequency signals where the steepness of the stopband rolloff is a significant constraint.

There is probably no specific programmatic way to choose the filter parameters. Signal processing is heuristic by definition, so you need to experiment to see what works best.

Looking at the spectrum, the fundamental signal is from 1 Hz to 2 Hz. The others appear to be harmonics (at approximately integer multiples of the fundamental frequency), and all are likely modulated by the sampling process and other system nonlinearities, so there is a progressive ‘spread’ in the bandwidth of the harmonics. (I don’t have your signal to experiment with.)

I would put the lowpass filter passband frequency at 1.75 Hz and the stopband frequency at 2 Hz, since this includes the fundamental frequency, the d-c (baseline) offset (if those are important), and any low-frequency baseline drift.

If you want to eliminate the d-c offset and low-frequency baseline drift, I would use a bandpass design with the low-frequency stopband at 1 Hz, low-frequency passband at 1.25 Hz, and the corresponding high-frequency parameters at 1.75 Hz and 2 Hz.

Adjust the frequencies to get a stable filter that produces the result you want.

[RJ]

When I use the suggested frequencies in a lowpassfilter, I lose the 'shoulder' of the graph, which is a 'real' thing.

When I create a bandpass filter using the suggested values, I get the following graph (which does not make sense to me):

I have added the data if you'd like to have a look (x=time, y=signal).

[Star Strider]

It helps to have your data.

I’m not certain what result you want. I kept increasing the filter cutoffs until I got something that looks appropriate. The strategy I used was to keep the stopband frequency at least 0.5 Hz higher than the passband, then just incremented. You can easily use the same idea with designfilt, and keep increasing them until you get the result you want. The code I used here made it easier to change the filter parameters. Note that for a Chebyshev Type II design, the passband ripple can be anything because the passband is actually flat. Designing with a higher passband ripple will create a more stable filter and allow a steeper passband-to-stopband rolloff.

The Code:

D = load('RJ data.mat');
t = D.x(:);
s = D.y(:);

Ts = mean(diff(t));
Fs = 1/Ts;
Fn = Fs/2;

Wp = 13.60/Fn;
Ws = 14.10/Fn;
Rp = 10;
Rs = 30;
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);
[b,a] = cheby2(n,Rs,Ws);
[SOS,G] = tf2sos(b,a);

figure(1)
freqz(SOS, 2^15, Fs)
set(subplot(2,1,1), 'XLim',[0 15])
set(subplot(2,1,2), 'XLim',[0 15])

sf = filtfilt(SOS,G,s);

figure(2)
plot(t, s)
hold on
plot(t, sf)
hold off
grid
axis([1.5  2.0    ylim])
legend('Original', 'Filtered')

The Plot:

[RJ]

Thank you again!

However, your code doesn't seem to work when I copy-paste it into Matlab. I'm getting numbers for sf in the 10^165!

How would things work for a bandpass-filter? I'm unable to get that to work with the designfilt-option. Maybe it's related to the issue I'm experiencing with your code as well (using Matlab 2016 here).

To give you a better understanding what I'm trying to do, have a look at this example: Doppler spectrum

The signal I record is an envelop of the doppler spectrum and I want to determine the (averaged) slope between 1 and 2 in Matlab.

[Star Strider]

My pleasure!

I abandoned the Chebyshev IIR design and went with a FIR filter in order to get the passband and stopband I want. The filter works.

Since you want to calculate the slopes, I used the findpeaks function on the filtered data to identify the peaks and valleys, then subtracted adjacent values of both vectors to get a version of the approximate slopes. I will leave that exact calculation to you. I plotted a histogram of the approximate slopes my code calculates to give you an idea of the range and distribution. You will probably want to tweak my code to produce the result you want.

The Code:

D = load('RJ data.mat');
t = D.x(:);
s = D.y(:);

Ts = mean(diff(t));
Fs = 1/Ts;
Fn = Fs/2;

fcuts = [0.5 1.0 13.5 14.5];                                    % Frequency Vector
mags =   [0 1 0];                                               % Magnitude (Defines Passbands & Stopbands)
devs = [0.05  0.01  0.05];                                      % Allowable Deviations
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);
n = n + rem(n,2);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'scale');

figure(1)
freqz(hh, 1, 2^15, Fs)
set(subplot(2,1,1), 'XLim',[0 15])
set(subplot(2,1,2), 'XLim',[0 15])

sf = filtfilt(hh,1,s);

figure(2)
plot(t, s)
hold on
plot(t, sf)
hold off
grid
axis([1.5  2.0    ylim])
legend('Original', 'Filtered')

figure(3)
plot(t, -s)
hold on
plot(t, -sf)
hold off
grid
% axis([1.5  2.0    ylim])
axis([1.5  12.0    ylim])
legend('Original', 'Filtered')

[pks,pkt] = findpeaks( sf, t, 'MinPeakHeight',25, 'MinPeakDist',0.7);
[vls,vlt] = findpeaks(-sf, t, 'MinPeakHeight',25, 'MinPeakDist',0.7);
maxlen = min(length(pks),length(vls));

apx_slope = (pks(1:maxlen)+vls(1:maxlen))./(pkt(1:maxlen)-vlt(1:maxlen));
apx_slope = apx_slope(apx_slope > 0);

figure(4)
plot(t, sf)
hold on
plot(pkt,  pks, '^r', 'MarkerFaceColor','r')
plot(vlt, -vls, 'vr', 'MarkerFaceColor','r')
hold off
grid

figure(5)
histogram(apx_slope, 25)
grid
xlabel('Approximate Slope')
ylabel('Count')

ts = pkt(apx_slope > 0);

figure(6)
plot(t, sf)
hold on
plot(ts, apx_slope)
hold off
grid
legend('Filtered Signal', 'Approximate Slope')

[RJ]

I'll have a look at your suggestion. It looks very promising! Thank you kindly :)

[Star Strider]

As always, my pleasure!

[RJ]

Dear Star Strider, Hope you don't mind, but I've got a follow-up question. I'm also interested in the derivative of the signal. I've used your filter with the following settings.

fcuts = [1.0 13.5 14.5];                                    % Frequency Vector
mags =   [1 0];                                               % Magnitude (Defines Passbands & Stopbands)
devs = [0.01  0.05];                                      % Allowable Deviations

I've split my signal in separate curves and subsequently used Matlab's diff()-function. However, the derivative is not smooth at all and has a lot of distortions. The filter seems to be the culprit. Do you have a suggestion to obtain a smooth derivative curve as well?

Below you can see 20 isolated curves superimposed to show you what I mean. I'd like to fix each curve.

[Star Strider]

I’m not certain what the problem could be with the derivative itself. I doubt the filter itself is the problem. I would use the gradient function rather than diff, since it uses a much more sophisticated algorithm than simply the element-wise differences such as diff returns.

Note that taking the derivative of any ‘real world’ signal will amplify the noise. If you want to use the derivative, you may have to re-filter the the derivative of the signal — with the same filter — to get rid of the noise. The spectra of the original signal and the derivative should be the same, so the same filter should work well.

EDIT — I posted my last Comment before your plot posted. My Comment remains unchanged. You have high-frequency noise that my original filter will eliminate.