Quy - 1 year ago 73

Python Question

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:

- 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)}$
- Perform a median filter: to make signal smoothly
- Butter-worth filter: my cutoff is from $0.3 Hz$ to $5.0 Hz$
- 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)

_accel = pd.read_csv('data.csv')

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)

`FS = 99`

Answer Source

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.