# Source code for signals

'''
Utilities for signal processing
'''

'''
Author: Thomas Haslwanter
Version: 1.2
Date: Nov-2013
'''

import numpy as np
import matplotlib.pyplot as plt

import math
from numpy import dot

[docs]def pSpect(data, rate):
'''
Power spectrum and frequency

Parameters
----------
data : array, shape (N,)
measurement data
rate : float
sampling rate [Hz]

Returns
-------
powerspectrum : array, shape (N,)
frequency : array, shape (N,)

Example
-------
>>> pxx, freq = pSpect(data, 1000)

'''

from scipy.fftpack import fft
nData = len(data)
window = np.hamming(nData)
fftData = fft(data*window)
PowerSpect = fftData * fftData.conj() / nData
freq = np.arange(nData) * float(rate) / nData
return (PowerSpect, freq)

[docs]def show_se(raw):
'''Show mean and standard error, of a dataset in column form.

Parameters
----------
raw : array (N,M)
input data, M sets of N data points

Returns
-------
avg : array (N,)
average value
se : array (N,)
standard error

Examples
--------
>>> t = np.arange(0,20,0.1)
>>> x = np.sin(t)
>>> data = []
>>> for ii in range(10):
>>>     data.append(x + np.random.randn(len(t)))
>>> show_se(np.array(data).T)

Notes
-----
.. image:: _static/show_se.png
:scale: 50%

'''

N = len(raw)

# Calculate mean and standard error
avg = np.mean(raw, axis=1)
std = np.std(raw, axis=1, ddof=1)
se = std/np.sqrt(N)

# Calculate upper and lower limit, for showing the standard error
upper = avg + se
lower = avg - se

# Plot the data
plt.fill_between(t, lower, upper, color='gray', alpha=0.5)
plt.hold(True)
plt.plot(t,avg)
plt.show()

return (avg, se)

[docs]def corrvis(x,y):
'''
Visualize correlation, by calculating the cross-correlation of two signals.
The aligned signals and the resulting cross correlation value are shown,
and advanced when the user hits a key or clicks with the mouse.

Parameters
----------
X : array (N,)
Comparison signal

Y : array (M,)
Reference signal

Examples
--------
>>> x = np.r_[0:2*pi:10j]
>>> y = np.sin(x)
>>> corrvis(y,y)

Notes
-----
Based on an idea from dpwe@ee.columbia.edu

'''

Nx = x.size
Ny = y.size
Nr = Nx + Ny -1

xmin = -(Nx - 1)
xmax = Ny + Nx -1

# First plot: Signal 1
ax1 = plt.subplot(311)
ax1.plot(range(Ny), y)
ax = ax1.axis()
ax1.axis([xmin, xmax, ax[2], ax[3]])
ax1.grid(True)
ax1.set_xticklabels(())
ax1.set_ylabel('Y[n]')

# Precalculate limits of correlation output
axr = [xmin, xmax, np.correlate(x,y,'full').min(), np.correlate(x,y,'full').max()]

# Make a version of y padded to the full extent of X's we'll shift
padY = np.r_[np.zeros(Nx-1), y, np.zeros(Nx-1)]
R = []

# Generate the cross-correlation, step-by-step
for p in range(Nr):

# Figure aligned X
ax2 = plt.subplot(312)
ax2.hold(False)
ax2.plot(np.arange(Nx)-Nx+p+1, x)
ax = ax2.axis()
ax2.axis([xmin, xmax, ax[2], ax[3]])
ax2.grid(True)
ax2.set_ylabel('X[n-l]')
ax2.set_xticklabels(())

# Calculate correlation
# Pad an X to the appropriate place

# Third plot: cross-correlation values
ax3 = plt.subplot(313)
ax3.hold(False)
ax3.plot(np.arange(len(R))-(Nx-1), R, linewidth=2)
ax3.axis(axr)
ax3.grid(True)
ax3.set_ylabel('Rxy[l]')

# Update the plot
plt.draw()
plt.waitforbuttonpress()

plt.show()

if __name__ == '__main__':
t = np.arange(0,10,0.1)
x = np.sin(t) + 0.2*np.random.randn(len(t))
smoothed = savgol(x, 11)
plt.plot(t, smoothed)
plt.show()
print('Done')