How to filter a low frequency movement from my data

Goodmorning,

For my thesis i'm analyzing wave data of flume tests. But it seems that there is some sort of 'long wave' present in the flume. Now i am trying to filter this long wave out of my data but i can't find a decent way to to this.

I already tried to filter my data using the Polyfit and polyval functions, this did filter the long wave out of the signal. But in my opinion this was more or less tweaking the function to get a good result. Secondly i do not want to use a moving avarage because this eliminates data. My supervisor told me to try the Filtfilt option. But now i am stuck wondering how to filter my signal with use of this FiltFilt.

See the figure below, there seems to be some longer wave present which i want to filter out for both the water level elevation and velocity. I was wondering if this is even possible or that my 'timeframe' is to short.

[Birdman]

Can you share this data in a mat file?

[Twevers]

Dear birdman,

added is the raw data of the signal. 1st colum is the time, 2nd column is the waterlevel elevation and the third is the velocity.

Answer 1

[Star Strider]

This uses a highpass Chebyshev Type II filter to eliminate the low frequency variation and the c-d (constant) offset from both signals. You can use the Fourier transform in figure(2) as a guide to fine-tuning the passband (here Wp) and stopband (here Ws) frequencies to eliminate the frequencies you want. (Since this is a highpass filter, Ws must be lower than Wp.) I designed a highpass filter rather than a bandpass filter since your signals seem to have very little high-frequency noise.

The Code —

D = load('Time elevation velocity.mat');
t = D.datavec(:,1);
wle = D.datavec(:,2);
vel = D.datavec(:,3);

figure(1)
plot(t, wle,    t, vel)
grid

Ts = mean(diff(t));                                         % Sampling Interval
Fs = 1/Ts;                                                  % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
L = size(t,1);                                              % Signal Length
FTwle = fft([wle vel])/L;                                   % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector

figure(2)
plot(Fv, abs(FTwle(Iv))*2)
grid
axis([0  3    ylim])

Wp = 0.25/Fn;                                               % Stopband Frequency (Normalised)
Ws = 0.20/Fn;                                               % Passband 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,'high');                           % Filter Design, Sepcify Bandstop
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability

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

wle_vel_filtered = filtfilt(sos, g, [wle vel]);             % Filter Signal

figure(4)
plot(t, wle_vel_filtered)
grid

[Twevers]

dear Star Strider,

many thanks for your help! I appreciate it alot.

Still i have two more questions regarding your answer.

1) In your explanation you say that Wp is the passband and Ws is the stopband. But in the script you say the opposite :-)

2) Figure 2 shows a clear peak at 0,42. As i am new to this kind of filtering i was wondering how to fine tune it such a way to eliminate this frequency. The reason is that i have more than 100 of these signals and each signal is different. As i am trying to understand the script such that i can fine tune each signal on its own with understanding what i did, instead of doing something to make the signal looks nice :-)

[Star Strider]

As always, my pleasure.

1) The explanation is correct. The code should be:

Wp = 0.25/Fn;                                               % Passband Frequency (Normalised)
Ws = 0.20/Fn;                                               % Stopband Frequency (Normalised)

2) You can change those values to:

Wp = 0.44/Fn;                                               % Passband Frequency (Normalised)
Ws = 0.43/Fn;                                               % Stopband Frequency (Normalised)

I chose the original values to preserve the relatively constant region between 0 and 12 (seconds). If you want to eliminate the frequencies between 0.22 and 0.45 Hz, change the respective assignments to:

Wp = [0.22 0.43]/Fn;                                       % Passband Frequency (Normalised)
Ws = [0.20 0.45]/Fn;                                       % Stopband Frequency (Normalised)

and:

[z,p,k] = cheby2(n,Rs,Ws,'stop');                           % Filter Design, Specify Bandstop

I am not quite sure what result you want. Feel free to experiment.

Answer 2

[Domanic]

Applying a filter with filtfilt removes the filter offset, so you'll probably find that combining filtfilt with a moving average (high pass) filter gives you the results that you're looking for.

Let's say that you want to filter out a sine wave oscillating at 1/10th the frequency of another sine wave. Then you could use:

X=0:0.01:10*pi;
sumsin = sin(X)+sin(10*X);
N=129; % The cut-off frequency -> 0 as N->inf
plot(X,sumsin)
ff = -ones(1,N)/N; % the 
ff(65)=1-ff(floor(N/2));
hold all;
plot(X,filtfilt(ff,1,sumsin))

Of course, better quality low pass filters exist, for example Butterworth ('butter' in Matlab), but moving average is a really easy place to start. If you want to visualise your filter, use

fvtool(ff,1)