High-pass filter

Hello,

I'm working on acceleration signals coming from an accelerometer that has a sensitivity change of 0% from 2 Hz to 10000 Hz. When I perform the FFT of one this accelerometric signals I often obtain strong low frequency contents close to 0.1 Hz. I would like to cut away all the frequencies before 2 Hz because I don't think they are correct. The code I used is:

Fca=2;  %cutoff frequency

[z,p,k]=butter(8,Fca/(Fsa/2),'high');
sos=zp2sos(z,p,k);
%fvtool(sos,'Analysis','freq')

xafiltfilt=filtfilt(sos,g,xa);
%figure;
%plot(ta,xa,ta,xafiltfilt);
%grid on
%title('xa (accelerometric signal)')
%xlabel('t (s)')
%ylabel('Acceleration (g)')
%legend({'Original xa Signal','Actual xa Signal (filtered and realigned signal)'});
%saveas(gcf,'x_acc_filt_comp','jpg')


Xaf=fft(detrend(xafiltfilt),NFFT2_acc)/Na;
figure;
plot(fa,2*abs(Xaf(1:NFFT2_acc/2+1))) 
grid on
title('Single-Sided Amplitude Spectrum of xa(t) filtered and realigned Signal')
xlabel('f (Hz)')
ylabel('Xa_filtered(f)')

I was wondering if this type of filter may be adequate.

Thank you very much,

Guglielmo

Answer 1

[Star Strider]

It could work, however a better option could be:

Fca = 2;
Fsa = 1000;                                                                 % Create Parameter
xa = randn(1E+4,1);                                                         % Create Signal
[xafiltfilt,hpdf] = highpass(xa, Fca, Fsa, 'ImpulseResponse','iir');

figure
freqz(hpdf.Coefficients, 2^16, Fsa)
set(subplot(2,1,1), 'XLim',[0 10])                                          % Zoom Frequency Axis
set(subplot(2,1,2), 'XLim',[0 10])                                          % Zoom Frequency Axis

Compare to the Butterworth design:

[z,p,k]=butter(8,Fca/(Fsa/2),'high');
sos=zp2sos(z,p,k);

figure
freqz(sos, 2^16, Fsa)
set(subplot(2,1,1), 'XLim',[0 10])                                          % Zoom Frequency Axis
set(subplot(2,1,2), 'XLim',[0 10])                                          % Zoom Frequency Axis

.

[William Rose]

@Guglielmo Giambartolomei,

@Star Strider has a good idea. Since you use filtfilt(), the phase difference will be zero at all frequencies, which is a significant difference from @Star Strider's suggestion. Likewise, your attenuation (in dB) will be double at all frequencies, compared to if you had not used filtfilt(). I am always on the watch for filter instability with filtfilt(), when the initial filter order is greater than four. But if it works, great.

[Star Strider]

I do not understand that filtfilt is in any way different from what I suggested. The filtfilt funciton has its occasional problems, however I routinely advocate its use, since it is phase-neutral. My purpose here was to compare the Butterworth and elliptic filter designs.

[William Rose]

@Guglielmo Giambartolomei, @Star Strider is right that filtfilt() behaves as he described. I misinterpreted his suggestion to be advocating the use of an elliptic filter instead of "Butterworth plus filtfilt". In fact he is suggesting using "elliptic plus filtfilt" instead of "Butterworth plus filtfilt". filtfilt() does not appear in his code, but it is called internally by highpass(), which he calls. The transfer functions shown have nonzero phase, which ones does not expect when using filfilt(). That is because they are the transfer functions before filtfilt(); after filtfilt(), the phase will be zero. highpass() , with the arguments used, finds the lowest order IIR filter that meets the (default) specs. That turns out to be an 8th order elliptic filter. (You can find the order of hpdf with filtord(hpdf).) An elliptic filter has a much narrower transition width than a Butterworth of equal order - which is good. Before the advent of forward-backward filtering, people often avoided elliptic filters, despite their narrow transition width, because of their highly nonlinear phase, which causes a lot more ringing and overshoot than a Butterworth. But forward-backward filtering (which filtfilt() does) eliminates the phase delay, so there is more reason to use an elliptic filter.

