How do I perform zero phase filtering in Matlab?

Hello All,

I am a student working on a tremor suppression project for those with parkinsons and I would really appreciate your expertise in Matlab in performing FFT's and applying filters.

What we need to do is filter out frequencies with a butterworth bandpass from 3-9 Hz and apply a phase shift to that data that will be sent to our neuro-stimulator.

I have code that properly filters our data between 3-9 Hz and I can use the FFT to verify that our filter is working, but I am having trouble getting usable phase information from our post-filtered data. I have used the filtfilt function in matlab, but for some reason our phase data is still not perfect to the original. I understand the nature of filters to change phase but we need some way to do this so that we can apply a phase shift to the data later.

I have attached my code here:

fs = 100;
% fs is Sampling frequency in hz
% Time vector - use colon operator to generate integer vector of sample
% numbers
t = 1:.01:7.01; %7 seconds of sampling
f1 = 6; %6 Hz freq
f2 = 40; %40 Hz freq
f3 = 45; %45 Hz freq
signal = 3*cos(2*pi*f1*t + 2) + 1*cos(2*pi*f2*t + 1) + 1*cos(2*pi*f3*t + 3);

fs=100; % sample frequency
fc1=2; %cutoff frequency low
fc2=15; %cutoff frequency high
[B,A]=butter(3,[fc1/(fs/2) fc2/(fs/2)] ,'bandpass'); 
filtered_signal=filtfilt(B,A,signal);

plot(signal)
hold on
plot(filtered_signal, 'r')
figure

X = fft(signal);
Z = fft(filtered_signal);

n=length(X);
X_mag = abs(X);
X_phase = angle(X);
w = fs/(length(t)-1);
binHzConv = ((fs)/(n) +.000001);
freq = 0:binHzConv:100;

p = unwrap(angle(X)); 
g = unwrap(angle(Z));

Any help at all would be much appreciated! Thank you.

Answer 1

[Daniel kiracofe]

First, your line "freq = 0:binHzConv:100;" is not correct. Frequency does not range between 0 and 100 Hz in your example, but rather -50 to 50. This is the so-called "two sided" spectrum. I have some more details on that here: http://www.mechanicalvibration.com/Making_matlab_s_fft_functio.html

The second thing, is that I suggest you plot both amplitude and phase together. I.e. using the "simpleFFT" script in the above web page:

[frq, amp1, phase1] = simpleFFT( signal, fs);
[frq, amp2, phase2] = simpleFFT( filtered_signal, fs);
figure;
subplot(2,1,1);
plot( frq, amp1, frq, amp2);
xlim([0 10])
subplot(2,1,2);
 plot( frq, phase1, frq, phase2);
 xlim([0 10])

You will see that your filtered signal only has a large amplitude near 6 Hz. For example, at 2 Hz, it's about 100 times smaller than at 6 Hz. This is of course exactly what you had intended by use of the filter. But the thing to keep in mind that the phase of a signal with zero amplitude is basically meaningless. It's like asking whether zero is a "positive zero" or a "negative zero". The phase at 6 Hz in your filtered signal does match exactly the phase at 6 Hz in your original signal. At other frequencies it does not match, just due to numerical roundoff, but you don't really care because the amplitude is zero and phase is meaningless there anyway.

It may help you to understand if you plot the real() and imag() components of X and Z instead of amplitude and phase.

[Chiraag Lathia]

Thank you very much for your help! I think I have a better understanding of the fft in matlab now.

Best regards.