Commit 17429fa0 authored by Federico Asara's avatar Federico Asara

Removed useless matplotlib import.

Signal class is as generic as possible.
SingleToneSignal reflects changes in Signal.
TwoTone signal computes IMD, but a mechanism to reduce incoherency is needed.
doubleSinefit4 fixed.
parent 4dd6d588
...@@ -38,55 +38,18 @@ class Signal(object): ...@@ -38,55 +38,18 @@ class Signal(object):
self.data = self.fulldata self.data = self.fulldata
self.nsamples = len(data) self.nsamples = len(data)
if self.fullnsamples > 0: if self.fullnsamples > 0:
# remove DC component # remove DC component
self.fulldata -= (max(self.fulldata) +min(self.fulldata))/2. self.fulldata -= (max(self.fulldata) +min(self.fulldata))/2.
self.data = self.fulldata
# calculate the |fft|
self.fulldft = abs(fft.fft(self.fulldata))
# useful names
N = len(data)
fdft = self.fulldft
# index of the biggest peak
first = 1. + argmax(fdft[1:N/2])
# index of the biggest peak nearest to `first`
# can only be first +-1.
second = first + (argmax(fdft[first-1:first+2:2])*2) -1
ratio = (fdft[second] / fdft[first])
# save first in self
self.first = first
# self.beta quantifies the sampling incoherency, defining the
# fraction of a period sampled in excess.
self.beta = N/pi * arctan(sin(pi/N)/(cos(pi/N)+1./ratio))
# the position the peak between first and second
self.w0index = first+ self.beta
# sampling frequency
freqSample = 2 * pi * self.rate
# initial frequency guess
w0 = freqSample * float(self.w0index)/self.nsamples
print "Frequency initial guess ->", w0
# fit the sine
self.w0, self.A, self.B, self.C = Sinefit.sinefit4(data, 1.0/rate, w0, 1e-7)
print "Frequency fit ->", self.w0
# limit data removing incoherency
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() self.precalculateAll()
def report(self):
for i in self.items():
print "%s: %s" % (i[0], i[1] % i[2])
def items(self): def items(self):
"""Create a list of tuples that describes the signal. """Create a list of tuples that describes the signal.
...@@ -101,9 +64,6 @@ class Signal(object): ...@@ -101,9 +64,6 @@ class Signal(object):
output.append(('Number of bits', '%d b', self.nbits)) output.append(('Number of bits', '%d b', self.nbits))
output.append(('Sampling rate', '%.2f Hz', self.rate)) output.append(('Sampling rate', '%.2f Hz', self.rate))
output.append(('Number of samples', "%d samples", self.nsamples)) 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 return output
......
...@@ -32,7 +32,7 @@ def doubleSinefit3(samples, sample_t, w0, w1): ...@@ -32,7 +32,7 @@ def doubleSinefit3(samples, sample_t, w0, w1):
x0 = (D0T * D0).I * D0T * matrix(samples).T x0 = (D0T * D0).I * D0T * matrix(samples).T
return array(x0).reshape((5,)) return array(x0).reshape((5,))
def doubleSinefit4matrix(samples, sample_t, w0, w1, tol=1e-7): def doubleSinefit4matrix(samples, sample_t, w0, w1, tol=1e-9):
"""fit a sampled wave to a sine of unknown frequency """fit a sampled wave to a sine of unknown frequency
This routine implements the algorithm of IEEE 1241, sect. 4.1.4.3., This routine implements the algorithm of IEEE 1241, sect. 4.1.4.3.,
...@@ -48,7 +48,7 @@ def doubleSinefit4matrix(samples, sample_t, w0, w1, tol=1e-7): ...@@ -48,7 +48,7 @@ def doubleSinefit4matrix(samples, sample_t, w0, w1, tol=1e-7):
The return value is the quadruple (w0, a0, b0, c0) The return value is the quadruple (w0, a0, b0, c0)
""" """
a0, b0, a1, b1, c0 = sinefit3(samples, sample_t, w0, w1) a0, b0, a1, b1, c0 = doubleSinefit3(samples, sample_t, w0, w1)
deltaw0 = 0 deltaw0 = 0
deltaw1 = 0 deltaw1 = 0
......
...@@ -50,6 +50,9 @@ class SingleToneSignal(Signal): ...@@ -50,6 +50,9 @@ class SingleToneSignal(Signal):
def items(self): def items(self):
output = super(SingleToneSignal, self).items() output = super(SingleToneSignal, self).items()
output.append(('Peak at', "%f", self.w0index))
output.append(('Input frequency', "%.6f Hz", self.w0/2/pi))
output.append(('Beta', "%f", self.beta))
output.append(('Max DNL', "%f ", self.maxDNL)) output.append(('Max DNL', "%f ", self.maxDNL))
output.append(('Max INL', "%f ", self.maxINL)) output.append(('Max INL', "%f ", self.maxINL))
output.append(('Theoretical SNR', "%.2f dB", self.thSNR)) output.append(('Theoretical SNR', "%.2f dB", self.thSNR))
...@@ -61,7 +64,56 @@ class SingleToneSignal(Signal): ...@@ -61,7 +64,56 @@ class SingleToneSignal(Signal):
def precalculateAll(self): def precalculateAll(self):
"""Evaluates all the parameters of the signal, and also call the """Evaluates all the parameters of the signal, and also call the
precalculate method for each window function we know.""" precalculate method for each window function we know."""
if self.fullnsamples > 0:
# remove DC component
# self.fulldata -= (max(self.fulldata) +min(self.fulldata))/2.
# calculate the |fft|
self.fulldft = abs(fft.fft(self.fulldata))
# useful names
data = self.data
rate = self.rate
N = len(data)
fdft = self.fulldft
# index of the biggest peak
first = 1. + argmax(fdft[1:N/2])
# index of the biggest peak nearest to `first`
# can only be first +-1.
second = first + (argmax(fdft[first-1:first+2:2])*2) -1
ratio = (fdft[second] / fdft[first])
# save first in self
self.first = first
# self.beta quantifies the sampling incoherency, defining the
# fraction of a period sampled in excess.
self.beta = N/pi * arctan(sin(pi/N)/(cos(pi/N)+1./ratio))
# the position the peak between first and second
self.w0index = first+ self.beta
# sampling frequency
freqSample = 2 * pi * self.rate
# initial frequency guess
w0 = freqSample * float(self.w0index)/self.nsamples
print "Frequency initial guess ->", w0
# fit the sine
self.w0, self.A, self.B, self.C = Sinefit.sinefit4(data, 1.0/rate, w0, 1e-7)
print "Frequency fit ->", self.w0
# limit data removing incoherency
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
# First of all evaluate the histograms # First of all evaluate the histograms
skip = False skip = False
......
...@@ -4,6 +4,10 @@ class TwoToneSignal(Signal): ...@@ -4,6 +4,10 @@ class TwoToneSignal(Signal):
"""A class that represent a time-domain sampled signal with two sinusoids """A class that represent a time-domain sampled signal with two sinusoids
with similiar frequenciess. with similiar frequenciess.
""" """
m1 = 0
m2 = 0
tow1 = 0.
tow2 = 0.
def __init__(self, nbits = 0, rate = 1, data = []): def __init__(self, nbits = 0, rate = 1, data = []):
"""initialize a signal object """initialize a signal object
...@@ -19,11 +23,99 @@ class TwoToneSignal(Signal): ...@@ -19,11 +23,99 @@ class TwoToneSignal(Signal):
def items(self): def items(self):
output = super(TwoToneSignal, self).items() output = super(TwoToneSignal, self).items()
output.append(('The cake is a lie', "%s ", True)) output.append(('Peak 1 at', "%d bin", self.m1))
output.append(('Peak 2 at', "%d bin", self.m2))
output.append(('Guess for wave 1', "%f Hz", self.tow1/2./pi))
output.append(('Guess for wave 2', "%f Hz", self.tow2/2./pi))
output.append(('Sine 1 frequency', "%f Hz", self.w1/2./pi))
output.append(('Sine 1 A', "%f", self.a1))
output.append(('Sine 1 B', "%f", self.b1))
output.append(('Sine 2 frequency', "%f Hz", self.w2/2./pi))
output.append(('Sine 2 A', "%f", self.a2))
output.append(('Sine 2 B', "%f", self.b2))
output.append(('DC component', "%f", self.c0))
output.append(('IMD', "%f dBc", self.imd))
return output return output
def toi(self, x):
return self.nsamples*x/2/pi/self.rate
def precalculateAll(self): def precalculateAll(self):
"""Evaluates all the parameters of the signal, and also call the """Evaluates all the parameters of the signal, and also call the
precalculate method for each window function we know.""" precalculate method for each window function we know."""
# ok, compute the FFT
ff = abs(fft.fft(self.fulldata))
lbnd = max(ff) * 10e-12
self.fft = ff = where(ff < lbnd, 10e-12, ff)
hff = ff[:len(ff)/2]
self.lfft = lff = 10*log10(ff)
self.m1 = m1 = argmax(hff);
self.m2 = m2 = argmax(hstack([hff[:m1-1], array([0, 0, 0]), hff[m1+2:]]))
self.tow1 = tow1 = 2*pi*self.rate*float(m1)/self.nsamples
self.tow2 = tow2 = 2*pi*self.rate*float(m2)/self.nsamples
(self.w1, self.a1, self.b1), (self.w2, self.a2, self.b2), self.c0 = \
Sinefit.doubleSinefit4matrix(self.fulldata, self.rate**-1, tow1, tow2)
delta = abs(self.w1 - self.w2)
i1 = min([self.w1, self.w2]) - delta
i2 = max([self.w1, self.w2]) + delta
fw1, fw2 = self.fft[self.toi(self.w1)], self.fft[self.toi(self.w2)]
fi1, fi2 = max(self.fft[self.toi(i1)-1:self.toi(i1)+2]), max(self.fft[self.toi(i2)-1:self.toi(i2)+2])
print fw1, fw2
print fi1, fi2
meaningful = hypot(fw1, fw2)
interferences = hypot(fi1, fi2)
self.imd = 10*log10(meaningful/interferences)
# fix incoherency
self.report()
fd = abs(self.w1 + self.w2)/4./pi
print "fd is", fd
wModIndex = self.nsamples * fd / self.rate
print "wModIndex is", wModIndex
factor = (wModIndex -(wModIndex % 2))/wModIndex
limit = floor(self.nsamples*factor)
print "limit is", limit
return
self.data = self.fulldata[:limit]
self.nsamples = limit
ff = abs(fft.fft(self.data))
lbnd = max(ff) * 10e-12
self.fft = ff = where(ff < lbnd, 10e-12, ff)
hff = ff[:len(ff)/2]
self.lfft = lff = 10*log10(ff)
self.m1 = m1 = argmax(hff);
self.m2 = m2 = argmax(hstack([hff[:m1-1], array([0, 0, 0]), hff[m1+2:]]))
self.tow1 = tow1 = 2*pi*self.rate*float(m1)/self.nsamples
self.tow2 = tow2 = 2*pi*self.rate*float(m2)/self.nsamples
(self.w1, self.a1, self.b1), (self.w2, self.a2, self.b2), self.c0 = \
Sinefit.doubleSinefit4matrix(self.data, self.rate**-1, tow1, tow2)
fd = abs(self.w1 - self.w2)/2./pi
print "fd is", fd
return return
...@@ -2,6 +2,3 @@ __author__="federico" ...@@ -2,6 +2,3 @@ __author__="federico"
__date__ ="$Jul 11, 2011 2:39:08 PM$" __date__ ="$Jul 11, 2011 2:39:08 PM$"
import matplotlib
matplotlib.use('Qt4Agg')
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