1
$\begingroup$

I'm looking to a way to detect peaks in Matlab for noisy random signals generated by IMU's accelerometers. The signal could be approximate to a sinusoidal-like signal with zero mean and random noise on it. A lot of of noise.

I don't want to remove noise using pass-band filters to not add non idealities but the only use of findpeaks seams not sufficient. I'd like to remain in the time domain and act on RAW data.

Looking at documentation of findpeaks() I can set some parameters to modify the window it detects. In particular:

  • MinPeakHeight, I could use it to avoid to detect low amplitude signal. I was thinking something computed like mean(signal) + 2*std
  • MinPeakProminence, to avoid false positive of peaks too close to each other with high peaks. 1.5*std?
  • MinPeakDistance, as 0.005*Fs, the sampling frequency of my IMU (100 Hz).

But I think this is not so robust for at least three problems: mean close to zero (1e-6), code not robust if signal periodicity changes (this is the main one), SNR to low in some cases. Looking at the RAW signal shape I can tell visually where are the peaks. But is like doing a sort of filtering and taking a moving average of the signal. I don't want to do it for now (if possible).

[peaks, locs] = findpeaks(signal, time, 'MinPeakHeight',hThr, 'MinPeakProminence',pThr,'MinPeakDistance',dPts);

This is an example of signal. Any way to do the job? Increase the SNR? How?

enter image description here

Blue is RAW, red is a linear detrended version to (try) remove IMU's drift.

Here the CSV signal with time and RAW data. A third column with a moving average to visually see the peaks.

enter image description here

Improve this question
$\endgroup$
8
  • $\begingroup$ Welcome to SE.SP! What do you mean by linear detrended version ? The red and blue plots appear to just differ by a constant. Matlab's detrend operation usually removes a linear (or possibly polynomial) trend. $\endgroup$ Commented Jun 16 at 11:12
  • 1
    $\begingroup$ Precisely what you said. I just did detrend(signal,1) to obtain the red curve. I used just to check if the shift in amplitude of the blue curve was due to drift in the IMU or not. And I guess so (they told me these IMUs have 10° drift over 1 hour sensing...a lot...) $\endgroup$ Commented Jun 16 at 11:56
  • $\begingroup$ @Shika93, Can you share the signal as CSV? Could you mark some peaks you'd think are good? If you do, I can try some algorithms on it. $\endgroup$ Commented Jun 16 at 12:50
  • $\begingroup$ OK, thanks for confirmation. The MinPeakHeight setting is not going to work well for a signal with a wandering baseline. I'm wondering if Threshold might be better to set to something like std? $\endgroup$ Commented Jun 16 at 12:51
  • 1
    $\begingroup$ In order to increase your SNR, you need to characterize your noise so that you can reduce its presence. We use the frequency domain to help perform this characterization because we often have some assumptions about what our signal looks like (e.g., narrow-band/sinusoidal). Time varying properties of the signal's frequency informs what window length to use in the analysis and whether to use a single FFT or a STFT. It's perfectly reasonable to use a spectral domain to better inform what time domain operations you should use, or how to develop a filter that minimizes signal distortion. $\endgroup$ Commented Jun 16 at 15:05

2 Answers 2

Reset to default
2
$\begingroup$

Assuming you wan to have the Oscillatory component decomposed from the trend you need to have a decomposition model which fits this assumption.

It could probably be achieved with classic decompositions such as DFT and Wavelets.

Yet there are specialized decompositions for this case.
2 common ones are:

  • STL: A Seasonal Trend Decomposition Procedure Based on Loess.
  • Singular Spectrum Analysis for Time Series.

I used the STL to decompose the signal to many periodic components which results in a good trend:

enter image description here

I took the trend component and removed it from the original signal which yielded:

enter image description here

On the DeTrend signal one could do the peak analysis.

Signal decomposition is relatively "hot" use case. There is even an approach based on Deep Learning (Using a Transformer), See: Additive Decomposition of One Dimensional Signals using Transformers.

