High-pass filter

68 次查看（过去 30 天）

显示更早的评论

Guglielmo Giambartolomei 2021-5-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/828375-high-pass-filter

评论： Star Strider 2021-7-22

采纳的回答： Star Strider

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider 2021-5-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/828375-high-pass-filter#answer_698280

在 MATLAB Online 中打开

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

28 个评论
显示 26更早的评论隐藏 26更早的评论

William Rose 2021-5-13

在 MATLAB Online 中打开

@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.

William Rose 2021-5-14

@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."

Star Strider 2021-5-18

在 MATLAB Online 中打开

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 2021-7-22

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 2021-7-22

As always, my pleasure!

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

请先登录，再进行评论。

类别

Radar Phased Array System Toolbox Waveform Design and Signal Synthesis Matched Filter and Ambiguity Function

在 Help Center 和 File Exchange 中查找有关 Matched Filter and Ambiguity Function 的更多信息

产品

MATLAB

版本

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by