sinefit - directory for testing

e8c7e7a0 · Juan David González Cobas · d632e645 · e8c7e7a0 · e8c7e7a0 · e8c7e7a0
Commit e8c7e7a0 authored Jul 17, 2011 by Juan David González Cobas
30 changed files
--- a/sam/00readme.txt
+++ b/sam/00readme.txt
+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
--- a/sam/complex.txt
+++ b/sam/complex.txt
--- a/sam/data.txt
+++ b/sam/data.txt
--- a/sam/data_coherent.txt
+++ b/sam/data_coherent.txt
--- a/sam/diff-doc.txt
+++ b/sam/diff-doc.txt
--- a/sam/filter.py~
+++ b/sam/filter.py~
+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()    
+ 
+
--- a/sam/gaussian_1000th.txt
+++ b/sam/gaussian_1000th.txt
--- a/sam/gaussian_100th.txt
+++ b/sam/gaussian_100th.txt
--- a/sam/incoherent.txt
+++ b/sam/incoherent.txt
--- a/sam/incoherent_1000gaussian.txt
+++ b/sam/incoherent_1000gaussian.txt
--- a/sam/log
+++ b/sam/log
+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")
--- a/sam/pure_sin.txt
+++ b/sam/pure_sin.txt
--- a/sam/pure_sin.txt~
+++ b/sam/pure_sin.txt~
--- a/sam/simple.txt
+++ b/sam/simple.txt
+[SIGNAL]
+nbits = 16
+rate = 123
+data = 1
+  2
+  3
+  4
+  5
+
--- a/sam/sinefit.py
+++ b/sam/sinefit.py
+#!  /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,))
+
+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
+    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)
+
+
+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)
--- a/sam/snr.txt
+++ b/sam/snr.txt
+	                                                                     
+                                                                     
+                                                                     
+                                             
+-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
--- a/sam/test/00readme.txt
+++ b/sam/test/00readme.txt
+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
--- a/sam/test/diff-doc.txt
+++ b/sam/test/diff-doc.txt
--- a/sam/test/gaussian_1000th.txt
+++ b/sam/test/gaussian_1000th.txt
--- a/sam/test/gaussian_100th.txt
+++ b/sam/test/gaussian_100th.txt
--- a/sam/test/incoherent.txt
+++ b/sam/test/incoherent.txt
--- a/sam/test/incoherent_1000gaussian.txt
+++ b/sam/test/incoherent_1000gaussian.txt
--- a/sam/test/log
+++ b/sam/test/log
+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")
--- a/sam/test/pure_sin.txt
+++ b/sam/test/pure_sin.txt
--- a/sam/test/pure_sin.txt~
+++ b/sam/test/pure_sin.txt~
--- a/sam/test/thd_0.1_0.08_n_0.005.txt
+++ b/sam/test/thd_0.1_0.08_n_0.005.txt
--- a/sam/test/thd_44_7_times.txt
+++ b/sam/test/thd_44_7_times.txt
--- a/sam/thd_0.1_0.08_n_0.005.txt
+++ b/sam/thd_0.1_0.08_n_0.005.txt
--- a/sam/thd_44_7_times.txt
+++ b/sam/thd_44_7_times.txt
--- a/sam/x
+++ b/sam/x