Low pass filter implementation Correct or wrong?

Question

I've been working on 2 sensors signals measuring vibrations of a rotating shaft. Since there is residual noise in the signal. I tried to filter it by detrending, zero padding and applying low pass filter. Below I'm attaching the graphs of the signal before and after filtering. There is a huge variation in the signal after filtering that makes me think if I'm really doing it in the right way.

enter image description here

My Matlab code is

X = xlsread(filename,'F:F');
Y = xlsread(filename,'G:G');

%Calculate frequency axis
fs = 1e6  ; % Sampling frequency (Hz)
NFFT = 2^nextpow2(length(X));  % Zero padding to nearest N power 2
df = fs/NFFT;
dt = 1/df;
%Frequency Axis defintion
f = (-(fs-df)/2:df:(fs-df)/2)';


X(2^ceil(log2(length(X))))=0;
Y(2^ceil(log2(length(Y))))=0;

%calculate time axis
T = (dt:dt:(length(X)*dt))';

subplot(2,2,1)
plot(T,X);
xlabel('Time(s)')
ylabel('X amplitude')
title('X signal before filtering')
subplot(2,2,2)
plot(T,Y);
xlabel('Time(s)')
ylabel('Y amplitude')
title('Y signal before filtering')

 X = detrend(X,0); % Removing DC Offset
 Y = detrend(Y,0); % Removing DC Offset


 % Filter parameters:
 M = length(X);    % signal length
 L = M;          % filter length   
 fc = 2*(38000/60);   % cutoff frequency


 % Design the filter using the window method:
 hsupp = (-(L-1)/2:(L-1)/2);
 hideal = (2*fc/fs)*sinc(2*fc*hsupp/fs);
 h = hamming(L)' .* hideal; % h is our filter



 % Zero pad the signal and impulse response:
 X(2^ceil(log2(M)))=0;
 xzp = X; 
 hzp = [ h zeros(1,NFFT-L) ];

 % Transform the signal and the filter:
 X = fft(xzp);
 H = fft(hzp)';
 X = X .* H;
 X = ifft(X);
 relrmserrX = norm(imag(X))/norm(X); % checked... this for zero
 X = real(X)';

 % Zero pad the signal and impulse response:
 Y(2^ceil(log2(M)))=0;
 xzp = Y; 
 hzp = [ h zeros(1,NFFT-L) ];

 % Transform the signal and the filter:
 Y = fft(xzp);
 H = fft(hzp)';
 Y = Y .* H;
 Y = ifft(Y);
 relrmserrY = norm(imag(Y))/norm(Y); % check... should be zero
 Y = real(Y)';

I plotted the after filtering and as you can see in the picture there is a clear deviation. I want to only filter noise but the signal seem to loose other components and I'm little confused if thats the right way of doing filtering. Any suggestion, hints or ideas would be helpful.

In the end I want to plot X vs Y to give orbits of the shaft vibration. Please also find below another picture of the unfiltered and filtered orbit. As you can see in the picture there is also change in the orbits from the original one (left image with lot of noise).

Unfiltered Orbit and Filtered Orbit

P.S.: I don't have DSP tool box

Thanks for specifying you don't have the DSP toolbox. Suggestions are much different as a result. — Juderb
Except for the phase and time direction being wrong (probably an artefact of the specific way matlab outputs fft data... it's a bit odd btw, check the docs), I can't see what's supposed to be wrong here. — Lucas
I haven't dug through your code, but the difference in orbits does show something isn't right, at least for your application. The ringing where the signal is cut off and the smoothness of the orbit plot look great - a good sign that it is definitely filtering out the high frequency components. Does this filter have a linear phase response? You'll want that for your application. A symmetric impulse response should ensure a linear phase. — Katie
After staring at the graphs a bit, I'm more convinced that its the matlab fft output format catching you out. — Lucas
Please check the way you are defining H. It has a time-shift with respect to how you defined X, which would result in a final time-shift. — Lazarus

Lazarus Lazarus · Accepted Answer · 2015-03-16T16:19:36

There is no problem with your FFT and IFFT.

You can test your code with a simple sine wave, with a very low frequency. The final output should be (almost) the same sine wave, since you are low-pass filtering it.

You've defined X (and Y) from 0 to some value. However, you've defined H from - (L-1)/2 to some positive value. Mathematically, this is fine, but you are simply taking an fft of H. Matlab thinks this has the same time-scale as X, when you multiply the ffts together!

So, in reality, you've taken the fft of XHfft(delta(t-d)), where d is the resulting time-shift. You can either undo this time-shift in the frequency domain, or the time-domain.

Low pass filter implementation Correct or wrong?

1 Answers