8
$\begingroup$

I have troubles doing a simple high pass filtering...

Downloading the test data https://pastebin.com/74J8t2YV as "test.dat", you can see that overimposed on the signal there are some low frequencies and an almost constant slope, that I'd like to get rid of:

 ListPlot[Import[NotebookDirectory[] <> "test.dat", Joined -> True]]

The problem is that using the following code, to have the signal oscillate around a more-or-less constant line I need to use hf > 2 pi / fsamp, but at frequencies so high the important details of the signal are lost.

Module[{dat, a, times, lf, hf = 1, fl, filt, avgy, fsamp},
 dat = Import[NotebookDirectory[] <> "test.dat"];
 times = dat\[Transpose][[1]];
 a = dat\[Transpose][[2]];
 avgy = Mean[a];
 a = a /. y_ -> y - avgy;  (* not interested in signal offset *)

 fsamp = 1/Mean[Differences[times]]; (* sampling slightly irregular *)

  hf = 2 \[Pi] 0.9/fsamp;
 filt = HighpassFilter[a, hf];
 ListPlot[{{times, a}\[Transpose], {times, filt}\[Transpose]}, 
  Joined -> True]]

original and filtered data

I'm not an expert at all in frequency analysis, but I did these kind of stuff blindly using Microcal Origin - and there the highpass filtering does not suffer of this drawback... what can be done?

Here is the results in Origin (the signal is the same, I just forgot to subtract the average) Origin high pass filter in action

unfortunately the help page doesnt enter into details about the algorithm, where you only choose the cutoff frequency

$\endgroup$
  • $\begingroup$ I would recommend taking the first order discrete derivative of the data. $\endgroup$ – Ihsan Aug 9 '17 at 15:39
  • $\begingroup$ Relevant for constructing filters 'mathematica.stackexchange.com/a/153421/12558' $\endgroup$ – Hugh Aug 9 '17 at 17:05
  • $\begingroup$ The samples you provided don't have a constant time step. You may have to resample or just assume a constant time step (like the average time step, which is 0.0610). Furthermore, how about just fitting a low-degree polynomial and subtracting it to detrend the signal. I would subtract the mean first. Then the trend doesn't look quite linear so maybe fit a quadratic and subtract it. $\endgroup$ – Fixed Point Aug 9 '17 at 18:58
7
$\begingroup$

The default kernel length of HighpassFilter seems to be too small, after modifying it to Round[Length@a/10], the result is almost the same as that of Origin:

times = dat\[Transpose][[1]];
a = dat\[Transpose][[2]];
avgy = Mean[a];
a = a - avgy;
(* Alternative method for obtaining a: *)
(*
   a = Standardize[dat\[Transpose][[2]], Mean, 1 &];
 *)

fsamp = 1/Mean@Differences@times;(*sampling slightly irregular*)
hf = 0.3;(*Found by trial and error *)

filt = HighpassFilter[#, hf, Round[Length@#/10], SampleRate -> Round@fsamp] &@a;

ListLinePlot[{{times, a}\[Transpose], {times, filt}\[Transpose]}, PlotRange -> All]

Mathematica graphics

$\endgroup$
  • $\begingroup$ wonderful, thank you. I notice you added also SampleRate, which I didnt use $\endgroup$ – alessandro Aug 9 '17 at 12:23
  • 1
    $\begingroup$ @alessandro AFAIK, setting SampleRate is the standard way to tell HighpassFilter etc. how many sample is taken every one second. Your approach i.e. directly scaling the cutoff frequency is also valid, of course :) $\endgroup$ – xzczd Aug 9 '17 at 13:03
8
$\begingroup$

Here is my take on how to do filtering for a time series. The HighpassFilter you used is I think for image processing.

Import the data and find the sample rate

data = ToExpression[Partition[StringSplit[Import["74J8t2YV.txt"]], 2]];
ListLinePlot[data]
sr = 1/(data[[2, 1]] - data[[1, 1]])

