add signal processing code and new GUI

Remove copyrighted docs and maciej program
parent e8a3e159
# 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 simple store the output
signal of an ADC with some other useful informations. """
from Utilities import *
from numpy import *
import Sinefit
class Signal(object):
"""A class that represent a time-domain sampled signal.
"""
def __init__(self, nbits = 0, rate = 1, 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.fulldata = data = array(data, dtype=float)
self.fullnsamples = len(data)
self.data = self.fulldata
self.nsamples = len(data)
# remove DC component
if self.fullnsamples > 0:
self.fulldata -= (max(self.fulldata) +min(self.fulldata))/2.
self.fulldft = abs(fft.fft(self.fulldata))
# useful names
N = len(data)
fdft = self.fulldft
first = 1. + argmax(fdft[1:N/2])
second = first + (argmax(fdft[first-1:first+2:2])*2) -1
ratio = (fdft[second] / fdft[first])
self.first = first
self.beta = N/pi * arctan(sin(pi/N)/(cos(pi/N)+1./ratio))
self.w0index = first+ self.beta
freqSample = 2 * pi * self.rate
w0 = freqSample * float(self.w0index)/self.nsamples
print "Frequency initial guess ->", w0
self.w0, self.A, self.B, self.C = Sinefit.sinefit4(data, 1.0/rate, w0, 1e-7)
print "Frequency fit ->", self.w0
# limit data
self.w0index = self.w0 /freqSample * self.nsamples
self.limit = floor(0.5 + N*int(self.w0index)/self.w0index)
self.data = data[:self.limit]
self.nsamples = len(self.data)
print "limit is:", self.limit
self.precalculateAll()
def items(self):
"""Create a list of tuples that describes the signal.
The structure of a tuple is"""
output = []
output.append(('Number of bits', '%d b', self.nbits))
output.append(('Sampling rate', '%.2f Hz', self.rate))
output.append(('Number of samples', "%d samples", self.nsamples))
output.append(('Peak at', "%f", self.w0index))
output.append(('Input frequency', "%.6f Hz", self.w0/2/pi))
output.append(('Beta', "%f", self.beta))
return output
def precalculateAll(self):
return
......@@ -31,16 +31,47 @@ def sinefit3(samples, sample_t, w0):
D0 = D0T.T
x0 = (D0T * D0).I * D0T * matrix(samples).T
return array(x0).reshape((3,))
def sinefit4matrix(samples, sample_t, w0, tol=1e-7):
"""fit a sampled wave to a sine of unknown frequency
def sinefit4(data, Ts, w0):
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
while True:
w0 += deltaw0
n = samples.size
t = arange(n) * sample_t
w0t = w0 * t
D0T = matrix(vstack([cos(w0t), sin(w0t), ones(n),
-a0*t*sin(w0t) + b0*t*cos(w0t)]))
D0 = D0T.T
x0 = (D0T * D0).I * D0T * matrix(samples).T
x0 = array(x0).reshape((4,))
a0, b0, c0, deltaw0 = x0
if abs(deltaw0/w0) < tol:
return (w0+deltaw0, a0, b0, c0)
def sinefit4scipy(data, Ts, w0, *args, **kwargs):
# 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: p[0]*cos(p[3]*x) +p[1]*sin(p[3]*x) +p[2] - y
errfunc = lambda p, x, y: y - (p[0]*cos(p[3]*x) +p[1]*sin(p[3]*x) +p[2])
# Initial guess for the parameters
A, B, C = sinefit3(data, Ts, w0)
A, B, C = sinefit3(data, Ts, w0)
p0 = [A, B, C, w0];
# Time serie
......@@ -49,7 +80,7 @@ def sinefit4(data, Ts, w0):
p1, success = optimize.leastsq(errfunc, p0[:], args=(Tx, data))
return p1
return p1[:4]
def makesine(samples, periods, bits, amplitude=1, noise=0):
t = arange(samples)/float(samples)
......@@ -60,6 +91,8 @@ def makesine(samples, periods, bits, amplitude=1, noise=0):
place(sine, sine <= -2**bits, -2**bits)
return sine
sinefit4 = sinefit4matrix
if __name__ == "__main__":
sine = makesine(20000, 20, 12, 443, 10)
print sine
......
from Signal import *
class TwoToneSignal(Signal):
"""A class that represent a time-domain sampled signal with two sinusoids
with similiar frequenciess.
"""
def __init__(self, nbits = 0, rate = 1, 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)
"""
super(TwoToneSignal, self).__init__(nbits, rate, data)
def items(self):
output = super(TwoToneSignal, self).items()
output.append(('The cake is a lie', "%s ", True))
return output
def precalculateAll(self):
"""Evaluates all the parameters of the signal, and also call the
precalculate method for each window function we know."""
return
from ConfigParser import RawConfigParser
from numpy import max, abs, sum, sqrt, log10
from time import time
# timeit decoratr
def timeit(method):
def timed(*args, **kw):
ts = time()
result = method(*args, **kw)
te = time()
print '%r (%r, %r) %2.3f sec' % (method.__name__, args, kw, te-ts)
return result
return timed
# db converter
def dB(x): return 20*log10(x)
def norm(v):
m = max(abs(v))
w = v/m
return m * sqrt(sum(w*w))
# empty generator
def fake(x, y):
return (i for i in [])
@timeit
def readFile(path):
config = RawConfigParser()
config.read(path)
nbits = config.getint('SIGNAL', 'nbits')
rate = config.getint('SIGNAL', 'rate')
dataString = config.get('SIGNAL', 'data').split('\n')
data = map(int, dataString)
return nbits, rate, data
......@@ -50,4 +50,4 @@ windows = {
if __name__ == "__main__":
for name in windows:
print "Using window", name
print windows[name][15]
\ No newline at end of file
print windows[name][15]
from numpy import *
class WindowedSignal:
data = array([])
th = array([])
nsamples = 0
dft = array([])
udft = array([])
ldft = array([])
ludft = array([])
ratio = 0
maxPeak = 0
w0 = 0
w0index = 0
rate = 0
def adjust(self, data, fs):
while data >= fs:
data -= fs
if data >= fs/2.:
data = fs -data;
return data
def harmonicPeaksGenerator(self, start, end):
indexes = (self.adjust(x * self.w0index, self.nsamples -1) for x in xrange(start , end))
return ((i, self.ratio*i, self.dft[i]) for i in indexes)
def logHarmonicPeaksGenerator(self, start, end):
indexes = (self.adjust(x * self.w0index, self.nsamples -1) for x in xrange(start , end))
return ((i, self.ratio*i, self.ldft[i]) for i in indexes)
THD = 0
noiseFloor = 0
signalPower = 0
noisePower = 0
SNR = 0
SINAD = 0
ENOB = 0
maxPeak = 0
secondPeak = 0
SFDR = 0
def items(self):
output = []
output.append(('Input frequency', '%.5f Hz', self.w0/(2*pi)))
output.append(('Peak', '%d', self.w0index))
output.append(('THD', '%.2f dB', self.THD))
output.append(('Noise floor', "%g dB", self.noiseFloor))
output.append(('Signal power', "%.f ", self.signalPower))
output.append(('Noise power', "%.f ", self.noisePower))
output.append(('SNR', "%.2f dB", self.SNR))
output.append(('SINAD', "%.2f dB", self.SINAD))
output.append(('ENOB', "%.2f b", self.ENOB))
output.append(('SFDR', "%.2f dBc", self.SFDR))
return output
def report(self):
for i in self.items():
print "%s: %s" % (i[0], i[1] % i[2])
__author__="federico"
__date__ ="$Jul 11, 2011 2:39:08 PM$"
import matplotlib
matplotlib.use('Qt4Agg')
# 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[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
__author__="federico"
__date__ ="$Jul 11, 2011 2:39:08 PM$"
\ No newline at end of file
from numpy import*
from math import sin, pi, atan
import matplotlib.pyplot as p
def sine_samples(n):