Remark: With tuning of the decomposition you may get the periodic component you're after directly. Probably can be done using DFT / STFT as well.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Detrending the signal in this way or calculating a moving average and subtract it from the RAW is different? Because I thought to detrend through movmean and find the peaks though islocalmax instead of `findpeaks. $\endgroup$ Commented Jun 21 at 12:41
  • $\begingroup$ @Shika93, I thought you tried moving average, did you get the same result? $\endgroup$ Commented Jun 21 at 18:37
  • 1
    $\begingroup$ @Shika93, At the end, you need some kind of abstract High Pass Filter. The question whether it is defined only by frequency (Classic HPF) or adaptive to data which some decompositions allow. $\endgroup$ Commented Jun 22 at 5:20
0
$\begingroup$

You need a MOVING AVERAGE FILTER, not a conventional filter.

You have only acquired what looks like 1/4 of a cycle only, while taking 1 time sample per second, what apparently looks like many time samples, but the sampling frequency is not correct.

Actually there are not enough time samples because the sampling frequency is too low compared to the frequency of the sought tone.

It doesn't matter taking thousands of samples if they are all constrained to 1/4 of the basic time cycle of the sought signal.

You correctly point out that a conventional filter is not going to help.

Because the sampling frequency is way below the bare minimum of 2 samples per cycle, the smallest available frequency with fft is going to be the basic frequency step df, not the sought frequency, which is way below df.

findpeaks on the time signal doesn't work, too much noise.

You will spend considerable time attempting to set up the right parameters of findpeaks to find out that a simple frequency shift in the input signal, more noise, less noise, will render the selected parameters pointless and you will have to start all over again trial and error attempting to configure findpeaks.

But findpeaks is in this context useful along the frequency domain.

1.- CONVENTIONAL FILTERING DOES NOT WORK

Following an example of a basic moving average filter with ma=100 and ma=300 , ma being half the size of the time window used to shift along the available time samples

clear all; close all;clc % T=readtable('s1.csv'); % save T.mat T; load('T.mat') t=T.Time_s_; s=T.RAWSignal; L=numel(t); dt=(abs(t(1)-t(end))/L); Fs=1/dt; % s=s-mean(s); figure(1) plot(t,s) grid on title('s(t)') xlabel('t [s]')

enter image description here

abS=abs(fft(s)); angS=angle(fft(s)); abS_ = abS/L; abS1 = abS_(1:L/2+1); abS1(2:end-1) = 2*abS1(2:end-1); f = Fs*(0:(L/2))/L; % f1=f(f>.5); % DO NOT remove LOW FREQ components % n1=find(f==f1(1)); % S1=abS1(n1:end); % find Spectrum max peak [maxS1,nfmax] = findpeaks(abS1,'NPeaks',1,'SortStr','descend'); % Expected tone of the signal fmax=f(nfmax) % [Hz] T=1/fmax % [s] figure(2) plot(f,abS1) grid on title('S(f)') xlabel('f') hold on plot(fmax,maxS1,'or') plot([fmax fmax],[0 maxS1+1]) hold off

enter image description here

enter image description here

% frequency gating, instead of applying standard filters % abS1(1:n1-1)=0; % abS1(nfmax+10)=0; % abS1(1:n1-3)=0; % abS1(nfmax+3)=0; abS1(1:nfmax-1)=0; abS1(nfmax+1)=0; [S2re,S2im]=pol2cart(angS(1:L/2+1),abS1); S2=S2re+1j*S2im; s2=L*ifft(S2,numel(t)); figure(3) plot(t,s-s2) grid on title('s2(t)') xlabel('t [s]');

enter image description here

The apparent fundtamental frequency with the marker on the spectrum is in fact the 1st harmonic.

2.- MOVING AVERAGE FILTER

