Butterworth filtfilt and fft/ifft problem

I have to filter strong motion data using a bandpass n=4 butterworth filter with cut-off frequencies of 0.01 and 40Hz. I have the original time series x(t) which I included in a txt file below. However I'm new with signal processing and I'm not sure if my filter is right and the fft functions are used correctly. Here is the code I have :

 Fs=1/Ts; %200Hz 
fn=1/(2*Ts);
fhc=0.8*fn/(Fs/2); %40Hz
flc=0.01/(Fs/2); %0.01Hz
Wn=[flc fhc];
f=(0:Nf-1)*Fs/Nf;
X=fft(x,Nf);
[b,a]=butter(1,Wn,'bandpass');
Y=filtfilt(b,a,X);
y=ifft(Y,'symmetric');
figure()
loglog(f,abs(X*Ts),'b',f,abs(Y3*Ts),'r');
hfvt=fvtool(b,a,'FrequencyScale','log','Fs',Fs);

When I plot the fourier's amplitude spectrum I get a downward shift in the spectrum. I'm not sure why it does that as the magnitudes plot using fvtool in between the cut-off frequencies are equal to 0. See below :

<<

>>

Thanks in advance.

Answer 1

[Star Strider]

First, you are specifying an order 1 filter, not an order 4 filter.

Second, your unfiltered signal appears only to have an upper frequency limit of 40 Hz, so your filter is only going to affect its lower frequencies. Because of that, I plotter the unfiltered and filtered signals separately, since they overplot each other.

I would do this:

filename = 'acc time series.txt';
fidi = fopen(filename);
D = textscan(fidi, '%s%s');
t = str2double(strrep(D{1}, ',','.'));                          % Replace Decimal Separator
X = str2double(strrep(D{2}, ',','.'));                          % Replace Decimal Separator
L = numel(t);

Ts = mean(diff(t));
Fs = 1/Ts;
Fn = Fs/2;
fhc = 40/Fn;
flc = 0.01/Fn;
Wn=[flc fhc];
n = 4;
[b,a] = butter(n,Wn,'bandpass');

Y = filtfilt(b, a, X);

S = [X Y];
FTS = fft(S)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
subplot(2,1,1)
plot(Fv, abs(FTS(Iv,1))*2)
grid
set(gca, 'XScale','log')
title('Unfiltered signal')
subplot(2,1,2)
plot(Fv, abs(FTS(Iv,2))*2)
grid
set(gca, 'XScale','log')
title('Filtered signal')

figure
freqz(b, a, 2^14, Fs)
set(subplot(2,1,1), 'XScale','log')
set(subplot(2,1,2), 'XScale','log')
title(subplot(2,1,1),'Filter Bode Plot')

Experiment to get the result you want.

[Eliza M.]

Great! Thank you. I put n=1 because I thought that the bandpass filter and filtfilt() each generated a 2n order effect.

Otherwise, the only problem I'm still having is that I have the FFT of the unfiltered signal computed from another software and when I compare it with X=fft(x,L), the spectrums are different. Note that with my original code the 2 unfiltered signal computed were identical. Any idea why this could happen? See below :

[Star Strider]

My pleasure.

The filtfilt function creates an effective doubling of the order because it is a forward-backward filter. I believe it then corrects for that in the output.

The second argument to fft is the length of the Fourier transform, and since it defaults to the length of the signal, ‘L’ has no effect here:

X=fft(x,L);

I have no idea what you are doing or how you created the third spectrum, so I cannot address that problem.

[Eliza M.]

Ok thanks.

I experienced with the version of your code and I'm still having some issues. I'm wondering if the filter seems right. When I use the ifft to convert it back in the time domain I get peaks at the beginning/end, the amplitude of acc is significantly reduced and at 150s it's almost zero and than goes back up a little. I tried multiple formulations and always seem to get this response as well as with other signals. Do you know what could be the cause of that? Here's my code :