Look at the effect of your filter in the time domain as well as the frequency domain. I have measured cranium accelerations during soccer heading. I have measured boot, shank, and sacrum accelerations during figure skater jump landings. In both cases, there is a sharp acceleration impulse which we wanted to characterize (peak, rate of change, etc.). The highpass filters we are comparing here, 8th order Butterworth+filtfilt and 8th order elliptical+filtfilt, alter the shape of a 100 millisecond square impulse a lot - see plots. Notice that the peak positive acceleration is a lot lower after filtering (although -peak to +peak is well preserved), and the filters cause significant ringing for more than 500 msec before and after. We would not have used either of these filters, due to the significant shape distortion which they cause.

Code to make figure above (same filters as @Star Strider):

>> x=[zeros(1,4950),ones(1,100),zeros(1,4950)];
>> fc=2; fs=1000;
>> [z,p,k]=butter(8,fc/(fs/2),'high');
>> [sos,g]=zp2sos(z,p,k);
>> yb=filtfilt(sos,g,x);
>> [ye,hpdf]=highpass(x,fc,fs,'ImpulseResponse','iir');
>> t=(0:length(x)-1)/fs;
>> figure; subplot(2,1,1);
>> plot(t,x,'k.-',t,yb,'r.-',t,ye,'b.-');
>> legend('Raw','Butter','Ellip.');
>> subplot(2,1,2);
>> plot(t,x,'k.-',t,yb,'r.-',t,ye,'b.-');
>> legend('Raw','Butter','Ellip.');
>> xlim([4.5,5.5]);

Good luck @Guglielmo Giambartolomei.

[Guglielmo Giambartolomei]

Thank you very much @Star Strider and @William Rose for your precious contribution! I carefully read what you wrote and I appreciate the effort you put in describing the behaviour of the two filters. However I didn't understand what are the substantial differences between an elliptic filter and a Butterworth filter. Why the former is more suitable for my application? I'm sorry but I'm not an expert and I would like to understand what I'm doing.

Thank you very much!

Best regards,

Guglielmo

[Star Strider]

@Guglielmo Giambartolomei — It depends on the result you want, and the characteristics of the signal. The Butterworth design is maximally flat in the passband, however it has a transition region that is much more gradual. The elliptic filter has ripple in both the stopband and passband (the pasband ripple is characteristically about 1 dB, so not large, however is much closer to the ‘ideal’ filter, with a very steep transition region.

@William Rose — Thank you for the interesting explanation!

Note that you can actually run your code onlilne here with the Run app (use the right-pointing green arrowhead in the toolbar) —

x=[zeros(1,4950),ones(1,100),zeros(1,4950)];
fc=2; fs=1000;
[z,p,k]=butter(8,fc/(fs/2),'high');
[sos,g]=zp2sos(z,p,k);
yb=filtfilt(sos,g,x);
[ye,hpdf]=highpass(x,fc,fs,'ImpulseResponse','iir');
t=(0:length(x)-1)/fs;
figure; subplot(2,1,1);
plot(t,x,'k.-',t,yb,'r.-',t,ye,'b.-');
legend('Raw','Butter','Ellip.');
subplot(2,1,2);
plot(t,x,'k.-',t,yb,'r.-',t,ye,'b.-');
legend('Raw','Butter','Ellip.');
xlim([4.5,5.5]);

.

[Guglielmo Giambartolomei]

Thank you again @Star Strider and @William Rose! I wouldn't want to abuse your kindness but I would like to sk what could be a good filter for the signals coming from this accelerometer. The sampling frequency was 600 Hz and it has a perfectly flat response (sensitivity change = 0%) between 2 Hz and 10000 Hz.

Thank you very much!

Best regards,

Guglielmo

[Star Strider]

For a perfectly flat response in the passband, use a Butterworth design.

The problem is this:

‘The sampling frequency was 600 Hz and it has a perfectly flat response (sensitivity change = 0%) between 2 Hz and 10000 Hz.’

With a sampling frequency of 600 Hz, the highest usable frequency is half that (the Nyquist frequency), 300 Hz.

This is one option for a highpass filter with a 2 Hz corner frequency —

Fs = 600;
Fn = Fs/2;
Wp = 2.0/Fn;
Ws = 1.5/Fn;
Rp =  1;
Rs = 60;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn,'high');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[0 10])                              % Zoom Frequency Axis
set(subplot(2,1,2), 'XLim',[0 10])                              % Zoom Frequency Axis

