In [45]:
import pandas as pd
file_path = "./Physionet_ECGs/sel_100_tabsep.csv"
data = pd.read_csv(file_path, sep = '\t', index_col= 0, header=None)
data.head()
Out[45]:
Select the column you want to smooth
In [110]:
import numpy as np
import matplotlib.pyplot as plt
col = 2 #choose 1 or 2
values = data[col].values
times = data[col].index
Out[110]:
Fast Fourier Transform and plot the data
In [95]:
from numpy import fft
signal = np.array(values)
transformed = fft.fft(signal)
real = transformed.real
plt.plot(range(len(real)),real)
plt.ylim(1000,-1000)
plt.show()
In [117]:
ps = np.abs(transformed)**2
timestep = times[1]-times[0]
freqs = np.fft.fftfreq(signal.size,timestep)
idx = np.argsort(freqs)
print ps
plt.plot(freqs[idx], ps[idx])
#plt.plot(freqs, real)
plt.show()
In [135]:
#from glowing-python: http://glowingpython.blogspot.com/2011/08/how-to-plot-frequency-spectrum-with.html
from numpy import sin, linspace, pi
from pylab import plot, show, title, xlabel, ylabel, subplot
from scipy import fft, arange
def plotSpectrum(y,Fs=None):
"""
Plots a Single-Sided Amplitude Spectrum of y(t)
calcs the sampling rate from the first two elements of the array.
"""
Fs = Fs or 1/(signal[1]-signal[0]) #sampling rate
n = int(len(y)) # length of the signal
print n
k = arange(n)
T = n/Fs
frq = k/T # two sides frequency range
frq = frq[range(int(n/2))] # one side frequency range
Y = fft(y)/n # fft computing and normalization
Y = Y[range(int(n/2))]
plot(frq,abs(Y),'r') # plotting the spectrum
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
#subplot(2,1,1)
#plot(times,signal)
#xlabel('Time')
#ylabel('Amplitude')
plotSpectrum(signal[0:1000],Fs)
plt.ylim(0,.03)
show()
Mask the noise in the FFT'ed data and plot the signal with the mask
In [58]:
def mask_fft(transformed, masks=[(54000,171000)], ylims=(1000,-1000)):
'''
The default settings work well for column 1
'''
masked = transformed
for i in masks:
mask_min = i[0]
mask_max = i[1]
masked[mask_min:mask_max] = 0 + 0j
real2 = transformed.real
plt.plot(range(len(real2)),real2)
plt.ylim(*ylims)
return masked
masked = mask_fft(transformed, [(40000,190000)])
plt.gca()
plt.show()
reverse FFT the signal and plot the smoothed data
In [55]:
ifft = fft.ifft(transformed)
plt.plot(times, ifft.real, "r", label = 'transformed')
plt.plot(times, signal, 'b', label="original")
plt.legend()
plt.xlim(0,1000)
plt.show()
Save!
In [ ]:
In [6]:
%load_ext rmagic
In [70]:
%%R -i values -o values_spect
values_spect = spectrum(values)
#print(summary(values_spect))
In [94]:
print values_spect.names
for i in values_spect.items():
print i
print
In [38]:
values = ifft
In [39]:
%%R -i values -o sg_smoothed
sg5 = c(-3, 12, 17, 12, -3)/35
sg_smoothed = filter(values,sg5)
In [40]:
df_sg = pd.DataFrame(data = np.array(sg_smoothed), index = times, columns=[col]).dropna()
plt.plot(df_sg.index, df_sg.values, lw=2)
plt.plot(times, values, 'r', lw=1)
#plt.plot(times,ifft.real,"g", lw=1)
plt.show()
In [41]:
%%R -i values -o smoothed
ma5 = c(1, 1, 1, 1, 1)/5;
smoothed =filter(values,ma5)
print(smoothed[0:100])
In [44]:
df = pd.DataFrame(data = np.array(smoothed), index = times, columns=[col])
df_clean = df.dropna()
plt.plot(df_clean.index, df_clean.values, lw=3, label="ifft MA Smoothed")
plt.plot(times, values, 'r', lw=2, label="original")
plt.plot(times,ifft.real,"g", lw=1, label="ifft")
plt.legend()
plt.show()
In [9]:
len(df_clean.values), len(df_clean.index)
Out[9]:
In [10]:
df_clean.to_csv('./Physionet_ECGs/sel_100_tabsep_mvsmooth.csv', header=False, sep="\t")
In [ ]:
def calc_slopes(X, Y):
dY = (np.roll(Y, -1, axis=1) - Y)[:-1]
dX = (np.roll(X, -1, axis=0) - X)[:-1]
slopes = dY/dX
return slopes
calc_slopes(df_clean.index, df_clean.values)
In [62]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import RadioButtons
t = np.arange(0.0, 2.0, 0.01)
s0 = np.sin(2*np.pi*t)
s1 = np.sin(4*np.pi*t)
s2 = np.sin(8*np.pi*t)
fig, ax = plt.subplots()
l, = ax.plot(t, s0, lw=2, color='red')
plt.subplots_adjust(left=0.3)
axcolor = 'lightgoldenrodyellow'
rax = plt.axes([0.05, 0.7, 0.15, 0.15], axisbg=axcolor)
radio = RadioButtons(rax, ('2 Hz', '4 Hz', '8 Hz'))
def hzfunc(label):
hzdict = {'2 Hz':s0, '4 Hz':s1, '8 Hz':s2}
ydata = hzdict[label]
l.set_ydata(ydata)
plt.draw()
radio.on_clicked(hzfunc)
rax = plt.axes([0.05, 0.4, 0.15, 0.15], axisbg=axcolor)
radio2 = RadioButtons(rax, ('red', 'blue', 'green'))
def colorfunc(label):
l.set_color(label)
plt.draw()
radio2.on_clicked(colorfunc)
rax = plt.axes([0.05, 0.1, 0.15, 0.15], axisbg=axcolor)
radio3 = RadioButtons(rax, ('-', '--', '-.', 'steps', ':'))
def stylefunc(label):
l.set_linestyle(label)
plt.draw()
radio3.on_clicked(stylefunc)
plt.show()
In [ ]: