1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
|
import core.modules
import core.modules.module_registry
from core.modules.vistrails_module import Module, ModuleError
from Matrix import *
from Array import *
import scipy
import scipy.signal
from scipy import fftpack
import numpy
class DSPModule(object):
my_namespace = 'scipy|signals'
class SignalGenerator(DSPModule, Module):
my_namespace = 'scipy|signals|generator'
def compute(self):
samples = self.getInputFromPort("Samples")
periods = self.getInputFromPort("Periods")
freqs = self.getInputListFromPort("Frequencies")
ar = numpy.linspace(0., float(periods) * 2. * scipy.pi, periods * samples)
out_ar = numpy.zeros(periods * samples)
for f in freqs:
out_ar += scipy.sin(f * ar)
out = NDArray()
out.set_array(out_ar)
self.setResult("Output", out)
@classmethod
def register(cls, reg, basic):
reg.add_module(cls, namespace=cls.my_namespace)
reg.add_input_port(cls, "Samples", (basic.Integer, "Sampling Rate"))
reg.add_input_port(cls, "Periods", (basic.Integer, "Signal Length"))
reg.add_input_port(cls, "Frequencies", (basic.Float, "Additive Frequency"))
reg.add_output_port(cls, "Output", (NDArray, "Output Signal"))
class FFT(DSPModule, Module):
__doc__ = """ Calculate the discrete Fourier transform of the arbitrary
sequence presented on the Signal port. This is done using
SciPy's FFTPack module.\n\n"""
__doc__ += """From fftpack.fft:\n\t"""
__doc__ += fftpack.fft.__doc__
my_namespace = 'scipy|signals|fourier'
def compute(self):
sig_array = self.getInputFromPort("Signals")
# If there is no input on the samples port,
# use the number of samples in an array row for
# the number of fft points.
if self.hasInputFromPort("Samples"):
pts = self.getInputFromPort("Samples")
else:
try:
pts = sig_array.get_shape()[1]
except:
pts = sig_array.get_shape()[0]
sh = sig_array.get_shape()
if len(sh) < 2:
shp = (1, sh[0])
sig_array.reshape(shp)
(num_sigs, num_samps) = sig_array.get_shape()
phasors = fftpack.fft(sig_array.get_row_range(0,0), pts)
out_ar = phasors
for i in xrange(1,num_sigs):
phasors = fftpack.fft(sig_array.get_row_range(i,i), pts)
out_ar = numpy.vstack([out_ar, phasors])
out = NDArray()
out.set_array(out_ar)
self.setResult("FFT Output", out)
@classmethod
def register(cls, reg, basic):
reg.add_module(cls, namespace=cls.my_namespace)
reg.add_input_port(cls, "Signals", (NDArray, 'Input Signal Array'))
reg.add_input_port(cls, "Samples", (basic.Integer, 'FFT Samples'))
reg.add_output_port(cls, "FFT Output", (NDArray, 'FFT Output'))
class FFTN(DSPModule, Module):
__doc__ = """ Calculate the discrete Fourier transform of the arbitrary
sequence presented on the Signal port. This is done using
SciPy's FFTPack module.\n\n"""
__doc__ += """From fftpack.fftn:\n\t"""
__doc__ += fftpack.fftn.__doc__
my_namespace = 'scipy|signals|fourier'
def compute(self):
sig_array = self.getInputFromPort("Signals")
# If there is no input on the samples port,
# use the number of samples in an array row for
# the number of fft points.
if self.hasInputFromPort("Samples"):
pts = self.getInputFromPort("Samples")
else:
pts = sig_array.get_shape()[1]
sh = (sig_array.get_shape()[0], pts)
phasors = fftpack.fftn(sig_array.get_array(), shape=sh)
out = NDArray()
out.set_array(phasors)
self.setResult("FFT Output", out)
@classmethod
def register(cls, reg, basic):
reg.add_module(cls, namespace=cls.my_namespace)
reg.add_input_port(cls, "Signals", (NDArray, 'Input Signal Array'))
reg.add_input_port(cls, "Samples", (basic.Integer, 'FFT Samples'))
reg.add_output_port(cls, "FFT Output", (NDArray, 'FFT Output'))
class ShortTimeFourierTransform(DSPModule, Module):
""" Calculate the short time Fourier transform of the
sequence presented on the Signal port. This is done using
SciPy's FFTPack fft module in conjuction with an input window.
If a window is not specified, a Hamming window of the specified
size is used. """
my_namespace = 'scipy|signals|fourier'
def get_signal(self, sigs, window, offset, size):
win = scipy.zeros(sigs.shape[0])
win[offset:offset+size] = window
part = sigs * win
return part
def compute(self):
sigs = self.getInputFromPort("Signals")
sr = self.getInputFromPort("SamplingRate")
out_vol = None
if self.hasInputFromPort("Window"):
window = self.getInputFromPort("Window").get_array()
win_size = window.shape[0]
else:
win_size = self.getInputFromPort("Window Size")
window = scipy.signal.hamming(win_size)
if self.hasInputFromPort("Stride"):
stride = self.getInputFromPort("Stride")
else:
stride = int(win_size / 2)
sh = sigs.get_shape()
if len(sh) < 2:
shp = (1, sh[0])
sigs.reshape(shp)
(num_sigs, num_samps) = sigs.get_shape()
for i in xrange(num_sigs):
offset = 0
signal = sigs.get_array()[i]
# We need to do the first window here so that we
# can have something to call vstack on.
sig = self.get_signal(signal, window, offset, win_size)
im_array = fftpack.fft(sig)
offset += stride
while 1:
try:
sig = self.get_signal(signal, window, offset, win_size)
phasors = fftpack.fft(sig)
offset += stride
im_array = numpy.vstack([im_array, phasors.ravel()])
except:
break
# STFT of one signal is done. Clean up the output
(slices, freqs) = im_array.shape
ar = im_array[0:,0:sr*2]
ar = ar[0:,::-1]
if out_vol == None:
out_vol = ar
ovshape = out_vol.shape
out_vol.shape = 1, ovshape[0], ovshape[1]
else:
arshape = ar.shape
ar.shape = 1, arshape[0], arshape[1]
out_vol = numpy.vstack([out_vol, ar])
# All signals have been processed and are in the volume.
out = NDArray()
out.set_array(out_vol)
self.setResult("FFT Output", out)
@classmethod
def register(cls, reg, basic):
reg.add_module(cls, namespace=cls.my_namespace)
reg.add_input_port(cls, "Signals", (NDArray, 'Signal Array'))
reg.add_input_port(cls, "SamplingRate", (basic.Integer, 'Sampling Rate'))
reg.add_input_port(cls, "Window", (NDArray, 'Windowing Function'))
reg.add_input_port(cls, "Window Size", (basic.Integer, 'Window Size'))
reg.add_input_port(cls, "Stride", (basic.Integer, 'Stride'))
reg.add_output_port(cls, "FFT Output", (NDArray, 'FFT Output'))
class SignalSmoothing(DSPModule, Module):
"""
Documentation
"""
def compute(self):
window = self.getInputFromPort("Window").get_array()
in_signal = self.getInputFromPort("Signal").get_array()
to_conv = window/window.sum() # Make sure the window is normalized
if in_signal.ndim > 1:
out_ar = numpy.zeros(in_signal.shape)
else:
out_ar = numpy.zeros(1,in_signal.shape[0])
in_signal.shape = (1, in_signal.shape[0])
for row in xrange(in_signal.shape[0]):
out_ar[row] = numpy.convolve(to_conv, in_signal[row], mode='same')
out = NDArray()
out.set_array(out_ar)
self.setResult("Output Array", out)
@classmethod
def register(cls, reg, basic):
reg.add_module(cls, namespace=cls.my_namespace)
reg.add_input_port(cls, "Signal", (NDArray, "Input Signals"))
reg.add_input_port(cls, "Window", (NDArray, "Smoothing Filter"))
reg.add_output_port(cls, "Output Array", (NDArray, 'Smoothed Signals'))
# class SingleTrialPhaseLocking(DSPModule, Module):
# """
# Documentation
# """
# def get_time_indexes(self, t0, time_window):
# if time_window % 2:
# # odd number of samples: t0 +/- (window-1)/2
# tw = (time_window - 1) / 2
# else:
# tw = time_window / 2
# return (t0 - tw, t0 + tw)
# def calc_pli(self, f_n_ar, f_m_ar):
# phasors = numpy.concatenate((f_n_ar, f_m_ar))
# norm_c = numpy.sqrt(phasors.real*phasors.real + phasors.imag*phasors.imag)
# phasors /= norm_c
# mean_phasor = phasors.mean()
# pli = numpy.sqrt(mean_phasor.real*mean_phasor.real + mean_phasor.imag*mean_phasor.imag)
# return pli
# def compute(self):
# phasors = self.getInputFromPort("Phasor Array").get_array()
# time_window = self.getInputFromPort("Time Window")
# time_step = self.forceGetInputFromPort("Time Step")
# if time_step == None:
# time_step = 1
# ndims = phasors.ndim
# if ndims == 2:
# phasors.shape = (1, phasors.shape[0], phasors.shape[1])
# elif ndims == 1:
# phasors.shape = (1, phasors.shape[0], 1)
# else:
# raise ModuleError("Cannot Process Phasor set of dimension " + str(ndims))
# num_freqs = phasors[0].shape[0]
# num_times = phasors[0].shape[1]
# num_times /= time_step
# out_ar = numpy.zeros((phasors.shape[0], num_freqs, num_freqs, num_times))
# for channel in xrange(phasors.shape[0]):
# tfr = phasors[channel,:,:].squeeze()
# for f_m in xrange(tfr.shape[0]):
# f_m_row = tfr[f_m,:]
# for f_n in xrange(f_m+1, tfr.shape[0], 1):
# f_n_row = tfr[f_n,:]
# t0 = 0
# tn = 0
# (start_i, end_i) = self.get_time_indexes(t0, time_window)
# while t0 < f_m_row.shape[0]:
# f_m_range = f_m_row[max(0,start_i):min(end_i,f_m_row.shape[0]-1)]
# f_n_range = f_n_row[max(0,start_i):min(end_i,f_n_row.shape[0]-1)]
# pli = self.calc_pli(f_m_range, f_n_range)
# out_ar[channel, f_m, f_n, tn] = pli
# out_ar[channel, f_n, f_m, tn] = pli
# out_ar[channel, f_m, f_m, tn] = 1.0
# tn += 1
# t0 += time_step
# start_i += time_step
# end_i += time_step
# out = NDArray()
# out.set_array(out_ar)
# self.setResult("Output Array", out)
# @classmethod
# def register(cls, reg, basic):
# reg.add_module(cls, namespace=cls.my_namespace)
# reg.add_input_port(cls, "Phasor Array", (NDArray, 'Phasor Array'))
# reg.add_input_port(cls, "Time Step", (basic.Integer, 'Stride in the Time Domain'))
# reg.add_input_port(cls, "Time Window", (basic.Integer, 'Samples per Timeslice'))
# reg.add_output_port(cls, "Output Array", (NDArray, 'Result set'))
# class CalculatePhaseLocking(DSPModule, Module):
# """
# documentation
# """
# def Phi(self, p):
# return scipy.arctan2(p.real, p.imag)
# def compute(self):
# phasors = self.getInputFromPort("Phasor Array").get_array()
# if phasors.ndim != 3:
# raise ModuleError("Cannot handle phasor array with less than 3 dimensions")
# (trials, times, frequencies) = phasors.shape
# lowestF = self.getInputFromPort("Lowest Freq")
# # highestF = self.getInputFromPort("Highest Freq")
# Phi = self.Phi(phasors)
# gamma_ar = numpy.zeros((times, frequencies, frequencies))
# for t in range(times):
# n = lowestF
# for fn_i in range(frequencies):
# n += fn_i
# fn = fn_i + lowestF
# m = lowestF
# for fm_i in range(fn_i, frequencies, 1):
# m += fm_i
# fm = fm_i + lowestF
# DeltaPhi = (float((n+m)/(2*m))Phi[:,t,fm_i] - float((n+m)/(2*n))Phi[:,t,fn_i]) % (2. * scipy.pi)
# Gamma = exp(complex(0.,DeltaPhi))
# Gamma = Gamma.sum()
# Gamma = numpy.sqrt(Gamma * Gamma.conjugate())
# gamma_ar[t, fn_i, fm_i] = Gamma
# gamma_ar[t, fm_i, fn_i] = Gamma
# m += 1
# n += 1
# out = NDArray()
# out.set_array(gamma_ar)
# self.setResults("Gamma", out)
# @classmethod
# def register(cls, reg, basic):
# reg.add_module(cls, namespace=cls.my_namespace)
# reg.add_input_port(cls, "Phasor Array", (NDArray, 'Phasor Array'))
# reg.add_input_port(cls, "Lowest Freq", (basic.Integer, 'Lowest Frequency Phasor'))
# reg.add_output_port(cls, "Gamma", (NDArray, 'Phase Locking Volume'))
# class DifferentialPhaseLocking(DSPModule, Module):
# """
# documentation
# """
# def compute(self):
# phasors = self.getInputFromPort("Phasor Array").get_array()
# if phasors.ndim != 2:
# raise ModuleError("Cannot handle phasor array with more than 2 dimensions")
# mag = phasors.real*phasors.real + phasors.imag*phasors.imag
# mag = numpy.sqrt(mag)
# normalized = phasors / mag
|