x=Af;
t=T;
Fs=1/Ts;
Fn=1/(2*Ts);
fhc = 0.8*Fn/Fn; 
flc = 0.01/Fn;
Wn=[flc fhc];
n=2;
X=fft(x);
[b,a] = butter(n,Wn,'bandpass');
Y = filtfilt(b, a, X);
FTS = [X;Y];
Fv=(0:Nf/2-1)*Fs/Nf;
df=Fv(2)-Fv(1);
Iv = 1:numel(Fv);
y=ifft(Y,'symmetric');
Xplot=abs(X(1:Nf/2))/Ts; %graphs
Xplot=Xplot/max(Xplot);
Yplot=abs(Y(1:Nf/2))/Ts;
Yplot=Yplot/max(Yplot);

Thank you again.

[Star Strider]

As always, my pleasure.

I would simply plot the original unfiltered and filtered signal to look at both in the time domain, rather than taking the inverse Fourier transform of the Fourier transform of the filtered signal.

I do not see the second instance of ‘noise’ in the original filtered signal, so my conclusion is that what you are seeing is an artefact of taking the inverse Fourier transform of the filtered signal.

My complete code for your order 2 filter and a passband of [0.01 0.80] Hz is:

filename = 'acc time series.txt';
fidi = fopen(filename);
D = textscan(fidi, '%s%s');
t = str2double(strrep(D{1}, ',','.'));                          % Replace Decimal Separator
X = str2double(strrep(D{2}, ',','.'));
L = numel(t);

Ts = mean(diff(t));
Fs = 1/Ts;
Fn = Fs/2;
fhc = 0.8/Fn;
flc = 0.01/Fn;
Wn=[flc fhc];
n = 2;
[b,a] = butter(n,Wn,'bandpass');

Y = filtfilt(b, a, X);

S = [X Y];
FTS = fft(S)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
subplot(2,1,1)
plot(Fv, abs(FTS(Iv,1))*2)
grid
set(gca, 'XScale','log')
title('Unfiltered signal')
subplot(2,1,2)
plot(Fv, abs(FTS(Iv,2))*2)
grid
set(gca, 'XScale','log')
title('Filtered signal')

figure
freqz(b, a, 2^14, Fs)
set(subplot(2,1,1), 'XScale','log')
set(subplot(2,1,2), 'XScale','log')
title(subplot(2,1,1),'Filter Bode Plot')

figure
subplot(2,1,1)
plot(t, X)
grid
yl = ylim;
title('Original Signal')
subplot(2,1,2)
plot(t, Y)
grid
ylim(yl)
title('Filtered Signal')

Also, the sampling frequency of your signal is 100 Hz. You also initially said that you wanted the upper passband frequency to be 40 Hz, although you are using 0.8 Hz in your code. I get an entirely different result with an upper passband frequency of 0.8 Hz and with an order 2 filter:

The purpose of the Fourier transform is to see the effect of the filter on your signal. Once you have done that, and seen the result, nothing further is required of the Fourier-transformed signal.

[Eliza M.]

Yes. Not using the inverse Fourier transform seems to be the solution. Thank you for your help.

[Star Strider]

As always, my pleasure.

[Eliza M.]

Hi Star Strider! I was wondering if you have any experience with the Stockwell Transform ( for example the st.m file from : https://www.mathworks.com/matlabcentral/fileexchange/51808-time-frequency-distribution-of-a-signal-using-s-transform-stockwell-transform?focused=5484487&tab=function). If you can help me, I posted a recent question here : https://www.mathworks.com/matlabcentral/answers/418942-help-with-stockwell-transform-s-transform. Thank you in advance.

[Star Strider]

Hi Eliza!

I am not familiar with it. I had to look it up: S transform (link).

There also seems to be another implementation at CodeForge: st.m (link) that appears to have been written by Stockwell himself, posted in 2010.

[Eliza M.]

Yes this is the one I am using. Thank you anyway!

[Star Strider]

As always, my pleasure!