python处理fft

import numpy as np 
from scipy.fftpack import fft, fftshift, ifft
from scipy.fftpack import fftfreq
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import os

def read_txt(f):
    with open(f) as f:
#         print(f)
        data = f.read().strip().splitlines()
        data = [d[1:-1] for d in data]
        data = [d.split(',')[1:] for d in data]
        data = [[d[:4], d[4:]] for d in data]
        data = np.array(data).astype(np.float64)
        data = data.reshape(data.shape[0]*data.shape[1], data.shape[2])
        return data  
    
def do_fft(data, fs):
    Y = fft(data)
    # the positive part of fft, get from fft
    pos_Y_from_fft = Y[:Y.size//2]
    pos_Y_from_fft = np.abs(pos_Y_from_fft)
    num = len(pos_Y_from_fft)
    x = np.linspace(0,fs/2, num)
    return x,pos_Y_from_fft

两种切片操作 []

（1）给定范围型

a = np.array([10,11,12,13,14,15])
print(a[:3])   # 10,11,12
print(a[2:-1]) # 12, 13, 14

 （2）给定array型

给定array又分为两种

(A) 下标array

a = np.array([10,11,12,13,14,15])
print(a[np.array([0,3,5])]) # 10,13,15
print(a[np.array([2,4])]) # 12, 14 

 (B) True, False array

a = np.array([10,11,12,13,14,15])
print(a[np.array([True, True, False, False, False, False])]) # 10, 11
print(a[np.array([False, False, True, True, True, True])]) # 12, 13, 14, 15

 根据切片操作，可以做一个相对复杂点的例子：

a = np.array([2.1, 0.31,0.12,0.13,0.08,0.15, 0.16,0.12, 0.10,0.03, 0.08, 0.07])
b = np.array([0,1,2,3,4,5,6,7,8,9,10,11])
c = np.sum(a[np.where((b>7) & (b<13))])
print(c)  # 0.28 （=a[8]+a[9]+a[10]+a[11]））

 这里用到了np.where() 操作：

a = np.array([10,11,12,13,14,15])

np.where(a>13) # [False, False, False, False, True, True]

再来一个切片索引的示例：

FREQ_BANDS = {"delta": [0.5 , 4.5],
              "theta": [4.5 , 8.5],
              "alpha": [8.5 , 11.5],
              "sigma": [11.5, 15.5],
              "beta" : [15.5, 30],
              "gamma": [30  , 100] }

Fs = 200 # sampling rate


def psd_c1(array, Fs=200): # return delta,theta, alpha, sigma, beta, gamma, R for every channel.
    
    rows = len(array)
    freq = fftfreq(rows, 1/Fs)
    
#     fft_value = np.square(np.abs(fft(array)[:int(rows/2+0.5)]))
    fft_value = np.abs(fft(data)[:int(rows/2+0.5)])
    fft_norm = fft_value/(np.sum(fft_value))
    
    # delta
    cond_delta = np.where((freq>=FREQ_BANDS['delta'][0]) &
                          (freq<=FREQ_BANDS['delta'][1]))
    delta = np.sum(fft_norm[cond_delta])
    # theta
    cond_theta = np.where((freq>=FREQ_BANDS['theta'][0]) &
                          (freq<=FREQ_BANDS['theta'][1]))
    theta = np.sum(fft_norm[cond_theta])
    # alpha
    cond_alpha = np.where((freq>=FREQ_BANDS['alpha'][0]) &
                          (freq<=FREQ_BANDS['alpha'][1]))
    alpha = np.sum(fft_norm[cond_alpha])
    # sigma
    cond_sigma = np.where((freq>=FREQ_BANDS['sigma'][0]) &
                          (freq<=FREQ_BANDS['sigma'][1]))
    sigma = np.sum(fft_norm[cond_sigma])
    # beta
    cond_beta = np.where((freq>=FREQ_BANDS['beta'][0]) &
                          (freq<=FREQ_BANDS['beta'][1]))
    beta = np.sum(fft_norm[cond_beta])
    # gamma
    cond_gamma = np.where((freq>=FREQ_BANDS['gamma'][0]) &
                          (freq<=FREQ_BANDS['gamma'][1]))
    gamma = np.sum(fft_norm[cond_gamma])
    
    # R
    R = (alpha + sigma) / beta
    
    return np.array([delta, theta, alpha, sigma, beta, gamma, R])

python画图

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(21, 9))

ax1 = fig.add_subplot(511) 
ax2 = fig.add_subplot(512)  
ax3 = fig.add_subplot(513) 
ax4 = fig.add_subplot(514) 
ax5 = fig.add_subplot(515) 

ax1.plot(data[:,0])   # ax1.imshow() for image.
ax2.plot(data[:,1], color='b')
ax3.plot(data[:,2], color='r')
ax4.plot(data[:,3], color='g')
ax5.plot(data[:,3], color='c')


plt.show()

滤波

from scipy.fftpack import fft, fftfreq
from scipy import signal
from scipy.ndimage import gaussian_filter
from scipy.ndimage.filters import gaussian_filter1d

def preprocess(array):
    # butter filter
    b, a = signal.butter(8, [0.05, 0.5], 'bandpass') 
    data = signal.filtfilt(b, a, array, axis=0)
    # gaussian filter
    data = gaussian_filter1d(data, sigma=1.0, axis=0)
    return data

列表表达式

array = [np.mean(ch1[i-5:i+5]) for i in range(5,len(ch1)-5,10)]

peak检测

import numpy as np
from scipy.signal import find_peaks


def txt2array_from_mobile(filename):
    with open(filename) as f:
        data = f.read().strip().splitlines()
        data = [d.strip() for d in data] #去掉行首行尾空格
#         data = [d[1:-1] for d in data]
        data = [d.split(',')[2:] for d in data]
        return np.array(data, dtype=np.float64)



filename = './data/10.txt'
array = txt2array_from_mobile(filename)
# x = electrocardiogram()[2000:4000]
x1 = array[:,0]
x2 = array[:,1]
x3 = array[:,2]
x4 = array[:,3]
x1 -= np.mean(x1)
x2 -= np.mean(x2)
x3 -= np.mean(x3)
x4 -= np.mean(x4)
# peaks, _ = find_peaks(x, height=0, distance=5, prominence=(None, 0.6))
peaks1, r1 = find_peaks(x1, height=0.05, distance=10, width=2.5)
peaks2, r2 = find_peaks(x2, height=0.05, distance=10, width=2.5)
peaks3, r3 = find_peaks(x3, height=0.05, distance=10, width=2.5)
peaks4, r4 = find_peaks(x4, height=0.05, distance=10, width=2.5)


fig = plt.figure(figsize=(15,16))

ax1 = fig.add_subplot(4,1,1)
ax2 = fig.add_subplot(4,1,2)
ax3 = fig.add_subplot(4,1,3)
ax4 = fig.add_subplot(4,1,4)

ax1.plot(x1)
ax1.plot(peaks1, x1[peaks1], "x")
ax1.plot(np.zeros_like(x1), "--", color="gray")

ax2.plot(x2)
ax2.plot(peaks2, x2[peaks2], "x")
ax2.plot(np.zeros_like(x2), "--", color="gray")

ax3.plot(x3)
ax3.plot(peaks3, x3[peaks3], "x")
ax3.plot(np.zeros_like(x3), "--", color="gray")


ax4.plot(x4)
ax4.plot(peaks4, x4[peaks4], "x")
ax4.plot(np.zeros_like(x4), "--", color="gray")

fig.show()

peak检测2

def blink_detect(array):
    b, a = signal.butter(8, [0.025, 0.35], 'bandpass') 
    array = signal.filtfilt(b, a, array, axis=0)
    array = array - np.median(array, axis=0, keepdims=True)
    x0 = np.mean(array, axis=1)
    x1 = -x0
    peaks0, r0 = find_peaks(x0, height=[0.02,0.5], prominence=0.02, distance=20, width=4)
    peaks1, r1 = find_peaks(x1, height=[0.02,0.5], prominence=0.02, distance=20, width=4)
    peaks0 = peaks0[(r0['prominences']*r0['widths']>1.1)]
    peaks1 = peaks1[(r1['prominences']*r1['widths']>1.1)]
    return len(peaks0)+len(peaks1)