.

[William Rose]

@Guglielmo Giambartolomei, @Star Strider makes a very good point about the 600 Hz samplng frequency, which limits your spectrum to 300 Hz. In practice, I do not have great confidence in spectral estimates at frequencies higher than 1/4 of the sampling rate. This is because if you think about characterizing a sinusoidal wave - its frequency, ampltude, and phase - from only two or three samples in a period, you have a lot of room for error. If you have four or more.samples in a period, you have a better chance of getting a reasonably good estimate. Therefore, if my sampling rate is 600 Hz, I would have less confidence in the spectral estimates above 150 Hz.

As for what filter to use: it depends very much on what you want to do with the signal. It is very good that the response is flat from 2 to 10,000 Hz. However, flatness is not the same as lack of noise. There may still be noise over that range. If so, is the noise a problem for you? The answer will affect your choice of filter. You can record the noise spectrum while the accelerometer is at rest.

Test your candidate filter with simulated data to see how it affects the quantities you intend to measure, and choose a filter that you deem acceptable in that regard. Add simulated noise with a spectrum similar to what you measure when the accelerometer is quiet, if that noise is significant. In my previous posting, I simulated a 100 msec impulse and viewed it before and after filtering.

In that case, if my goal were to measure the low-to-high change in acceleration during an impulse, I would conclude that both filters are acceptable, and it doesn't matter which one I choose, because both filters preserve the height of the impulse, if measured from bottom to top. The conclusion "either filter is fine" is true surprisingly often, in my experience. I kind of feel guilty when I tell students that - but I remind them that they need to understand why. As I have gained experience (i.e as I have gotten old), I have become less particular about exactly what filter must be used, when I review a manuscript, as long as it is within the bounds of reason, for the purpose at hand.

If my goal were to measure pulse width, I would conclude that neither filter is acceptable, because they both cause great widening of the pulse. Although my example did not show it, the widening would have been just as great for a 1 msec or a 10 msec pulse as for a 100 msec pulse. Thereore both filters made it impossible to tell the difference between a 1 and a 10 and a 100 msec pulse.

If my goal were to measure total power within a frequency band, then I would choose a filter that is flat over that band - the Butterworth, as @Star Stridersaid.

If my goal were calculate someone's position by double-integration of their acceleration, I would just give up. All of my accelerometers have too much low-frequency drift. We tried it, for fun. We discovered that by the end of the practice, our skater (who was still in the rink) had moved by tens of kilometers, according to our double integration. Inertial navigation is hard to do well!

I'm sorry that I am not giving you a short simple answer to your question. As Blaise Pascal wrote (translated), "I have made this longer than usual because I have not had time to make it shorter."

[Guglielmo Giambartolomei]

Thank you @William Rose for trying to explain the situation to me in a funny way. I appreciated!

[Guglielmo Giambartolomei]

Dear @Star Strider, I am trying to use your precious suggestion (https://it.mathworks.com/matlabcentral/answers/828375-high-pass-filter#comment_1518918) but even if I read the help for "buttord" I did not understand the sintax. In particular, where is the corner frequency of 2 Hz? It seems to be W_s but I don't see it. Could you please quote all the passages in order to let me understand them?

Thank you very much!

Best regards,

Guglielmo

[Guglielmo Giambartolomei]

Dear @Star Strider,

I use the code you wrote:

Fs = 600;
Fn = Fs/2;
Wp = 2.0/Fn;
Ws = 1.5/Fn;
Rp =  1;
Rs = 60;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn,'high');
[sos,g] = zp2sos(z,p,k);

ya_filt_CT4H24X1_start = filtfilt(sos,g,ya_CT4H24X1_start );

figure(181);
plot(ta_CT4H24X1_start,ya_CT4H24X1_start,ta_CT4H24X1_start,ya_filt_CT4H24X1_start);
grid on
title('CT4H24X1 (8.23÷30 kPa)')
xlabel('t (s)')
ylabel('Acceleration (g)')
legend({'Original Signal','Actual Signal (filtered and realigned signal)'});
saveas(gcf,'CT4H24X1_start_filt_comp','jpg')

I plotted the 2 signals (original and filtered) and it seems there are some problems:

Where did I go wrong?

Thank you!

Regards,

Guglielmo

[Star Strider]

I did not see your earlier comment until just now. There are several terms for filter parameters. The ‘corner’ or ‘passband’ frequency is ‘Ws’ in the documentation.

Replying to your later Comment —

I doubt anything is actually wrong. The filter is stable, and highpass filters amplify noise, and this may be what you are seeing. Decreasing the stopband frequency seems to work to eliminate the ‘ringing’ at the end. Try that.

Using that code with a synthesized signal —

Fs = 600;
Fn = Fs/2;
Wp = 2.0/Fn;
Ws = 1.0/Fn;                                                    % Decrease 'Ws' From 1.5 Hz To 1.0 Hz
Rp =  1;
Rs = 60;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn,'high');
[sos,g] = zp2sos(z,p,k);

ta_CT4H24X1_start = linspace(0, 33*600, 33*600)/Fs;             % Create Time Vector
ya_CT4H24X1_start = randn(1, numel(ta_CT4H24X1_start))*0.05;    % Create Signal Vector

ya_filt_CT4H24X1_start = filtfilt(sos,g,ya_CT4H24X1_start );

figure(181);
plot(ta_CT4H24X1_start,ya_CT4H24X1_start,ta_CT4H24X1_start,ya_filt_CT4H24X1_start);
grid on
title('CT4H24X1 (8.23÷30 kPa)')
xlabel('t (s)')
ylabel('Acceleration (g)')
legend({'Original Signal','Actual Signal (filtered and realigned signal)'});
saveas(gcf,'CT4H24X1_start_filt_comp','jpg')

.

[Guglielmo Giambartolomei]

You are very kind and I am grateful to you!

I changed Ws = 1.5/Fn with Ws = 1.0/Fn and I found a more coherent signal:

It seems that the applied filter aligned the signal to zero (it seems it shifted the filter towards lower values of the acceleration). What can I say in physical terms about this filter? What does it do?

Now I want to try to perform some FFTs of the filtered signal.

Thank you so much!

Best regards,

Guglielmo

[Star Strider]

Thank you!

As always, my pleasure!

‘It seems that the applied filter aligned the signal to zero ...’

A highpass filter will eliminate any constant offset, and baseline variations below the cutoff frerquency.

‘What can I say in physical terms about this filter? What does it do?’

That is difficult to describe. A highpass filter passes only the frequencies above the cutoff frequency, up to the Nyquist frequency (for a discrete filter) or infinity (for an analog filter). Highpass filters can also amplify noise, as this one also appears to be doing, at least in relation to the frequencies it eliminates.

[Guglielmo Giambartolomei]

Hi again, just for curiosity I performed the FFTs of the original signal and of the filtered one:

Ya_filt_CT4H24X1_start=fft(detrend(ya_filt_CT4H24X1_start),NFFT2_CT4H24X1_acc_start)/Na_CT4H24X1_start;

figure(182);
plot(fa_CT4H24X1_start,2*abs(Ya_CT4H24X1_start(1:NFFT2_CT4H24X1_acc_start/2+1)),fa_CT4H24X1_start,2*abs(Ya_filt_CT4H24X1_start(1:NFFT2_CT4H24X1_acc_start/2+1)));
grid on
title('FFT of the CT4H24X1 (8.23÷30 kPa) acceleration signal')
xlabel('Frequency (Hz)')
ylabel('Amplitude(f)')
legend({'Original Signal','Actual Signal (filtered and realigned signal)'});
saveas(gcf,'CT4H24X1_start_filt_comp_FFT','jpg')

And the result seems to be good:

After about 2 Hz the contents in frequency are almost the same (the 2 plots are overlapped), while before 2 Hz there are no contents (orange plot). I hided all the unknow low frequency contents. I would like to know to what physically correspond those low frequency contents but I think I can't.

Thank you very very much!!!

You are very kind!

Best regards,

Guglielmo

[Star Strider]