ma=[100 300] for k1=1:1:numel(ma) s2=[]; for k=ma(k1):1:L-ma(k1) s2=[s2 sum(s(k-ma(k1)+1:1:k+ma(k1)))/(2*ma(k1))]; end fgk=figure(k1+3) plot(s2) title(['moving average window : ' num2str(ma(k1)) ' samples']); grid on xlabel('t [s]') saveas(fgk,['00' num2str(k1+3) '.jpg']) end

enter image description here

enter image description here

This is the the accelerometer signal.

3.- NEITHER DECIMATING NOR INCREASING FREQUENCY POINTS HELP

Because too many time samples on a quarter of cycle not not a full cycle acquired, removing time samples does not help either :

nt=numel(t); % amount f time samples

default fft uses same amount of frequency points as input time samples

log2(nt)

% let's jump to, not the next power of 2, that'd be 2^13, let's make it 2^15

% floor(log2(nt))+3 n3 = floor(log2(nt))+3 nf=2^n3; abS=abs(fft(s,nf)); angS=angle(fft(s,nf)); abS_ = abS/nf; abS1 = abS_(1:nf/2+1); abS1(2:end-1) = 2*abS1(2:end-1); f = Fs*(0:(nf/2))/nf; f1=f(f>.5); % remove DC n1=find(f==f1(1)); S1=abS1(n1:end); % find Spectrum max peak [maxS1,nfmax] = findpeaks(S1,'NPeaks',1,'SortStr','descend'); % Expected tone of the signal fmax=f1(nfmax) % [Hz] T=1/fmax % [s] figure(6) plot(f1,S1) grid on title('S(f)') xlabel('f') hold on plot(fmax,maxS1,'or') plot([fmax fmax],[0 maxS1+1]) hold off

enter image description here

% frequency gating, instead of applying standard filters % abS1(1:n1-1)=0; % abS1(nfmax+10)=0; % abS1(1:n1-3)=0; % abS1(nfmax+3)=0; abS1(1:nfmax-1)=0; abS1(nfmax+1)=0; [S2re,S2im]=pol2cart(angS(1:nf/2+1),abS1); S2=S2re+1j*S2im; s2=L*ifft(S2,numel(t)); figure(7) plot(t,s-s2) grid on title('s2(t)') xlabel('t [s]');

enter image description here

And this is decimating the available time samples. It doesn't work either :

enter image description here

4.- DECIMATION LEADING TO 1ST HARMONIC

Decimating the input signal is tempting: Less samples, less processing required. And as shown in the following lines, it may lead to an apparently perfect tone, the 1st harmonic, that again may be considered as the fundamental.

% preceeding not shown .. nd=17; s=decimate(s,nd); t=decimate(t,nd); L=numel(t); figure(6) plot(t,s) grid on title('s(t)') xlabel('t [s]') abS=abs(fft(s)); angS=angle(fft(s)); abS_ = abS/L; abS1 = abS_(1:L/2+1); abS1(2:end-1) = 2*abS1(2:end-1); f = Fs*(0:(L/2))/L; f1=f(f>.5); % remove DC n1=find(f==f1(1)); S1=abS1(n1:end); % find Spectrum max peak [maxS1,nfmax] = findpeaks(S1,'NPeaks',1,'SortStr','descend'); % This is the main tone of the signal fmax=f1(nfmax) % [Hz] T=1/fmax % [s] figure(7) plot(f1,S1) grid on title('S(f)') xlabel('f') hold on plot(fmax,maxS1,'or') plot([fmax fmax],[0 maxS1+1]) hold off % frequency gating, instead of applying standard filters % abS1(1:n1-1)=0; % abS1(nfmax+10)=0; % abS1(1:n1-3)=0; % abS1(nfmax+3)=0; abS1(1:nfmax-10)=0; abS1(nfmax+10)=0; [S2re,S2im]=pol2cart(angS(1:L/2+1),abS1); S2=S2re+1j*S2im; s2=L*ifft(S2,numel(t)); figure(8) plot(t,s-s2) grid on title('s2(t)') xlabel('t [s]');

enter image description here

This is NOT the accelerometer signal.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.