Mathematica graphics

The sample rate is 16.9492 and from the data it looks like you have a main frequency at about 0.2 Hz. The data starts at about 3400 units. I am assuming you wish to get rid of the rising drift. Now check the frequency content by doing a Fourier transform.

ft = Fourier[data[[All, 2]], FourierParameters -> {-1, -1}];
ff = Table[(n - 1) sr/Length[ft], {n, Length[ft]}];
ListLogPlot[Transpose[{ff, Abs@ft}][[1 ;; 500]], PlotRange -> All, 
 Joined -> True]

Mathematica graphics

Looking at the first peak this confirms that there is a frequency of about 0.2 Hz.

Now I build a high pass filter. I start with an analogue Butterworth filter of 4th order and with cut off frequency fp which I set to 0.1Hz. You may wish to play with this value to see what frequency is best for you.

fp = 0.1; (* filter cut off frequency *)
filter = ToDiscreteTimeModel[ 
  ButterworthFilterModel[{"Highpass", 4, fp 2 π}], 1/sr];

Now apply the filter and plot the results.

fd = RecurrenceFilter[filter, data[[All, 2]] - data[[1, 2]]];
ListLinePlot[Transpose[{data[[All, 1]], fd}], PlotRange -> All]

Mathematica graphics

I have taken away the first point of the data so that the filter is not hit by a big step. This value could be adjusted a little and may take out the low frequency initial oscillation.

Hope that helps.

Edit

Following a comment from xzczd I have changed the code for generating the filter by removing division by the sample rate. This is correct there is no need to divide by the sample rate here. Unfortunately the original data has gone so I can't retest the code.

$\endgroup$
  • $\begingroup$ very detailed, thank you! I'll try to read it in detail to understand it better... $\endgroup$ – alessandro Aug 9 '17 at 12:53
  • 1
    $\begingroup$ "The HighpassFilter you used is I think for image processing. " Hmm… I think image processing, signal processing, time series processing, etc. have no essential difference. Actually most examples in the document of HighpassFilter are not related to image, and HighpassFilter is listed in the guide for time series processing. Anyway, it's interesting to know how to build a filter in a relatively low level way, +1 :) $\endgroup$ – xzczd Aug 9 '17 at 13:10
  • 1
    $\begingroup$ @xzczd If you read the post I linked to I analyse the Lowpass, Highpass and Bandpass filters. The Bandpass seems to do nothing. I agree that filtering is filtering with the purpose of removing frequencies from the data. However, if you want features designed-in like causality, and precise cut off values (I should have prewarped above for more precision) then there is a vast literature on this and it is good that it is supported by Mathematica. Thanks for the +1. $\endgroup$ – Hugh Aug 9 '17 at 13:25
  • $\begingroup$ +1 for considering investigating the fourier transform. $\endgroup$ – mathreadler Aug 9 '17 at 18:08
  • $\begingroup$ Just had a second look on the answer… Is the fp 2 π/sr part wrong actually? The value of fp is determined based on the fact that there's a peak at about 0.2 Hz in Periodogram, right? Then it should be fp 2 π because the cutoff frequency in ButterworthFilterModel is the "true" frequency of the signal, and sample rate already has its place in ToDiscreteTimeModel. (In this case the result won't have significant change of course. ) $\endgroup$ – xzczd Jun 21 '18 at 6:46
7
$\begingroup$

How about LowpassFiltering the signal and subtracting the result from the original?

ListLinePlot[
 Transpose[{
   dat[[All, 1]],
   # - LowpassFilter[#, 0.005] &@dat[[All, 2]]
  }]
 , AspectRatio -> 1/5
 , ImageSize -> 1000
 ]

enter image description here

$\endgroup$
  • $\begingroup$ it works this way too! Possibly LowpassFilter has a larger kernel length? $\endgroup$ – alessandro Aug 9 '17 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.