The low frequencies are variations in the acceleration that exist at those frequencies. I have no idea what the acceleration signal actually represents (what acceleration it actually measures or how the accelerometer is oriented), so I cannot suggest the physical process that could have caused them. So just guessing, they could be starting or stopping transients (that would be expected to be gradual and therefore slow and of low frequency, for example with a vehicle), while the other parts of the signal would be accelerations caused by other influences of higher frequencies (such as road roughness, wind variations affecting light aircraft, or rough sea conditions).

Also, I thought a bit more about physical analogs to various filters, and for a highpass filter, that would be a capacitor in series, an inductor in parallel, or mechanically, perhaps a spring. A lowpass filter would be an inductor in series, a capacitor in parallel, or mechanically a mass. (My background is primarily electrical engineering, so there could be more appropriate mechanical analogs.)

[Guglielmo Giambartolomei]

Given that before 2 Hz the accelerometer sensitivity changes a lot, couldn't it be noise?

Thank you!

[Star Strider]

As always, my pleasure!

It certainly could be noise. The ‘Sensitivity Change’ may not reflect the device transfer function (it could be proportional to the derivative of the transfer function), so (given that I have essentially no experience with accelerometers in my own work) I have no idea.

[Guglielmo Giambartolomei]

Thank you again!!!

[Star Strider]

As always, my pleasure!

[Guglielmo Giambartolomei]

Good morning @Star Strider, sorry to bother you again, I have a doubt. With respect to the code:

Fs = 600;
Fn = Fs/2;
Wp = 2.0/Fn;
Ws = 1.0/Fn;
Rp =  1;
Rs = 60;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn,'high');
[sos,g] = zp2sos(z,p,k);

Ws is the stopband corner frequency and it should be a scalar|two elements vector. We used a scalar (0.0333). I can't understand how this scalar is linked to 2 Hz. Then I can't understand why the parameters Wp and Ws are related to Fn.

Thank you very much,

best regards,

Guglielmo

[Star Strider]

I do not understand your question.

The passband and stopband frequencies are scaled to the Nyquist frequency (‘Fn’), since that is the highest identifiable frequency (equal to π radians) in a sampled signal, and must be on the interval [0,1]. All frequencies must be scaled this way.

Please see the documentation for butter for a full explanation. (The buttord function calculates the optimal order for the filter.)

.

[Guglielmo Giambartolomei]

Good morning and thank you for your reply @Star Strider. I read the documentation but I still can't figure out how to set the stop band corner frequency. We set it at two Hz in the previous example and Ws turned out to be 0.0333. If i wanted to change it to 3 Hz how could I do?

Thank you very much,

best regards,

Guglielmo

[Star Strider]

The corner frequency is correct.

The 2 Hz design frequency becomes 0.0333 π radians/s when normalised by the Nyquist frequency.

.

[Guglielmo Giambartolomei]

I'm sorry @Star Strider but I did not understand. Is there a formula to calculate the stopband corner frequency normalised by the Nyquist frequency? In simple words, if I have a Fs (sampling frequency) and I want a stop band corner frequncy of "t" Hz, how can I write Ws as a function of Fs/2=Fn?

Sorry to bother you but I really would like to understand this passage.

Thank you very very much,

Guglielmo

[Star Strider]

‘Is there a formula to calculate the stopband corner frequency normalised by the Nyquist frequency?’

Yes.

‘In simple words, if I have a Fs (sampling frequency) and I want a stop band corner frequncy of "t" Hz, how can I write Ws as a function of Fs/2=Fn?’

Fs = ...;                                       % Sampling Frequency (Hz)
Fn = Fs/2;                                      % Nyquist Frequency (Hz)
t = ...;                                        % Stopband Frequency (Hz)
Ws = t/Fn;                                      % Stopband Frequency (Normalised)

No worries!

I hope this clears up any confusion over those calcualations.

.

[Guglielmo Giambartolomei]

Very easy! You are the best! I will always be grateful for your help! Thank you very much!

I was confused because in the previous example we wrote Ws=1.0/Fn and we said that the stopband corner frequency was 2 Hz. As a matter of fact is 1 Hz.

Best regards,

Guglielmo

[Star Strider]

As always, my pleasure!

I do not see where I typed that, however I do not always catch all my typographical errors.

.