Quy - 8 months ago 42
Python Question

# Derived Result from Butter Filter and FFT doesn't change over time

I am very new in Signal Processing, have met a situation that I am not sure if it is correct or not. Please correct me then I will update more details.

My data is here

I acquired an accelerometer signal taken from my cellphone (Samsung Galaxy Note 2, sampling rate $\approx 99 Hz$). I would like to analyze frequency from $0.3 Hz$ to $5.0 Hz$

My procedure is following steps:

1. combination: let say a sensor yields 3 channels $x$, $y$, $z$. The combination is to produce a new channel $v = \sqrt{(x * x + y * y + z * z)}$

2. Perform a median filter: to make signal smoothly

3. Butter-worth filter: my cutoff is from $0.3 Hz$ to $5.0 Hz$

4. FFT

Below image is my demonstrate with segment of 120 time-points with 4 steps: (can explore more at my video)

I observed result of step 3 and 4 do not change while signal varies over time

My question is if there is anything I can make sure this result is correct or not? Thanks in advance

Below is my code was used for applying filters

from __future__ import division
import numpy as np
from numpy.fft import rfft, rfftfreq
from numpy import absolute
import matplotlib.pyplot as plt
from scipy.signal import medfilt, hilbert
import pandas as pd

chunk = 120
LOW_CUT = 0.3
HIGH_CUT = 5.0
FS = 99
freqs = rfftfreq(chunk, 1 / FS)
for k, g in _accel.groupby(np.arange(len(_accel)) // chunk):
_v = g['v'].values
_v = medfilt(_v, 7)
_v = butter_bandpass_filter(_v, LOW_CUT, HIGH_CUT, FS, order=4)
v = 1 / chunk * absolute(rfft(_v))
plt.stem(freqs, v)

update 2 updated sampling rate in code
FS = 99

update 3 increased chunk size to 512, plotted data again. Made a video of result without bandpass

I give the problem a quick try, and below is my snippet

data = _accel['v'].tolist()
Fs = 99
# remove the DC part, to help the plotting later
data = data - np.mean(data)
# Perform FFT for real data, on the whole 6000 samples,
# using 4096 discrete frequencies, which is dense enough to capture
# the frequency information within 0.3-5 Hz.
fdata = rfft(data,4096)

# the frequencies we are looking at in the FFT
freqs = map(lambda x: float(x)*Fs/4096, range(0,4097))

# Plot
plt.plot(freqs[0:2049],fdata)
plt.xlabel('Frequency')
plt.show()

The resulting plot does contain information in the band you are interested in. Plot of frequency magnitude

I guess your problem is in choosing chunk too small. The resolution in the frequency domain is Fs/N, with N is the number of points to perform FFT (and usually the length of the signal vector in the time domain). So, if you want to capture information in the range of 0.3-5Hz, I assume you would need a resolution of about 0.2Hz, which means N should be at least 500. Your choice of 120 for window length is obviously not enough.