Commit d632e645 authored by Federico Asara's avatar Federico Asara

initial commit with all the stuff, even the irrelevant ones

parent fd3e10b4
from numpy import *
import sinefit, ScipySineFit
class WindowedSignal(object):
pass
class FFTSignal(object):
"""a representation of a frequency-domain signal
"""
window_function = {
'No window' : [ones, []],
'BARTLETT' : [bartlett, []],
'BLACKMAN' : [blackman, []],
'HAMMING' : [hamming, []],
'HANNING' : [hanning, []],
'KAISER' : [kaiser, [1]],
}
WINDOW_TYPES = window_function.keys()
selectedWindow = 'No window'
data = {}
def __init__(self, origin, dB = True):
"""initialize an FFTSignal object
origin: its counterpart in time domain
dB: true when in dB units, false if raw amplitude
"""
# TODO make it possible to create a fd signal without its td counterpart
self.origin = origin
self.dB = dB
self.precalculateAll()
def get(self, what = None):
if what is None:
what = self.selectedWindow
# just to be sure
if what not in self.WINDOW_TYPES:
what = self.WINDOW_TYPES[0]
print "Get request: %s --> %s" % (what, self.data[what])
return self.data[what]
def precalculateAll(self):
self.data = {}
for i in self.WINDOW_TYPES:
print "Precalculating fft data for", i
self.data[i] = self.precalculate(i)
def precalculate(self, windowName):
# windows are ok right now
# wrong values: THD
# wrong formatting output: THD
output = WindowedSignal()
window = self.window_function[windowName][0]
# DFT of the signal scaled by window
windowParams = [self.origin.nsamples] + self.window_function[windowName][1]
output.dft = 10*log10(abs(fft.rfft(self.origin.data * window(*windowParams))))
# m = len(output.dft)
# sinusoid frequency detection
# we need a first approximation, let's search for the maximum in
# abs(output.dft) -- we need its index
print "Calculating frequency: rate = %d, nsamples = %d" % (self.origin.rate, self.origin.nsamples)
w0index = argmax(output.dft)
w0 = w0index * (2 * pi * self.origin.rate) / self.origin.nsamples
freqSample = 2*pi*(self.origin.rate)**-1
args = (array(self.origin.data), freqSample, w0)
print "index:", w0index, "array: %s\nFs [Hz]: %f\tw0 [rad/sec]:%f" % args
output.w0 = ScipySineFit.sinefit4(*args)#array(self.origin.data), freqSample, w0)
print "Frequency detected: %.6f" % output.w0
# let's go through harmonic peaks. we will produce a generator
# this way
def adjust(data, fs):
while abs(data -fs) > (fs/2.):
data += fs if data < fs else -fs
return data
def getGenerator(start, end):
return (adjust(x, freqSample) for x in xrange(start, end))
output.harmonicPeaksGenerator = getGenerator
# THD seems easy..
tenHarmonicsIndexes = getGenerator(0, 10)
tenHarmonics = array([output.dft[i] for i in tenHarmonicsIndexes])
output.THD = 10.0 * log10(sum(10 ** (tenHarmonics/2.0)))
# we need the avg of HDs
avgHarmonics = mean(tenHarmonics)
# we also need the noise floor: the average of all components below
# avgHarmonics
filteredNoise = where(output.dft >= avgHarmonics, output.dft, 0)
output.noiseFloor = mean(filteredNoise)
# now we can evaluate SNR = max component - noise floor - process gain
output.SNR = output.dft[w0index] -output.noiseFloor -self.origin.cachedProcessGain
# now it's time for SINAD
output.SINAD = -10 * log10(10**(-output.SNR/10) + 10**(-abs(output.THD)/10))
# missing ENOB, how can I get FSR & RSR?
output.ENOB = 0
# missing SFDR
output.SFDR = 0
return output
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 8, 2011 2:51:20 PM$"
from pylab import *
from scipy import *
import sinefit
from scipy import optimize
def sinefit4(data, Ts, w0):
# Target function
# fitfunc = lambda p, x: p[0]*cos(2*pi/p[1]*x+p[2]) + p[3]*x
fitfunc = lambda p, x: p[0]*cos(p[3]*x) +p[1]*sin(p[3]*x) +p[2]
# Distance to the target function
errfunc = lambda p, x, y: fitfunc(p, x) - y
# Initial guess for the parameters
# p0 = [-15., 0.8, 0., -1.]
A, B, C = sinefit.sinefit3(data, Ts, w0)
p0 = [A, B, C, w0];
# Time serie
n = data.size
Tx = arange(n) * Ts
p1, success = optimize.leastsq(errfunc, p0[:], args=(Tx, data))
return p1[3]
if __name__ == "__main__":
pass
from FFTSignal import FFTSignal
from numpy import *
__doc__ = """
Interface for the adc ADC characterization module
This is a preliminary spec of the interface to be expected from the
GUI of the ADC characterization environment
Some notes:
- A Signal object packs all the info we may need from a signal, to
wit,
* its raw data vector and, implicitly,
* raw data vector length == number of samples,
* number of bits (size of the integer samples in raw data vector)
* sampling rate
"""
class Signal(object):
"""a representation of a time-domain sampled signal
"""
MAXRES = 512 # maximum tolerable resolution for a histogram
cacheStatus = False
def __init__(self, nbits, rate, data):
"""initialize a signal object
nbits: bit width of the sample values
rate: sampling rate of sample production
data: an array of samples (usually nbits-long words, stored in
a numpy array)
"""
self.nbits = nbits
self.rate = rate
self.data = data
# Also known as Mt
self.nsamples = len(data)
self.precalculate()
def precalculate(self):
self.cachedIdealHistogram = self.ideal_histogram()
self.cachedRealHistogram = self.histogram()
self.cachedDNL, self.cachedDNLDiscard = self.newDNL()
self.cachedINL, self.cachedINLDiscard = self.INL()
self.cachedThSNR = 6.02 * self.nbits + 1.76
self.cachedProcessGain = 10.0 * log10(self.nbits / 2.0)
self.cacheStatus = True
def histogram_resolution(self):
bins = 2**self.nbits
return 512 if bins > self.MAXRES else bins
def histogram(self):
"""Compute histogram of a sampled signal
The number of bins in the histogram is implicitly given by the number
of bits of the samples: 2**signal.nbits bins.
returns: an array of 2**signal.nbits numbers (frequencies)
"""
# bins = self.histogram_resolution()
# hist, bins = histogram(array(self.data), bins)
bins = 2**self.nbits
hist, discard = histogram(array(self.data), bins)
# why it discards the first and the last values?
# maybe to get DNL easily?
return hist[1:-1]
# why?
def _ideal_histogram(self):
"""Produce an ideal vector of frequencies (histogram) for the
nsamples samples of a perfect nbits ADC. Mostly for auxiliary and
display purposes
returns: an array of 2**signal.nbits numbers (frequencies)
"""
Mt = len(self.data)
A = sin(pi/2 * Mt / (Mt + self.data[0] + self.data[-1]))
range = 2**self.nbits
midrange = range/2
n = arange(1, range-1)
p = arcsin(A/midrange * (n - midrange))
q = arcsin(A/midrange * (n - 1 - midrange))
p = (p - q) / pi
return Mt * p
def ideal_histogram(self):
"""Produce an ideal vector of frequencies (histogram) for the
nsamples samples of a perfect nbits ADC. Mostly for auxiliary and
display purposes
returns: an array of 2**signal.nbits numbers (frequencies)
"""
Mt = self.nsamples
A = sin(pi/2 * Mt / (Mt + self.data[0] + self.data[-1]))
range = 2**self.nbits
midrange = range/2
n = arange(1, range-1)
p = arcsin(A/midrange * (n - midrange))
q = arcsin(A/midrange * (n - 1 - midrange))
p = (p - q) / pi
# already commented out..
# t = linspace(-1, 1, 2**self.nbits)[1:-1]
# print sum(1/pi * 1/sqrt(1-t**2) / 2**self.nbits * Mt)
# return 1/pi * 1/sqrt(1-t**2) / 2**self.nbits * Mt
return Mt * p
def DNL(self):
"""Compute differential non-linearity vector for a given time-domain
signal
returns: a pair (dnl, total) where
- dnl is an array of 2**signal.nbits real values and
- total is a real value (computed from dnl)
"""
ideal = self.ideal_histogram()
real = self.histogram()
print size(ideal), size(real)
dnl = real/ideal - 1
return dnl, max(abs(dnl))
def newDNL(self):
dnl = (self.cachedRealHistogram/self.cachedIdealHistogram) -1
print "newDNL"
return dnl, max(abs(dnl))
def INL(self):
"""Compute integral non-linearity vector for a given time-domain signal
returns: a pair (inl, total) where
- inl is an array of 2**signal.nbits real values and
- total is a real (computed from inl)
"""
dnl, discard = self.DNL()
inl = cumsum(dnl)
return inl, max(abs(inl))
def newINL(self):
inl = cumsum(self.cachedINL)
print "newINL"
return inl, max(abs(inl))
def FFT(self, navg=1, window='No window'):
"""Compute the amplitudes (in dB) of the FFT of signal, averaging navg
slices of it and applying window to it
navg: number of signal slices to average
window: a value from a finite list of windows defined in window_types
returns: an FFTSignal object
"""
print 'using window ', window
# win = FFTSignal.window_function[window](self.nsamples)
# dft = 10*log10(abs(fft.rfft(self.data * win)))
output = FFTSignal(self, True)
output.selectedWindow = window
print "fft created"
return output
def makesine(samples, periods, bits, amplitude=1, noise=0):
t = arange(samples)/float(samples)
sine = amplitude * sin(2*pi*periods*t)
sine += noise * ((t % 0.02) / 0.02 - 0.01)
sine = (sine * 2**bits + 0.5).astype(int)
place(sine, sine >= 2**bits, 2**bits)
place(sine, sine <= -2**bits, -2**bits)
out = file('data', 'w')
for datum in sine:
out.write(str(datum) + '\n')
out.close()
return sine
if __name__ == '__main__':
from matplotlib import pyplot
bits = 12
makesine(20000, 20, bits, 1.1)
f = [ int(sample) for sample in file('data')]
s = Signal(nbits = bits, rate = 123, data = f)
ideal = s.ideal_histogram()
real = s.histogram()
pyplot.plot(ideal, label = 'ideal hist.')
pyplot.plot(real, label = 'real his.')
dnl = real/ideal-1
pyplot.plot(dnl, label = 'DNL')
pyplot.show()
pyplot.plot(cumsum(dnl), label = 'INL')
pyplot.show()
print dnl[0:5], dnl[-5:]
__author__="federico"
__date__ ="$Jul 11, 2011 5:21:40 PM$"
\ No newline at end of file
#! /bin/sh
# coding: utf-8
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
def ideal_adc(vals, N, Vmin, Vmax, mid_tread=False):
"""computes an ideal ADC characteristic function
vals is a vector of values to be converted by an
ideal ADC whose characteristic is defined by the Vmin and Vmax full
scale range, and the N number of bits. The optional mid_tread
parameter allow to use the mid_tread convention as per
IEEE 1241 sec. 1.2.
A vector of converted vals object is returned, taking into account
saturation at both ends of scale
"""
FS = float(Vmax - Vmin)
lsb = FS / 2**N
if mid_tread:
vals = vals + lsb/2
t = (vals - Vmin) / lsb
print vals, t
code = array(t, dtype=int)
code[code<0] = 0
code[code>(2**N-1)] = 2**N-1
return code
if __name__ == '__main__':
from numpy import *
from matplotlib import pyplot as pl
v = linspace(-1.0, 11, 3000)
print v
my_adc = ideal_adc(v, 4, 0, 10.0, mid_tread=1)
pl.plot(v, my_adc)
pl.show()
#! /bin/sh
# coding: utf-8
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
from numpy import linspace, pi, sin
def sampled_sine(A, w, phase, C, fs, t0=0.0, t1=1.0):
samples = (t1 - t0) * fs
time = linspace(t0, t1, samples, endpoint=False)
wave = A * sin(w * time + phase) + C
return wave
if __name__ == '__main__':
from matplotlib import pyplot as pl
s = sampled_sine(1, 2*pi, 0, 0, 8000)
pl.plot(s)
print 'samples = %d' % s.size
pl.show()
#! /bin/sh
# coding: utf-8
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
from numpy import *
from matplotlib import pyplot as pl
def sinefit3(samples, sample_t, w0):
"""fit a sampled wave to a sine of kwnown frequency
This routine implements the algoritm of IEEE 1241, sect. 4.1.4.1.,
fitting a sine of given frequency w0 to the samples given, which are
assumed equally spaced by a sample period of sample_t
a0 cos(w0 t + theta) + b0 sin(w0 t + theta) + c0
The return value is the triple (a0, b0, c0)
"""
n = samples.size
t = w0 * arange(n) * sample_t
D0T = matrix(vstack([cos(t), sin(t), ones(n)]))
D0 = D0T.T
x0 = (D0T * D0).I * D0T * matrix(samples).T
return array(x0).reshape((3,))
# page 27
def sinefit4(samples, sample_t, w0, tol=1e-4):
"""fit a sampled wave to a sine of unknown frequency
This routine implements the algorithm of IEEE 1241, sect. 4.1.4.3.,
fitting a sine wave the the samples given, which are assumed equally
spaced in time by a sample period sample_t.
a0 cos(w0 t + theta) + b0 sin(w0 t + theta) + c0
An initial frequency estimate w0 is required to start, and an
optional relative error in frequency is admitted (1e-4 by default)
The return value is the quadruple (w0, a0, b0, c0)
"""
a0, b0, c0 = sinefit3(samples, sample_t, w0)
deltaw0 = 0
print "Params: %f, %f" % (sample_t, w0)
print "a0\tb0\tc0\tdeltaw0"
print "%f\t%f\t%f\t%f" % (a0, b0, c0, deltaw0)
# we could move n and t definitions outside the while loop
while True:
# update freq
w0 += deltaw0
# get our array t1 -> tn
n = samples.size
t = arange(n) * sample_t
# multiply it with the frequency
w0t = w0 * t
# create the nX4 matrix and its transpose
D0T = matrix(vstack([cos(w0t), sin(w0t), ones(n),
-a0*t*sin(w0t) + b0*t*cos(w0t)]))
D0 = D0T.T
# get the solution
x0 = (D0T * D0).I * D0T * matrix(samples).T
x0 = array(x0).reshape((4,))
a0, b0, c0, deltaw0 = x0
print "%f\t%f\t%f\t%f" % (a0, b0, c0, deltaw0)
if abs(deltaw0/w0) < tol:
print "%f/%f < %f" % (deltaw0, w0, tol)
return (w0+deltaw0, a0, b0, c0)
if __name__ == '__main__':
a = 0.0
b = 10*pi
samples = 2000
w = 1
samp = sin(w*linspace(a, b, samples))
# fitting with known frequency
a0, b0, c0 = sinefit3(samp, (b-a)/samples, w)
ampl = hypot(a0, b0)
theta = arctan2(a0, b0)
print '--- %.2f * sin(%.5f t + %.2f) + %.2f' % (ampl, w, theta, c0)
# fitting with free frequency
err = 1+0.1*(random.random_sample() - 0.05)
w0, a0, b0, c0 = sinefit4(samp, (b-a)/samples, w*err)
ampl = hypot(a0, b0)
theta = arctan2(a0, b0)
print '%.2f * sin(%.5f t + %.2f) + %.2f' % (ampl, w0, theta, c0)
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="Federico Asara"
__date__ ="$Jul 11, 2011 2:39:38 PM$"
__doc__= """This module offers the Signal class, which calculate a variety of
parameters of an ADC output signal. """
# We use its array for data representation, also we use its FFT implementation
from numpy import *
from matplotlib import pyplot
# For frequency detection and unit-testin
import Sinefit
# We need all the window functions
import WindowFunction
def dB(x): return 20*log10(x)
class WindowedSignal:
"""Container object that will hold information about a DFT of a signal
multiplied by some window function"""
pass
class Signal(object):
"""A class that represent a time-domain sampled signal, and also holds many
WindowedSignal that represent the signal in frequency domain.
"""
# maximum tolerable resolution for a histogram
MAXRES = 512
# indicates the selected window
selectedWindow = "No window"
windowed = {}
def __getitem__(self, what):
"""Get a windowed signal item for a particular window function"""
if what is None:
what = self.selectedWindow
# just to be sure
if what not in WindowFunction.windows.keys():
what = self.WindowFunction.windows.keys()[0]
return self.windowed[what]
def __init__(self, nbits, rate, data):
"""initialize a signal object
nbits: bit width of the sample values
rate: sampling rate of sample production
data: an array of samples (usually nbits-long words, stored in
a numpy array)
"""
self.nbits = nbits
self.rate = rate
self.data = array(data, dtype=float) # just to be sure
print mean(self.data)
self.data -= (max(self.data) +min(self.data))/2.
self.nsamples = len(data)
print self.data
self.precalculateAll()
def precalculateAll(self):
"""Evaluates all the parameters of the signal, and also call the
precalculate method for each window function we know."""
# First of all evaluate the histograms
self.idealHistogram = [0]#self._ideal_histogram()
self.realHistogram = [0]#self._histogram()
# Then evaluate DNL and INL
self.DNL, self.maxDNL = [0], 0#self._DNL()
self.INL, self.maxINL = [0], 0#self._INL()
# ..theoretical SNR and process gain
self.thSNR = 6.02 * self.nbits + 1.76
self.processGain = 10.0 * log10(self.nbits / 2.0)
# This dictionary will hold all the WindowedSignal objects
self.windowed = {}
for i in WindowFunction.windows.keys():
print
print "Precalculating fft data for", i
print
self.windowed[i] = self.precalculate(i)
def precalculate(self, windowName):
"""Evaluates all the parameters for a particular window function.
windowName: a string containing the name of the window function
Returns a WindowFunction object, or None if windowName is not valid."""
# windows are ok right now
# wrong values: THD
# wrong formatting output: THD
if windowName not in WindowFunction.windows.keys():
return None # and maybe launch an exception
output = WindowedSignal()
window = WindowFunction.windows[windowName]
print "Filtering data"
output.data = self.data * window[self.nsamples]
print "FFT-ing data"
output.dft = abs(fft.rfft(output.data))
print output.dft
output.udft = output.dft[0:len(output.dft):10]
print "log10-ing data"
output.ldft = 10*log10(output.dft)
output.ludft = 10*log10(output.udft)
print "Guessing w0"
# sinusoid frequency detection
w0, w0index = Sinefit.sineguess(output.ldft, self.rate, self.nsamples)
freqSample = 2*pi*self.rate
ratio = w0 / w0index
print "Samples/period:",self.nsamples / w0index
under = self.data[::10]
# A, B, C, output.w0 = Sinefit.sinefit4(self.data[:500], freqSample, w0)
A, B, C, output.w0 = Sinefit.sinefit4(under, freqSample, w0)
output.w0 /= 10
amplitude = sqrt(A**2 + B**2)
phase = arctan(A/B)
# let's go through harmonic peaks. we will produce a generator
# this way
def adjust(data, fs):
while data >= fs:
data -= fs
if data >= fs/2.:
data = fs -data;
return data
def getGenerator(start, end):
indexes = (adjust(x * w0index, self.nsamples -1) for x in xrange(start , end))
return ((i, ratio*i, output.dft[i]) for i in indexes)
def getLogGenerator(start, end):
indexes = (adjust(x * w0index, self.nsamples -1) for x in xrange(start , end))
return ((i, ratio*i, 10*log10(output.dft[i])) for i in indexes)
output.harmonicPeaksGenerator = getGenerator
output.logHarmonicPeaksGenerator = getLogGenerator
# THD seems easy..
tenHarmonics = list(getGenerator(2, 10))
thindex = vstack(map(lambda x: x[0], tenHarmonics))
thvalues = vstack(map(lambda x: x[2], tenHarmonics))
tenHarmonicsValues = array(map(lambda x: x[2], tenHarmonics))
rssHarmonics = sqrt(sum(tenHarmonicsValues**2))
output.THD = dB(output.dft[w0index]/rssHarmonics)
# we need the avg of HDs
avgHarmonics = mean(tenHarmonicsValues)
# we also need the noise floor: the average of all components below
# avgHarmonics
filteredNoise = where(output.dft < avgHarmonics, output.dft, 0)
output.noiseFloor = dB(mean(filteredNoise))
output.signalPower = sum(abs(output.dft)**2)/self.nsamples
#sum(where(output.dft < output.noiseFloor, output.dft, 0)**2)/(self.nsamples**2 )
# thSin = C + Sinefit.makesine(self.nsamples, w0index, self.nbits, 1 )
time = arange(0, self.nsamples, dtype=float)/self.nsamples # linspace(0, 1, self.nsamples, endpoint=False )
print A, B, C, phase, amplitude
thSin = C + amplitude * sin(w0index*2*pi*time + phase)
tmp = file("/tmp/something", 'w')
tmp.writelines("%f\n" % i for i in thSin)
tmp.close()
noise = self.data - thSin
output.noisePower = sum(abs(fft.rfft(noise))**2)/self.nsamples
# now we can evaluate SNR = max component - noise floor - process gain
# output.SNR = output.dft[w0index] -output.noiseFloor -self.processGain
output.SNR = dB(output.signalPower/output.noisePower)
# now it's time for SINAD
clean = array(output.dft)
clean[0] = 0
clean[w0index] = 0
output.SINAD = dB(output.dft[w0index]/sqrt(sum(clean**2)))
# ENOB
fsr = 2**(self.nbits -1)
ra = (max(output.data) - min(output.data))/2.0
factor = dB(fsr/ra)
output.ENOB = (output.SINAD - 1.76 + factor) / 6.02
# SFDR - should I use dB(x) instead of 10log10 ?
output.maxPeak = 10*log10(output.dft[w0index])
secondPeak = max(thvalues)
output.secondPeak = 10*log10(secondPeak)
output.SFDR = 10*log10(output.dft[w0index] - secondPeak)[0] #10*log10
print "Sampling frequency = %.6f rad/s" % freqSample
print "Input frequency detected = %.6f rad/s" % output.w0
print "THD = %g dB" % output.THD
print "Noise floor = %.6f dB" % output.noiseFloor
print "Signal power = %.f W" % output.signalPower
print "Noise power = %.f W" % output.noisePower
print "SNR = %.6f dB" % output.SNR
print "SINAD = %.6f dB" % output.SINAD
print "ENOB = %.f b" % output.ENOB
print "SFDR = %.6f dBc" % output.SFDR
return output
def histogram_resolution(self):
bins = 2**self.nbits
return 512 if bins > self.MAXRES else bins
def _histogram(self):
"""Compute histogram of a sampled signal
The number of bins in the histogram is implicitly given by the number
of bits of the samples: 2**signal.nbits bins.
returns: an array of 2**signal.nbits numbers (frequencies)
"""
bins = 2**self.nbits
hist, discard = histogram(self.data, bins)
return hist[1:-1]
def _ideal_histogram(self):
"""Produce an ideal vector of frequencies (histogram) for the
nsamples samples of a perfect nbits ADC. Mostly for auxiliary and
display purposes
returns: an array of 2**signal.nbits numbers (frequencies)
"""
Mt = self.nsamples
A = sin(pi/2 * Mt / (Mt + self.data[0] + self.data[-1]))
range = 2**self.nbits
midrange = range/2
n = arange(1, range-1)
p = arcsin(A/midrange * (n - midrange))
q = arcsin(A/midrange * (n - 1 - midrange))
p = (p - q) / pi
return Mt * p
def _DNL(self):
dnl = (self.realHistogram/self.idealHistogram) -1
return dnl, max(abs(dnl))
def _INL(self):
inl = cumsum(self.DNL)
return inl, max(abs(inl))
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="Federico Asara"
__date__ ="$Jul 11, 2011 2:39:38 PM$"
__doc__= """This module offers the Signal class, which calculate a variety of
parameters of an ADC output signal. """
# We use its array for data representation, also we use its FFT implementation
from numpy import *
# For frequency detection and unit-testin
import Sinefit
# We need all the window functions
import WindowFunction
def dB(x): return 20*log10(x)
class WindowedSignal:
"""Container object that will hold information about a DFT of a signal
multiplied by some window function"""
pass
class Signal(object):
"""A class that represent a time-domain sampled signal, and also holds many
WindowedSignal that represent the signal in frequency domain.
"""
# maximum tolerable resolution for a histogram
MAXRES = 512
# indicates the selected window
selectedWindow = "No window"
def __getitem__(self, what):
"""Get a windowed signal item for a particular window function"""
if what is None:
what = self.selectedWindow
# just to be sure
if what not in WindowFunction.windows.keys():
what = self.WindowFunction.windows.keys()[0]
return self.windowed[what]
def __init__(self, nbits, rate, data):
"""initialize a signal object
nbits: bit width of the sample values
rate: sampling rate of sample production
data: an array of samples (usually nbits-long words, stored in
a numpy array)
"""
self.nbits = nbits
self.rate = rate
self.data = array(data) # just to be sure
self.nsamples = len(data)
self.precalculateAll()
def precalculateAll(self):
"""Evaluates all the parameters of the signal, and also call the
precalculate method for each window function we know."""
# First of all evaluate the histograms
self.idealHistogram = self._ideal_histogram()
self.realHistogram = self._histogram()
# Then evaluate DNL and INL
self.DNL, self.maxDNL = self._DNL()
self.INL, self.maxINL = self._INL()
# ..theoretical SNR and process gain
self.thSNR = 6.02 * self.nbits + 1.76
self.processGain = 10.0 * log10(self.nbits / 2.0)
# This dictionary will hold all the WindowedSignal objects
self.windowed = {}
for i in WindowFunction.windows.keys():
print
print "Precalculating fft data for", i
print
self.windowed[i] = self.precalculate(i)
def precalculate(self, windowName):
"""Evaluates all the parameters for a particular window function.
windowName: a string containing the name of the window function
Returns a WindowFunction object, or None if windowName is not valid."""
# windows are ok right now
# wrong values: THD
# wrong formatting output: THD
if windowName not in WindowFunction.windows.keys():
return None # and maybe launch an exception
output = WindowedSignal()
window = WindowFunction.windows[windowName]
output.data = self.data * window[self.nsamples]
output.dft = abs(fft.rfft(output.data))
output.ldft = 10*log10(output.dft)
# sinusoid frequency detection
w0, w0index = Sinefit.sineguess(output.ldft, self.rate, self.nsamples)
freqSample = 2*pi*self.rate
ratio = w0 / w0index
output.w0 = Sinefit.sinefit4(self.data, freqSample, w0)
# let's go through harmonic peaks. we will produce a generator
# this way
def adjust(data, fs):
while data >= fs:
data -= fs
if data >= fs/2.:
data = fs -data;
return data
def getGenerator(start, end):
indexes = (adjust(x * w0index, self.nsamples -1) for x in xrange(start , end))
return ((i, ratio*i, output.dft[i]) for i in indexes)
def getLogGenerator(start, end):
indexes = (adjust(x * w0index, self.nsamples -1) for x in xrange(start , end))
return ((i, ratio*i, 10*log10(output.dft[i])) for i in indexes)
output.harmonicPeaksGenerator = getGenerator
output.logHarmonicPeaksGenerator = getLogGenerator
# THD seems easy..
tenHarmonics = list(getGenerator(2, 10))
thindex = vstack(map(lambda x: x[0], tenHarmonics))
thvalues = vstack(map(lambda x: x[2], tenHarmonics))
tenHarmonicsValues = array(map(lambda x: x[2], tenHarmonics))
rssHarmonics = sqrt(sum(tenHarmonicsValues**2))
output.THD = dB(output.dft[w0index]/rssHarmonics)
# we need the avg of HDs
avgHarmonics = mean(tenHarmonicsValues)
# we also need the noise floor: the average of all components below
# avgHarmonics
filteredNoise = where(output.dft < avgHarmonics, output.dft, 0)
output.noiseFloor = dB(mean(filteredNoise))
output.signalPower = sum(abs(output.dft)**2)/self.nsamples
#sum(where(output.dft < output.noiseFloor, output.dft, 0)**2)/(self.nsamples**2 )
avg = mean(self.data)
maxval = max(self.data)
minval = min(self.data)
amplitude = 1
periods = w0index
thSin = avg + Sinefit.makesine(self.nsamples, periods, self.nbits, amplitude)
noise = self.data - thSin
output.noisePower = sum(abs(fft.rfft(noise))**2)/self.nsamples
# now we can evaluate SNR = max component - noise floor - process gain
# output.SNR = output.dft[w0index] -output.noiseFloor -self.processGain
output.SNR = dB(output.signalPower/output.noisePower)
# now it's time for SINAD
clean = array(output.dft)
clean[0] = 0
clean[w0index] = 0
output.SINAD = dB(output.dft[w0index]/sqrt(sum(clean**2)))
# ENOB
fsr = 2**(self.nbits -1)
ra = (max(output.data) - min(output.data))/2.0
factor = dB(fsr/ra)
output.ENOB = (output.SINAD - 1.76 + factor) / 6.02
# SFDR - should I use dB(x) instead of 10log10 ?
output.maxPeak = 10*log10(output.dft[w0index])
secondPeak = max(thvalues)
output.secondPeak = 10*log10(secondPeak)
output.SFDR = 10*log10(output.dft[w0index] - secondPeak) #10*log10
print "Sampling frequency = %.6f rad/s" % freqSample
print "Input frequency detected = %.6f rad/s" % output.w0
print "THD = %g dB" % output.THD
print "Noise floor = %.6f dB" % output.noiseFloor
print "Signal power = %.f W" % output.signalPower
print "Noise power = %.f W" % output.noisePower
print "SNR = %.6f dB" % output.SNR
print "SINAD = %.6f dB" % output.SINAD
print "ENOB = %.f b" % output.ENOB
print "SFDR = %.6f dBc" % output.SFDR
return output
def histogram_resolution(self):
bins = 2**self.nbits
return 512 if bins > self.MAXRES else bins
def _histogram(self):
"""Compute histogram of a sampled signal
The number of bins in the histogram is implicitly given by the number
of bits of the samples: 2**signal.nbits bins.
returns: an array of 2**signal.nbits numbers (frequencies)
"""
bins = 2**self.nbits
hist, discard = histogram(self.data, bins)
return hist[1:-1]
def _ideal_histogram(self):
"""Produce an ideal vector of frequencies (histogram) for the
nsamples samples of a perfect nbits ADC. Mostly for auxiliary and
display purposes
returns: an array of 2**signal.nbits numbers (frequencies)
"""
Mt = self.nsamples
A = sin(pi/2 * Mt / (Mt + self.data[0] + self.data[-1]))
range = 2**self.nbits
midrange = range/2
n = arange(1, range-1)
p = arcsin(A/midrange * (n - midrange))
q = arcsin(A/midrange * (n - 1 - midrange))
p = (p - q) / pi
return Mt * p
def _DNL(self):
dnl = (self.realHistogram/self.idealHistogram) -1
return dnl, max(abs(dnl))
def _INL(self):
inl = cumsum(self.DNL)
return inl, max(abs(inl))
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 11, 2011 2:40:18 PM$"
from numpy import *
from scipy import optimize
def sineguess(dft, rate, nsamples):
w0index = argmax(dft)
w0 = w0index * (2 * pi * rate) / nsamples
return w0, w0index
def sinefit3(samples, sample_t, w0):
"""fit a sampled wave to a sine of kwnown frequency
This routine implements the algoritm of IEEE 1241, sect. 4.1.4.1.,
fitting a sine of given frequency w0 to the samples given, which are
assumed equally spaced by a sample period of sample_t
a0 cos(w0 t + theta) + b0 sin(w0 t + theta) + c0
The return value is the triple (a0, b0, c0)
"""
n = samples.size
t = w0 * arange(n) * sample_t
D0T = matrix(vstack([cos(t), sin(t), ones(n)]))
D0 = D0T.T
x0 = (D0T * D0).I * D0T * matrix(samples).T
return array(x0).reshape((3,))
def sinefit4(data, Ts, w0):
# Target function
fitfunc = lambda p, x: p[0]*cos(p[3]*x) +p[1]*sin(p[3]*x) +p[2]
# Distance to the target function
errfunc = lambda p, x, y: fitfunc(p, x) - y
# Initial guess for the parameters
A, B, C = sinefit3(data, Ts, w0)
p0 = [A, B, C, w0];
# Time serie
n = data.size
Tx = arange(n) * Ts
p1, success = optimize.leastsq(errfunc, p0[:], args=(Tx, data))
return p1
def makesine(samples, periods, bits, amplitude=1, noise=0):
t = arange(samples)/float(samples)
sine = amplitude * sin(2*pi*periods*t)
sine += noise * ((t % 0.02) / 0.02 - 0.01)
sine = (sine * 2**bits + 0.5).astype(int)
place(sine, sine >= 2**bits, 2**bits)
place(sine, sine <= -2**bits, -2**bits)
return sine
if __name__ == "__main__":
sine = makesine(20000, 20, 12, 443, 10)
print sine
dft = 10*log10(abs(fft.rfft(sine)))
w0, w0index = sineguess(dft, 133, 20000)
out = sinefit4(sine, 2*pi*(133**-1), w0)
print w0, w0index, out
orig = makesine(20000, 20, 12, 443, 0);
noise = sine -orig
ndft = 10*log10(abs(fft.rfft(noise)))
print noise
snr = sum(dft**2)/sum(ndft**2)
print "SNR:", snr
\ No newline at end of file
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 11, 2011 2:40:18 PM$"
from numpy import *
from scipy import optimize
def sineguess(dft, rate, nsamples):
w0index = argmax(dft)
w0 = w0index * (2 * pi * rate) / nsamples
return w0, w0index
def sinefit3(samples, sample_t, w0):
"""fit a sampled wave to a sine of kwnown frequency
This routine implements the algoritm of IEEE 1241, sect. 4.1.4.1.,
fitting a sine of given frequency w0 to the samples given, which are
assumed equally spaced by a sample period of sample_t
a0 cos(w0 t + theta) + b0 sin(w0 t + theta) + c0
The return value is the triple (a0, b0, c0)
"""
n = samples.size
t = w0 * arange(n) * sample_t
D0T = matrix(vstack([cos(t), sin(t), ones(n)]))
D0 = D0T.T
x0 = (D0T * D0).I * D0T * matrix(samples).T
return array(x0).reshape((3,))
def sinefit4(data, Ts, w0):
# Target function
fitfunc = lambda p, x: p[0]*cos(p[3]*x) +p[1]*sin(p[3]*x) +p[2]
# Distance to the target function
errfunc = lambda p, x, y: fitfunc(p, x) - y
# Initial guess for the parameters
A, B, C = sinefit3(data, Ts, w0)
p0 = [A, B, C, w0];
# Time serie
n = data.size
Tx = arange(n) * Ts
p1, success = optimize.leastsq(errfunc, p0[:], args=(Tx, data))
return p1[3]
def makesine(samples, periods, bits, amplitude=1, noise=0):
t = arange(samples)/float(samples)
sine = amplitude * sin(2*pi*periods*t)
sine += noise * ((t % 0.02) / 0.02 - 0.01)
sine = (sine * 2**bits + 0.5).astype(int)
place(sine, sine >= 2**bits, 2**bits)
place(sine, sine <= -2**bits, -2**bits)
return sine
if __name__ == "__main__":
sine = makesine(20000, 20, 12, 443, 10)
print sine
dft = 10*log10(abs(fft.rfft(sine)))
w0, w0index = sineguess(dft, 133, 20000)
out = sinefit4(sine, 2*pi*(133**-1), w0)
print w0, w0index, out
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 11, 2011 2:42:12 PM$"
from numpy import ones, bartlett, blackman, hamming, hanning, kaiser
class WindowFunction(object):
def __init__(self, name, function, parameters):
self.name = name
self.function = function
self.parameters = parameters
def get(self, key, *parameters):
if len(parameters) != len(self.parameters):
parameters = self.dump()
return self.function(key, *parameters)
def set(self, what, how):
for i in xrange(len(self.parameters)):
if self.parameters[i] == what:
self.parameters[i +1] = how
return
def get(self, what):
for i in xrange(len(self.parameters)):
if self.parameters[i] == what:
return self.parameters[i +1]
def dump(self):
return [self.parameters[i] for i in xrange(len(self.parameters)) if (i % 2) == 1]
def names(self):
return [self.parameters[i] for i in xrange(len(self.parameters)) if (i % 2) == 0]
def __getitem__(self, key):
return self.function(key, *self.dump())
windows = {
'No window' : WindowFunction('No window', ones, [])#,
#'Bartlett' : WindowFunction('Bartlett', bartlett, []),
#'Blackman' : WindowFunction('Blackman', blackman, []),
#'Hamming' : WindowFunction('Hamming', hamming, []),
#'Hanning' : WindowFunction('Hanning', hanning, []),
#'Kaiser' : WindowFunction('Kaiser', kaiser, ["Beta", 14])
}
if __name__ == "__main__":
for name in windows:
print "Using window", name
print windows[name][15]
\ No newline at end of file
__author__="federico"
__date__ ="$Jul 11, 2011 2:39:08 PM$"
\ No newline at end of file
import wx
from wx.lib.pubsub import Publisher as pub
from Model import Model
from View import View
class Controller:
""" a 'middleman' between the View (visual aspects) and the Model (information) of the application.
It ensures decoupling between both.
"""
def __init__(self, app):
# initialize the model and view
# * The model handles all the data, and signal-related operations
# * The view handles all the data visualization
self.model = Model()
self.view = View()
# subscribe to messages sent by the view
pub.subscribe(self.parse_file, "FILE PATH CHANGED")
pub.subscribe(self.reprocess_fft, "FFT CONTROLS CHANGED")
# subscribe to messages sent by the model
pub.subscribe(self.signal_changed, "SIGNAL CHANGED")
pub.subscribe(self.signal_changed, "FFT CHANGED")
self.view.Show()
def parse_file(self, message):
"""
Handles "FILE PATH CHANGED" messages, send by the View. It tells the model to parse a new file.
message.data should contain the path of the new file
"""
#try:
self.model.parse_file(message.data)
#except Exception as exception:
# self.view.show_exception('Error reading file', 'The following error happened while reading the file:\n%s' % str(exception))
# raise exception
def reprocess_fft(self, message):
"""
Handler "FFT CONTROLS CHANGED" messages from the View. It tells the model to re-process the fft.
message.data should contain the array [window, slices, max_peaks]
"""
self.model.reprocess_fft(* message.data)
def signal_changed(self, message):
"""
Handles "SIGNAL CHANGED" messages sent by the model. Tells the view to update itself.
message is ignored
"""
self.view.signal_changed(self.model)
def fft_changed(self, message):
"""
Handles "FFT CHANGED" messages sent by the model. Tells the view to update itself.
message is ignored
"""
self.view.fft_changed(self.model)
#! /usr/bin/python
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 15, 2011 2:20:38 PM$"
import wx, os
from wx.lib.pubsub import Publisher as pub
class FileChooser(wx.Panel):
def __init__(self, parent):
wx.Panel.__init__(self, parent)
# buttons / controls
path_label = wx.StaticText(self, -1, "Choose a signal file:")
self.path_ctrl = wx.TextCtrl(self)
self.open_button = wx.Button(self, -1, "...", wx.DefaultPosition, wx.Size(30,30))
self.parse_button = wx.Button(self, -1, "Parse file", wx.DefaultPosition, wx.Size(120,30))
# Bind actions to the two buttons
self.open_button.Bind(wx.EVT_BUTTON, self.open_dialog)
self.parse_button.Bind(wx.EVT_BUTTON, self.send_path_changed_message)
# sizer for the textbox and buttons
ctrls_sizer = wx.BoxSizer(wx.HORIZONTAL)
ctrls_sizer.Add(self.path_ctrl, 1, wx.EXPAND)
ctrls_sizer.Add(self.open_button, 0, wx.EXPAND)
ctrls_sizer.Add(self.parse_button, 0, wx.EXPAND)
vsizer = wx.BoxSizer(wx.VERTICAL)
vsizer.Add(path_label, 0.5, wx.EXPAND)
vsizer.Add(ctrls_sizer, 0.5, wx.EXPAND)
self.SetSizer(vsizer)
# jdgc: pa abreviar
# fasara: using os.path facilities
defaultPath = os.path.join(os.getcwd(), 'samples', 'data.txt')
self.path_ctrl.SetValue(defaultPath)
self.send_path_changed_message()
def send_path_changed_message(self, evt=None):
""" Sends a message to the Controller saying that the file path has been changed and it's ready to be reloaded.
Called right after using the '...' button and choosing a file, and also when pressing the 'Parse file' button.
"""
pub.sendMessage("FILE PATH CHANGED", self.path_ctrl.GetValue())
def open_dialog(self, evt=None):
""" Shows and controls the File Opening dialog """
prevPath = self.path_ctrl.GetValue()
if len(prevPath) > 0:
defaultDir = os.path.dirname(prevPath)
else:
defaultDir = os.getcwd()
dlg = wx.FileDialog(
self, message="Choose a file",
defaultDir=defaultDir,
defaultFile="",
wildcard="All files (*.*)|*.*",
style=wx.OPEN
)
# Show the dialog and retrieve the user response. If it is the OK response,
# Parse the data.
if dlg.ShowModal() == wx.ID_OK:
path = dlg.GetPath()
self.path_ctrl.SetValue(path)
self.send_path_changed_message()
dlg.Destroy()
import ConfigParser
import string
from wx.lib.pubsub import Publisher as pub
from SignalsProcessing.Signal import Signal
class Model:
""" Holds all the pertinent information regarding the signals.
It also caches the information so it isn't re-calculated unnecesarily (for example when re-dimensioning windows)
It communicates with the Controller by emiting messages.
"""
FFT_METHODS = ['process_gain', 'noise_floor', 'SFDR', 'SINAD', 'THD', 'SNR', 'ENOB']
def __init__(self):
self.signal = None
def parse_file(self, path, max_peaks=0, dB = None, time_domain = None):
""" Invoked when a new file needs to be parsed. Will raise an exception if the file contains syntax errors
Emits the message 'SIGNAL CHANGED' if the file was cached succesfully
"""
self.signal = None
print "parse file"
#try:
if True:
config = ConfigParser.RawConfigParser()
config.read(path)
nbits = config.getint('SIGNAL', 'nbits')
rate = config.getint('SIGNAL', 'rate')
dataString = config.get('SIGNAL', 'data').split('\n')
data = map(string.atoi, dataString)
print "setting up signal"
self.signal = Signal(nbits, rate, data)
print "signal done"
if True:
#finally:
self.cache_signal()
pub.sendMessage("SIGNAL CHANGED")
def reprocess_fft(self, window, slices, max_peaks):
""" Invoked when the FFT controls on tab 3 have been changed. re-calculates fft
Emits the message 'FFT CHANGED' when done calculating
"""
#if self.signal is None:
# self.fft_signal = None
#else:
# print 'window = ', window
# self.fft_signal = self.signal.FFT(slices, window)
self.cache_fft_signal(max_peaks)
pub.sendMessage("FFT CHANGED")
def cache_signal(self):
""" Invokes all methods in Signal and caches their result for later use (mainly window resizing)
Resets values to safe values (0, []) if there was an error while processing the signal definition file
"""
print "do nothing"
return
if self.signal is None:
self.data = []
self.INL, self.max_INL = [], 0
self.DNL, self.max_DNL = [], 0
self.histogram, self.ideal_histogram = [], []
else:
self.data = self.signal.data
self.INL, self.max_INL = self.signal.INL()
self.DNL, self.max_DNL = self.signal.DNL()
self.histogram, self.ideal_histogram = self.signal.histogram(), self.signal.ideal_histogram()
def cache_fft_signal(self, max_peaks = 0):
""" Invokes all methods in FFTSignal and caches their result for later use (mainly window resizing)
Resets values to safe values (0, []) if there was an error while processing the signal definition file
"""
print "do nothing"
return
if self.fft_signal is None:
self.fft = []
self.harmonic_peaks = []
for name in Model.FFT_METHODS:
setattr(self, name, 0)
else:
self.fft = self.fft_signal.fft
self.harmonic_peaks = self.fft_signal.harmonic_peaks(max_peaks)
for name in Model.FFT_METHODS:
setattr( self, name, getattr(self.fft_signal, name)() )
import wx
from wx.lib.pubsub import Publisher as pub
import matplotlib
matplotlib.use('WX')
from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
from matplotlib.figure import Figure
#FIXME remove when unused
from numpy import arange, sin, pi, cos, hstack
__doc__ = """
This module contains the definitions of all the graphs used in the app.
It is the only one interfacing with matplotlib.
"""
class Signal(FigureCanvas):
""" Represents the raw data from a Signal, right after being loaded from a file (Tab 1)
"""
def __init__(self, parent):
self.figure = Figure((5,3)) # 5,3 is the initial figure
FigureCanvas.__init__(self, parent, -1, self.figure)
self.axes = self.figure.add_subplot(111)
""" Specifies what has to be done if the signal is changed (the data has to be re-plotted)
"""
def update(self, data):
self.axes.clear()
self.axes.plot(arange(0,len(data),1) ,data, '-')
self.draw()
class DNL(Signal):
""" Represents the DNL plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
class INL(Signal):
""" Represents the INL plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
class Histogram(Signal):
""" Represents the histogram plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
""" Specifies what has to be done if the signal is changed (the data has to be re-plotted)
"""
def update(self, histogram, ideal_histogram):
self.axes.clear()
self.axes.plot(arange(0,len(histogram),1), histogram, '.')
self.axes.plot(arange(0,len(ideal_histogram),1), ideal_histogram, '.')
self.draw()
class FFT(Signal):
""" Represents the FFT plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
def update(self, fft, harmonic_peaks, maxPeak, sfdrPeak, noiseFloor):
self.axes.clear()
self.axes.grid(color='black', linestyle='-', linewidth=1, alpha=0.5)
self.axes.vlines(arange(0,len(fft),1), 0, fft)
if len(harmonic_peaks) != 0:
thindex = hstack(map(lambda x: x[0], harmonic_peaks))
thvalues = hstack(map(lambda x: x[2], harmonic_peaks))
self.axes.plot(thindex , thvalues, '.')
# plot SFDR
self.axes.axhline(maxPeak)
self.axes.axhline(sfdrPeak)
self.axes.axhspan(sfdrPeak, maxPeak, facecolor='blue', alpha=0.15)
# noise floor
self.axes.axhline(noiseFloor)
self.axes.axhspan(0, noiseFloor, facecolor='gray', alpha=0.15)
self.draw()
import wx
from wx.lib.pubsub import Publisher as pub
import matplotlib
matplotlib.use('WX')
from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
from matplotlib.figure import Figure
#FIXME remove when unused
from numpy import arange, sin, pi, cos, hstack
__doc__ = """
This module contains the definitions of all the graphs used in the app.
It is the only one interfacing with matplotlib.
"""
class Signal(FigureCanvas):
""" Represents the raw data from a Signal, right after being loaded from a file (Tab 1)
"""
def __init__(self, parent):
self.figure = Figure((5,3)) # 5,3 is the initial figure
FigureCanvas.__init__(self, parent, -1, self.figure)
self.axes = self.figure.add_subplot(111)
""" Specifies what has to be done if the signal is changed (the data has to be re-plotted)
"""
def update(self, data):
self.axes.clear()
self.axes.plot(arange(0,len(data),1) ,data, '-')
self.draw()
class DNL(Signal):
""" Represents the DNL plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
class INL(Signal):
""" Represents the INL plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
class Histogram(Signal):
""" Represents the histogram plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
""" Specifies what has to be done if the signal is changed (the data has to be re-plotted)
"""
def update(self, histogram, ideal_histogram):
self.axes.clear()
self.axes.plot(arange(0,len(histogram),1), histogram, '.')
self.axes.plot(arange(0,len(ideal_histogram),1), ideal_histogram, '.')
self.draw()
class FFT(Signal):
""" Represents the FFT plot
"""
def __init__(self, parent):
Signal.__init__(self, parent)
def update(self, fft, harmonic_peaks, maxPeak, sfdrPeak, noiseFloor):
self.axes.clear()
self.axes.grid(color='black', linestyle='-', linewidth=1, alpha=0.5)
self.axes.vlines(arange(0,len(fft),1), 0, fft)
if len(harmonic_peaks) != 0:
thindex = hstack(map(lambda x: x[0], harmonic_peaks))
thvalues = hstack(map(lambda x: x[2], harmonic_peaks))
self.axes.plot(thindex , thvalues, '.')
# plot SFDR
self.axes.axhline(maxPeak)
self.axes.axhline(sfdrPeak)
self.axes.axhspan(sfdrPeak, maxPeak, facecolor='blue', alpha=0.15)
# noise floor
self.axes.axhline(noiseFloor)
self.axes.axhspan(0, noiseFloor, facecolor='gray', alpha=0.15)
self.draw()
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 8, 2011 3:48:20 PM$"
import wx
from wx.lib.pubsub import Publisher as pub
class SliceChooser(wx.Panel):
def __init__(self, windows):
wx.Panel.__init__(self, windows)
self.slices_label = wx.StaticText(self, -1, "Slices:")
self.slices_ctrl = wx.SpinCtrl(self, -1, "1", wx.DefaultPosition, (100,30))
self.slices_ctrl.SetRange(1,10)
self.vsizer = wx.BoxSizer(wx.VERTICAL)
self.vsizer.Add(self.slices_label, 0.5, wx.EXPAND)
self.vsizer.Add(self.slices_ctrl, 0.5, wx.EXPAND)
def get(self):
return self.slices_ctrl.Get_Value()
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 8, 2011 3:48:20 PM$"
import wx
from wx.lib.pubsub import Publisher as pub
from WindowParameter import WindowParameter
class WindowPicker(wx.Panel):
def __init__(self, windows):
wx.Panel.__init__(self, windows)
# create a DropBox
self.windowLabel = wx.StaticText( self, -1, "Window:")
self.windows = windows
self.allWindows = windows.keys()
self.currentWindow = self.allWindows[0]
self.windowBox = wx.ComboBox(self, 500, self.currentWindow, wx.DefaultPosition, (120,30),
allWindows, wx.CB_DROPDOWN | wx.CB_READONLY | wx.CB_SORT )
self.slices_label = wx.StaticText(self, -1, "Slices:")
self.slices_ctrl = wx.SpinCtrl(self, -1, "1", wx.DefaultPosition, (100,30))
self.slices_ctrl.SetRange(1,10)
self.spinboxes = []
self.vsizer = wx.BoxSizer(wx.VERTICAL)
self.vsizer.Add(self.windowLabel, 0.5, wx.EXPAND)
self.vsizer.Add(self.windowBox, 0.5, wx.EXPAND)
self.vsizer.Add(self.windowBox, 0.5, wx.EXPAND)
self.updateView()
def updateView(self):
if self.windowBox.GetValue() == self.currentWindow:
return
self.currentWindow = self.windowBox.GetValue()
window = self.windows[self.currentWindow]
# hide and show spinners accordingly
for i in self.spinboxes:
self.vsizer.Remove(i)
self.spinboxes = [WindowParameter(self, i, window) for i in window.parameters]
for i in self.spinboxes:
self.vsizer.Add(i, 0.5, wx.EXPAND)
def get(self):
window = self.windows[self.currentWindow]
parameters = [i.get() for i in self.spinboxes]
return window, parameters
\ No newline at end of file
# To change this template, choose Tools | Templates
# and open the template in the editor.
__author__="federico"
__date__ ="$Jul 8, 2011 4:42:35 PM$"
import wx
from wx.lib.pubsub import Publisher as pub
class WindowParameter(wx.Panel):
def __init__(self, parent, name, window):
wx.Panel.__init__(self, parent)
self.label = wx.StaticText(self, -1, "%s:" % name)
self.spin = wx.SpinCtrl(self, -1, "1", wx.DefaultPosition, (100,30))
self.vsizer = wx.BoxSizer(wx.VERTICAL)
self.vsizer.Add(self.label, 0.5, wx.EXPAND)
self.vsizer.Add(self.spin, 0.5, wx.EXPAND)
self.SetSizer(vsizer)
def get(self):
return self.spin.GetValue()
__author__="federico"
__date__ ="$Jul 11, 2011 5:22:02 PM$"
\ No newline at end of file
import wxversion
# We need at least version 2.8
wxversion.ensureMinimal('2.8')
import wx
from Controller import Controller
# We can start importing modules
import wx
import UI.Controller
# Create an application
app = wx.App(False)
controller = Controller(app)
# Create our controller and bind it to the app
controller = UI.Controller.Controller(app)
# Start the main loop
app.MainLoop()
* parser error display freezes compiz. Shorten it to 25 lines or less in
* parser error display freezes compiz. Shorten it to 25 lines or less in
the popup
* automatically parse? An Open button and the parse button is not needed
......
<?xml version="1.0" encoding="UTF-8"?>
<project-private xmlns="http://www.netbeans.org/ns/project-private/1">
<editor-bookmarks xmlns="http://www.netbeans.org/ns/editor-bookmarks/1"/>
</project-private>
adctest.dir=.
java.lib.path=
main.file=adctest.py
platform.active=Python_2.7.1+
python.lib.path=
<?xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://www.netbeans.org/ns/project/1">
<type>org.netbeans.modules.python.project</type>
<configuration>
<data xmlns="http://nbpython.dev.java.net/ns/php-project/1">
<name>NewPythonProject</name>
<sources>
<root id="adctest.dir"/>
</sources>
<tests/>
</data>
</configuration>
</project>
Dependencies:
wxpython >= 2.8
numpy
matplotlib
Possible UI optimization:
- multithreading support
UI freezes during long computations
- graphic improvemet
We could store graphics in SVG format somewhere (tempdir), then we could just
show that SVG without waiting for pyplot
This is strictly related to Caching
Possible data management classes optimizations:
- Trasparent cache
We should implement the caching mechanism inside Signal and FFT Signal. This
way Model should not worry about caching at all, which is a better design IMHO.
- FFT Signal progresses
Several functions are not implemented in FFTSignal class.
from sys import argv
def adj(x):
return round(x/5.0 * (2**15))
src = open(argv[1], "r")
lines = ['\t' + i for i in src.readlines()]
src.close()
dest = open(argv[2], "w")
dest.write("[SIGNAL]\nnbits = 24\nrate = 1000000\n=ndata = ")
dest.writelines(lines)
dest.close()
-0.0139388910499
0.0756613283641
0.128759127602
0.169575809817
0.244099439553
0.306779721623
0.375048880639
0.425027260066
0.484185545203
0.535421457098
0.591754384546
0.619414359581
0.687367305494
0.730445183176
0.773202705162
0.815470355326
0.844011499105
0.875280563655
0.900282819658
0.926638791334
0.952639501932
0.966977862275
0.969674255213
1.00198945358
0.981079126722
1.007817229
0.997527816021
0.994226242771
0.988549746779
0.968969798004
0.961382701302
0.940630319908
0.921477072425
0.866496012822
0.843734203894
0.817917693515
0.767138189045
0.724103628976
0.676820098275
0.636729241868
0.587238486772
0.547481768193
0.483829299129
0.42766303586
0.361268122452
0.302087025471
0.258784141815
0.182385331159
0.113858537263
0.0668136363639
0.00927897621256
-0.0701149564905
-0.135020225869
-0.18493204501
-0.246262058332
-0.315120604573
-0.379369268698
-0.426213780707
-0.485519591091
-0.537120876094
-0.597022415245
-0.627126718401
-0.691513911185
-0.716619564037
-0.76833225414
-0.802031287717
-0.85216960123
-0.882476561339
-0.891286115932
-0.927169692758
-0.937679159348
-0.961516883106
-0.977622840092
-1.00073607962
-0.997377280586
-0.99872923133
-0.992706267503
-0.996058135359
-0.987415419325
-0.977516994888
-0.953092683785
-0.927896757013
-0.921581038383
-0.873139853179
-0.84335529775
-0.819544059154
-0.777849859037
-0.725212422197
-0.679382760045
-0.64160346947
-0.588515246869
-0.536935656739
-0.484613091446
-0.432437173482
-0.372268373293
-0.300333618858
-0.256697655616
-0.189626702835
-0.123154588856
-0.0660641664756
pure_sin.txt Pure sinewave, A = 2.9, fa = 100, phi = 0, 1MHz, t=0.02s
gaussian_100th.txt Same with 0.01 rms gaussian noise (SNR = 2e2)
gaussian_1000th.txt Same with 0.001 rms gaussian noise (SNR = 2e3)
incoherent.txt Pure sinewave as above, incoherently sampled (t = 0.025s)
incoherent_1000gaussian.txt Same as above with 0.001 rms gaussian noise (SNR = 2e3)
thd_44_7_times.txt Same as above with 0.01 and 0.02 harmonics (SNR = 44.7)
thd_0.1_0.08_n_0.005.txt With 0.1 and 0.08 harmonics and 0.005 gaussian noise
log How these were made
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
1 : import syn
2 : syn.synth(2.9, 100, 0.0, 1000000, 0, 0.02)
3 : s = +
4 : s = _
5 : s
6 : plot s
7 : from matplotlib import pyplot
8 : pyplot.plot(s)
9 : pyplot.show()
10: _ip.magic("save ")
11: help(save)
12: _ip.magic("save ")
13: save(s, 'pure_sin.txt')
14: save(s, 'pure_sin.txt')
15: from numpy import *
16: save(s, 'pure_sin.txt')
17: save('pure_sin.txt', s)
18: _ip.magic("edit pure_sin.txt.npy")
19: help(save)
20: out = file('pure_sin.txt', 'w')
21:
for i in s:
print >>out, i
22: out.close()
23: _ip.system("view pure*")
24: _ip.system("rm pure_sin.txt.npy")
25:
def w(s, f):
out = file(f, 'w')
for i in s:
print >>out, i
out.close()
26:
27:
28:
29: s
30: t
31: syn.synth(2.9, 100, 0.0, 1000000, 0, 0.02)
32: random.randn(len(s))
33:
34:
35: r = 2.9 * 0.001 * random.randn(len(s))
36: s+r
37: sn = s+r
38: w(s, 'gaussian_1000th.txt')
39: plot s
40: w(sn, 'gaussian_1000th.txt')
41: plot s
42: plot = pyplot.plot
43: plot(s)
44: sn = s+r
45: r = 2.9 * 0.01 * random.randn(len(s))
46: sn = s+r
47: w(sn, 'gaussian_100th.txt')
48: plot(sn)
49: r = 2.9 * 0.001 * random.randn(len(s))
50: sn = s+r
51: w(sn, 'gaussian_100th.txt')
52: plot(sn)
53: s = syn.synth(2.9, 100, 0.0, 1000000, 0, 0.025)
54: sincoh = s
55: w(sincoh, 'incoherent.txt')
56: sincohn = s + 0.001 * 2.9 * random.randn(len(s))
57: w(sincohn, 'incoherent_1000gaussian.txt')
58: plot(sincohn)
59: s = syn.synth(2.9, 100, 0.0, 1000000, 0, 0.025)
60: s = s + 0.01 * syn.synth(2.9, 300, 0.0, 1000000, 0, 0.025)
61: s = s + 0.02 * syn.synth(2.9, 500, 0.0, 1000000, 0, 0.025)
62: plot(s)
63: plot(s)
64: plot(s)
65: 1 / sqrt(5)
66: w(s, 'thd_44_7_times.txt')
67: s = syn.synth(2.9, 100, 0.0, 1000000, 0, 0.025)
68: s = s + 0.1 syn.synth(2.9, 500, 0.0, 1000000, 0, 0.025)
69: s = s + 0.1 * syn.synth(2.9, 500, 0.0, 1000000, 0, 0.025)
70: s = s + 0.08 * syn.synth(2.9, 800, 0.0, 1000000, 0, 0.025)
71: s = s + 0.005 * random.randn(len(s))
72: w(s,'thd_0.1_0.08_n_0.005.txt)
73: w(s,'thd_0.1_0.08_n_0.005.txt')
74: help
75: help()
76: #?
77: #?save
78: #?history
79: histogram
80: _ip.magic("history ")
81: _ip.magic("history -h")
82: #?history -
83: _ip.magic("history -f")
84: #?history
85: _ip.magic("history -f log 1")
86: _ip.system(" view log")
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
pure_sin.txt Pure sinewave, A = 2.9, fa = 100, phi = 0, 1MHz, t=0.02s
gaussian_100th.txt Same with 0.01 rms gaussian noise (SNR = 2e2)
gaussian_1000th.txt Same with 0.001 rms gaussian noise (SNR = 2e3)
incoherent.txt Pure sinewave as above, incoherently sampled (t = 0.025s)
incoherent_1000gaussian.txt Same as above with 0.001 rms gaussian noise (SNR = 2e3)
thd_44_7_times.txt Same as above with 0.01 and 0.02 harmonics (SNR = 44.7)
thd_0.1_0.08_n_0.005.txt With 0.1 and 0.08 harmonics and 0.005 gaussian noise
log How these were made
This diff is collapsed.
from sys import argv
def adj(x):
return round(x/5.0 * (2**15))
src = open(argv[1], "r")
lines = ['\t' + i for i in src.readlines()]
src.close()
dest = open(argv[2], "w")
dest.write("[SIGNAL]\nnbits = 24\nrate = 1000000\n=ndata = ")
dest.writelines(lines)
dest.close()
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
1 : import syn
2 : syn.synth(2.9, 100, 0.0, 1000000, 0, 0.02)
3 : s = +
4 : s = _
5 : s
6 : plot s
7 : from matplotlib import pyplot
8 : pyplot.plot(s)
9 : pyplot.show()
10: _ip.magic("save ")
11: help(save)
12: _ip.magic("save ")
13: save(s, 'pure_sin.txt')
14: save(s, 'pure_sin.txt')
15: from numpy import *
16: save(s, 'pure_sin.txt')
17: save('pure_sin.txt', s)
18: _ip.magic("edit pure_sin.txt.npy")
19: help(save)
20: out = file('pure_sin.txt', 'w')
21:
for i in s:
print >>out, i
22: out.close()
23: _ip.system("view pure*")
24: _ip.system("rm pure_sin.txt.npy")
25:
def w(s, f):
out = file(f, 'w')
for i in s:
print >>out, i
out.close()
26:
27:
28:
29: s
30: t
31: syn.synth(2.9, 100, 0.0, 1000000, 0, 0.02)
32: random.randn(len(s))
33:
34:
35: r = 2.9 * 0.001 * random.randn(len(s))
36: s+r
37: sn = s+r
38: w(s, 'gaussian_1000th.txt')
39: plot s
40: w(sn, 'gaussian_1000th.txt')
41: plot s
42: plot = pyplot.plot
43: plot(s)
44: sn = s+r
45: r = 2.9 * 0.01 * random.randn(len(s))
46: sn = s+r
47: w(sn, 'gaussian_100th.txt')
48: plot(sn)
49: r = 2.9 * 0.001 * random.randn(len(s))
50: sn = s+r
51: w(sn, 'gaussian_100th.txt')
52: plot(sn)
53: s = syn.synth(2.9, 100, 0.0, 1000000, 0, 0.025)
54: sincoh = s
55: w(sincoh, 'incoherent.txt')
56: sincohn = s + 0.001 * 2.9 * random.randn(len(s))
57: w(sincohn, 'incoherent_1000gaussian.txt')
58: plot(sincohn)
59: s = syn.synth(2.9, 100, 0.0, 1000000, 0, 0.025)
60: s = s + 0.01 * syn.synth(2.9, 300, 0.0, 1000000, 0, 0.025)
61: s = s + 0.02 * syn.synth(2.9, 500, 0.0, 1000000, 0, 0.025)
62: plot(s)
63: plot(s)
64: plot(s)
65: 1 / sqrt(5)
66: w(s, 'thd_44_7_times.txt')
67: s = syn.synth(2.9, 100, 0.0, 1000000, 0, 0.025)
68: s = s + 0.1 syn.synth(2.9, 500, 0.0, 1000000, 0, 0.025)
69: s = s + 0.1 * syn.synth(2.9, 500, 0.0, 1000000, 0, 0.025)
70: s = s + 0.08 * syn.synth(2.9, 800, 0.0, 1000000, 0, 0.025)
71: s = s + 0.005 * random.randn(len(s))
72: w(s,'thd_0.1_0.08_n_0.005.txt)
73: w(s,'thd_0.1_0.08_n_0.005.txt')
74: help
75: help()
76: #?
77: #?save
78: #?history
79: histogram
80: _ip.magic("history ")
81: _ip.magic("history -h")
82: #?history -
83: _ip.magic("history -f")
84: #?history
85: _ip.magic("history -f log 1")
86: _ip.system(" view log")
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment