Skip to content
Snippets Groups Projects
Commit a0e40b1f authored by Thomas Bruns's avatar Thomas Bruns
Browse files

corrected 4-param fit, scipy -> numpy migration

parent 54bf8067
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 256 additions and 316 deletions
SineTools.py
+ 256
− 316
View file @ a0e40b1f
......@@ -8,22 +8,22 @@ auxiliary functions related to sine-approximation
"""
import scipy as sp
import numpy as np
from scipy import linalg as la
import matplotlib.pyplot as mp
def sampletimes(Fs, T): #
def sampletimes(Fs, T): #
"""
generate a t_i vector with \n
sample rate Fs \n
from 0 to T
"""
num = sp.ceil(T * Fs)
return sp.linspace(0, T, num, dtype=sp.float64)
num = np.ceil(T * Fs)
return np.linspace(0, T, num, dtype=np.float64)
# a = displacement, velocity or acceleration amplitude
def sinewave(f, a, phi, ti, offset=0, noise=0, absnoise=0, drift=0, ampdrift=0):
def sinewave(f,a,phi, ti, offset=0, noise=0, absnoise=0, drift=0, ampdrift=0):
"""
generate a sampled sine wave s_i = s(t_i) with \n
amplitude a \n
......@@ -38,25 +38,18 @@ def sinewave(f, a, phi, ti, offset=0, noise=0, absnoise=0, drift=0, ampdrift=0):
Tau = ti[-1] - ti[0]
n = 0
n = 0
if noise != 0:
n = a * noise * sp.randn(len(ti))
if noise != 0:
n = a*noise*sp.randn(len(ti))
if absnoise != 0:
n = n + absnoise * sp.randn(len(ti))
d = drift * a / Tau
s = (
a * (1 + ampdrift / Tau * ti) * sp.sin(2 * sp.pi * f * ti - phi)
+ n
+ d * ti
+ offset
)
n = n + absnoise*sp.randn(len(ti))
d = drift*a/Tau
s = a*(1+ampdrift/Tau*ti)*np.sin(2*np.pi*f * ti - phi) + n +d*ti + offset
return s
def fm_counter_sine(
fm, f, x, phi, ti, offset=0, noise=0, absnoise=0, drift=0, ampdrift=0, lamb=633.0e-9
):
def fm_counter_sine(fm, f, x, phi, ti, offset=0, noise=0, absnoise=0, drift=0, ampdrift=0,lamb=633.0e-9):
"""
calculate counter value of heterodyne signal at \n
carrier freq. fm\n
......@@ -72,30 +65,20 @@ def fm_counter_sine(
"""
Tau = ti[-1] - ti[0]
n = 0
if noise != 0:
n = x * noise * sp.randn(len(ti))
if noise != 0:
n = x*noise*sp.randn(len(ti))
if absnoise != 0:
n = n + absnoise * sp.randn(len(ti))
d = drift * x / Tau
s = (
1.0
/ lamb
* (
x * (1 + ampdrift / Tau * ti) * sp.sin(2 * sp.pi * f * ti - phi)
+ n
+ d * ti
+ offset
)
)
s = sp.floor(s + fm * ti)
n = n + absnoise*sp.randn(len(ti))
d = drift*x/Tau
s = 1.0/lamb * (x*(1+ampdrift/Tau*ti)*np.sin(2*np.pi*f * ti - phi) + n +d*ti + offset)
s = np.floor(s+fm*ti)
return s
# sine fit at known frequency
def threeparsinefit(y, t, f0):
def threeparsinefit(y,t,f0):
"""
sine-fit at a known frequency\n
y vector of sample values \n
......@@ -105,16 +88,15 @@ def threeparsinefit(y, t, f0):
returns a vector of coefficients [a,b,c]\n
for y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c
"""
w0 = 2 * sp.pi * f0
a = sp.array([sp.cos(w0 * t), sp.sin(w0 * t), sp.ones(t.size)])
abc = la.lstsq(a.transpose(), y)
return abc[0][0:3] ## fit vector a*sin+b*cos+c
w0 = 2*np.pi*f0
a = np.array([np.cos(w0*t),np.sin(w0*t),np.ones(t.size)])
abc = la.lstsq(a.transpose(),y)
return abc[0][0:3] ## fit vector a*sin+b*cos+c
# sine fit at known frequency and detrending
def threeparsinefit_lin(y, t, f0):
def threeparsinefit_lin(y,t,f0):
"""
sine-fit with detrending at a known frequency\n
y vector of sample values \n
......@@ -124,20 +106,20 @@ def threeparsinefit_lin(y, t, f0):
returns a vector of coefficients [a,b,c,d]\n
for y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c*t + d
"""
w0 = 2 * sp.pi * f0
a = sp.array([sp.cos(w0 * t), sp.sin(w0 * t), sp.ones(t.size), t, sp.ones(t.size)])
abc = la.lstsq(a.transpose(), y)
return abc[0][0:4] ## fit vector
w0 = 2*np.pi*f0
a = np.array([np.cos(w0*t),np.sin(w0*t),np.ones(t.size),t,np.ones(t.size)])
abc = la.lstsq(a.transpose(),y)
return abc[0][0:4] ## fit vector
def calc_threeparsine(abc, t, f0):
def calc_threeparsine(abc,t,f0):
"""
return y = abc[0]*sin(2*pi*f0*t) + abc[1]*cos(2*pi*f0*t) + abc[2]
"""
w0 = 2 * sp.pi * f0
return abc[0] * sp.cos(w0 * t) + abc[1] * sp.sin(w0 * t) + abc[2]
w0 =2*np.pi*f0
return abc[0]*np.cos(w0*t)+abc[1]*np.sin(w0*t)+abc[2]
def amplitude(abc):
......@@ -145,7 +127,7 @@ def amplitude(abc):
return the amplitude given the coefficients of\n
y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c
"""
return sp.absolute(abc[1] + 1j * abc[0])
return np.absolute(abc[1]+1j*abc[0])
def phase(abc, deg=False):
......@@ -154,10 +136,9 @@ def phase(abc, deg=False):
y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c \n
returns angle in rad by default, in degree if deg=True
"""
return sp.angle(abc[1] + 1j * abc[0], deg=deg)
def magnitude(A1, A2):
return np.angle(abc[1]+1j*abc[0], deg=deg)
def magnitude(A1,A2):
"""
return the magnitude of the complex ratio of sines A2/A1\n
given two sets of coefficients \n
......@@ -166,8 +147,7 @@ def magnitude(A1, A2):
"""
return amplitude(A2) / amplitude(A1)
def phase_delay(A1, A2, deg=False):
def phase_delay(A1,A2, deg=False):
"""
return the phase difference of the complex ratio of sines A2/A1\n
given two sets of coefficients \n
......@@ -177,9 +157,8 @@ def phase_delay(A1, A2, deg=False):
"""
return phase(A2, deg=deg) - phase(A1, deg=deg)
# periodical sinefit at known frequency
def seq_threeparsinefit(y, t, f0):
# periodical sinefit at known frequency
def seq_threeparsinefit(y,t,f0):
"""
period-wise sine-fit at a known frequency\n
y vector of sample values \n
......@@ -189,22 +168,21 @@ def seq_threeparsinefit(y, t, f0):
returns a (n,3)-matrix of coefficient-triplets [[a,b,c], ...]\n
for y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c
"""
Tau = 1.0 / f0
dt = t[1] - t[0]
N = int(Tau / dt) ## samples per section
M = int(sp.floor(t.size / N)) ## number of sections or periods
abc = sp.zeros((M, 3))
Tau = 1.0/f0
dt = t[1]-t[0]
N = int(Tau/dt) ## samples per section
M = int(np.floor(t.size/N)) ## number of sections or periods
abc = np.zeros((M,3))
for i in range(int(M)):
ti = t[i * N : (i + 1) * N]
yi = y[i * N : (i + 1) * N]
abc[i, :] = threeparsinefit(yi, ti, f0)
return abc ## matrix of all fit vectors per period
ti = t[i*N:(i+1)*N]
yi = y[i*N:(i+1)*N]
abc[i,:] = (threeparsinefit(yi,ti,f0))
return abc ## matrix of all fit vectors per period
# four parameter sine-fit (with frequency approximation)
def fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000):
def fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000):
"""
y sampled data values \n
t sample times of y \n
......@@ -214,39 +192,34 @@ def fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000):
\n
returns the vector [a, b, c, w] of a*sin(w*t)+b*cos(w*t)+c
"""
abcd = threeparsinefit(y, t, f0)
w = 2 * sp.pi * f0
abcd = threeparsinefit(y,t,f0)
w = 2*np.pi*f0
err = 1
i = 0
while (err > tol) and (i < nmax):
D = sp.array(
[
sp.cos(w * t),
sp.sin(w * t),
sp.ones(t.size),
(-1.0) * abcd[0] * t * sp.sin(w * t) + abcd[1] * t * sp.cos(w * t),
]
)
abcd = (la.lstsq(D.transpose(), y))[0]
while ((err > tol) and (i<nmax)):
D = np.array([np.cos(w*t),
np.sin(w*t),
np.ones(t.size),
(-1.0)*abcd[0]*t*np.sin(w*t)+abcd[1]*t*np.cos(w*t)
])
print(dw)
abcd = (la.lstsq(D.transpose(),y))[0]
dw = abcd[3]
w = w + 0.9 * dw
i += 1
err = sp.absolute(dw / w)
w = w+0.9*dw
i+=1
err = np.absolute(dw/w)
assert i < nmax, "iteration error"
return sp.hstack((abcd[0:3], w / (2 * sp.pi)))
def calc_fourparsine(abcf, t):
return np.hstack((abcd[0:3],w/(2*np.pi)))
def calc_fourparsine(abcf,t):
"""
return y = abc[0]*sin(2*pi*f0*t) + abc[1]*cos(2*pi*f0*t) + abc[2]
"""
w0 = 2 * sp.pi * abcf[3]
return abcf[0] * sp.cos(w0 * t) + abcf[1] * sp.sin(w0 * t) + abcf[2]
w0 =2*np.pi*abcf[3]
return abcf[0]*np.cos(w0*t)+abcf[1]*np.sin(w0*t)+abcf[2]
"""
from octave ...
function abcw = fourParSinefit(data,w0)
......@@ -271,8 +244,8 @@ function abcw = fourParSinefit(data,w0)
endfunction
"""
# periodical sinefit at known frequency
def seq_fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000, debug_plot=False):
# periodical sinefit at known frequency
def seq_fourparsinefit(y,t,f0, tol=1.0e-7, nmax=1000, debug_plot=False):
"""
period-wise sine-fit at a known frequency\n
y vector of sample values \n
......@@ -284,41 +257,39 @@ def seq_fourparsinefit(y, t, f0, tol=1.0e-7, nmax=1000, debug_plot=False):
returns a (n,3)-matrix of coefficient-triplets [[a,b,c], ...]\n
for y = a*sin(2*pi*f0*t) + b*cos(2*pi*f0*t) + c
"""
Tau = 1.0 / f0
dt = t[1] - t[0]
N = int(sp.floor(Tau / dt)) ## samples per section
M = int(sp.floor(t.size / N)) ## number of sections or periods
abcd = sp.zeros((M, 4))
Tau = 1.0/f0
dt = t[1]-t[0]
N = int(np.floor(Tau/dt)) ## samples per section
M = int(np.floor(t.size/N)) ## number of sections or periods
abcd = np.zeros((M,4))
for i in range(M):
ti = t[i * N : (i + 1) * N]
yi = y[i * N : (i + 1) * N]
abcd[i, :] = fourparsinefit(yi, ti, f0, tol=tol, nmax=nmax)
ti = t[i*N:(i+1)*N]
yi = y[i*N:(i+1)*N]
abcd[i,:] = (fourparsinefit(yi,ti,f0,tol=tol,nmax=nmax))
if debug_plot:
mp.ioff()
fig = mp.figure("seq_fourparsinefit")
fig.clear()
p1 = fig.add_subplot(211)
p2 = fig.add_subplot(212, sharex=p1)
p2 = fig.add_subplot(212,sharex=p1)
for i in range(M):
p1.plot(t[i * N : (i + 1) * N], y[i * N : (i + 1) * N], ".")
s = calc_fourparsine(
abcd[i, :], t[i * N : (i + 1) * N]
) # fitted data to plot
p1.plot(t[i * N : (i + 1) * N], s, "-")
r = y[i * N : (i + 1) * N] - s # residuals to plot
p2.plot(t[i * N : (i + 1) * N], r, ".")
yi = y[i * N : (i + 1) * N]
p1.plot(t[i*N:(i+1)*N],y[i*N:(i+1)*N],".")
s = calc_fourparsine(abcd[i,:],t[i*N:(i+1)*N]) # fitted data to plot
p1.plot(t[i*N:(i+1)*N],s,"-")
r = y[i*N:(i+1)*N]-s # residuals to plot
p2.plot(t[i*N:(i+1)*N],r,".")
yi = y[i*N:(i+1)*N]
mp.show()
return abcd ## matrix of all fit vectors per period
return abcd ## matrix of all fit vectors per period
# fitting a pseudo-random multi-sine signal with 2*Nf+1 parameters
def multi_threeparsinefit(y, t, f0): # fo vector of frequencies
def multi_threeparsinefit(y,t,f0): # f0 vector of frequencies
"""
fit a time series of a sum of sine-waveforms with a given set of frequencies\n
y vector of sample values \n
......@@ -328,37 +299,37 @@ def multi_threeparsinefit(y, t, f0): # fo vector of frequencies
returns a vector of coefficient-triplets [a,b,c] for the frequencies\n
for y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + c_i
"""
w0 = 2 * sp.pi * f0
w0 = 2 * np.pi * f0
D = np.ones((len(t),1)) # for the bias
# set up design matrix
a = sp.ones(len(t))
for w in w0:
a = sp.vstack((sp.vstack((sp.cos(w * t), sp.sin(w * t))), a))
for w in w0[::-1]:
D = np.hstack(( np.cos(w*t)[:,np.newaxis],
np.sin(w*t)[:,np.newaxis],
D))
abc = sp.linalg.lstsq(a.transpose(), y)
return abc[0] ## fit vector a*sin+b*cos+c
abc = np.linalg.lstsq(D, y)
return abc[0] ## fit vector a*cos+b*sin+c
def multi_amplitude(abc): # abc = [a1,b1 , a2,b2, ...,bias]
def multi_amplitude(abc): # abc = [a1,b1 , a2,b2, ...,bias]
"""
return the amplitudes given the coefficients of a multi-sine\n
abc = [a1,b1 , a2,b2, ...,bias] \n
y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + c_i
"""
x = abc[0::2] + 1j * abc[1::2]
return sp.absolute(x)
x = abc[:-1][0::2] + 1j*abc[:-1][1::2] # make complex without Bias (last value)
return np.absolute(x)
def multi_phase(abc, deg=False): # abc = [bias, a1,b1 , a2,b2, ...]
def multi_phase(abc, deg=False): # abc = [bias, a1,b1 , a2,b2, ...]
"""
return the initial phases given the coefficients of a multi-sine\n
abc = [a1,b1 , a2,b2, ...,bias] \n
y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + c_i
"""
x = abc[1::2] + 1j * abc[2::2]
return sp.angle(x, deg=deg)
x = abc[:-1][0::2] + 1j*abc[:-1][1::2] # make complex without Bias (last value)
return np.angle(x,deg=deg)
def multi_waveform_abc(f, abc, t):
def multi_waveform_abc(f,abc,t):
"""
generate a sample time series of a multi-sine from coefficients and frequencies\n
f vector of given frequencies \n
......@@ -368,18 +339,37 @@ def multi_waveform_abc(f, abc, t):
returns the vector \n
y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + bias
"""
ret = 0.0 * t + abc[-1] # bias
for fi, a, b in zip(f, abc[0::2], abc[1::2]):
ret = ret + a * sp.cos(2 * sp.pi * fi * t) + b * sp.sin(2 * sp.pi * fi * t)
return ret
ret = 0.0*t +abc[-1] # bias
for fi,a,b in zip(f,abc[0::2],abc[1::2]):
ret = ret+a*np.cos(2*np.pi*fi*t)+b*np.sin(2*np.pi*fi*t)
return ret
def multi_waveform_mp(f,m,p,t, bias=0, deg=True):
"""
generate a sample time series of a multi-sine from magnitude/phase and frequencies\n
f vector of given frequencies \n
m vector of magnitudes \n
p vector of phases \n
t vector of sample times t_i \n
bias scalar value for total bias
deg = boolean whether phase is in degree
\n
returns the vector \n
y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + bias
"""
ret = 0.0*t + bias # init and bias
if deg:
p = np.deg2rad(p)
for fi,m_i,p_i in zip(f,m,p):
ret = ret + m_i*np.sin(2*np.pi*fi*t + p_i)
return ret
##################################
# Counter based stuff
# periodical sinefit to the linearly increasing heterodyne counter
# version based on Blume
def seq_threeparcounterfit(y, t, f0, diff=False):
def seq_threeparcounterfit(y,t,f0, diff=False):
"""
period-wise (single-)sinefit to the linearly increasing heterodyne counter
version based on "Blume et al. "\n
......@@ -391,22 +381,21 @@ def seq_threeparcounterfit(y, t, f0, diff=False):
if diff=True use differentiation to remove carrier (c.f. source)
"""
Tau = 1.0 / f0
dt = t[1] - t[0]
N = int(sp.floor(Tau / dt)) ## samples per section
M = int(sp.floor(t.size / N)) ## number of sections or periods
Tau = 1.0/f0
dt = t[1]-t[0]
N = int(np.floor(Tau/dt)) ## samples per section
M = int(np.floor(t.size/N)) ## number of sections or periods
remove_counter_carrier(y, diff=diff)
abc = sp.zeros((M, 4))
abc = np.zeros((M,4))
for i in range(int(M)):
ti = t[i * N : (i + 1) * N]
yi = y[i * N : (i + 1) * N]
abc[i, :] = threeparsinefit_lin(yi, ti, f0)
return abc ## matrix of all fit vectors per period
ti = t[i*N:(i+1)*N]
yi = y[i*N:(i+1)*N]
abc[i,:] = threeparsinefit_lin(yi,ti,f0)
return abc ## matrix of all fit vectors per period
def remove_counter_carrier(y, diff=False):
"""
......@@ -414,21 +403,18 @@ def remove_counter_carrier(y, diff=False):
generated by the carrier frequency of a heterodyne signal\n
y vector of samples of the signal
"""
if diff:
d = sp.diff(y)
d = d - sp.mean(d)
y = sp.hstack((0, sp.cumsum(d)))
if diff :
d = np.diff(y)
d = d - np.mean(d)
y = np.hstack((0,np.cumsum(d)))
else:
slope = y[-1] - y[0] # slope of linear increment
y = y - slope * sp.linspace(
0.0, 1.0, len(y), endpoint=False
) # removal of linear increment
slope = (y[-1]-y[0]) # slope of linear increment
y = y-slope*np.linspace(0.0,1.0,len(y),endpoint=False) # removal of linear increment
return y
# calculate displacement and acceleration to the same analytical s(t)
# Bsp: fm = 2e7, f=10, s0=0.15, phi0=sp.pi/3, ti, drift=0.03, ampdrift=0.03,thd=[0,0.02,0,0.004]
def disp_acc_distorted(fm, f, s0, phi0, ti, drift=0, ampdrift=0, thd=0):
# Bsp: fm = 2e7, f=10, s0=0.15, phi0=np.pi/3, ti, drift=0.03, ampdrift=0.03,thd=[0,0.02,0,0.004]
def disp_acc_distorted(fm, f,s0,phi0, ti, drift=0, ampdrift=0, thd=0):
"""
calculate the respective (displacement-) counter and acceleration
for a parmeterized distorted sine-wave motion in order to compare accelerometry with interferometry \n
......@@ -440,59 +426,45 @@ def disp_acc_distorted(fm, f, s0, phi0, ti, drift=0, ampdrift=0, thd=0):
ampdrift is displacement amplitude druft\n
thd is vector of higher harmonic amplitudes (c.f. source)
"""
om = 2 * sp.pi * f
om2 = om ** 2
tau = ti[-1] - ti[0]
disp = sp.sin(om * ti + phi0)
if thd != 0:
om = 2*np.pi*f
om2 = om**2
tau = ti[-1]-ti[0]
disp = np.sin(om*ti+phi0)
if (thd != 0) :
i = 2
for h in thd:
disp = disp + h * sp.sin(i * om * ti + phi0)
i = i + 1
if ampdrift != 0:
disp = disp * (1 + ampdrift / tau * ti)
for h in thd :
disp = disp + h*np.sin(i*om*ti+phi0)
i = i+1
if ampdrift != 0 :
disp = disp *(1+ampdrift/tau*ti)
if drift != 0:
disp = disp + s0 * drift / tau * ti
disp = disp * s0
disp = sp.floor((disp * 2 / 633e-9) + fm * ti)
disp = disp + s0*drift/tau*ti
disp = disp*s0
disp = np.floor((disp*2/633e-9)+fm*ti)
acc = -s0 * om2 * (1 + ampdrift / tau * ti) * sp.sin(om * ti + phi0)
acc = -s0*om2*(1+ampdrift/tau*ti)*np.sin(om*ti+phi0)
if ampdrift != 0:
acc = acc + (2 * ampdrift * s0 * om * sp.cos(om * ti + phi0)) / tau
acc = acc+ (2*ampdrift*s0*om*np.cos(om*ti+phi0))/tau
if thd != 0:
i = 2
for h in thd:
acc = acc - s0 * h * om2 * (1 + ampdrift / tau * ti) * i ** 2 * sp.sin(
i * om * ti + phi0
)
for h in thd :
acc = acc - s0*h*om2*(1+ampdrift/tau*ti)*i**2*np.sin(i*om*ti+phi0)
if ampdrift != 0:
acc = (
acc
+ (2 * ampdrift * s0 * om * i * h * sp.cos(om * ti + phi0)) / tau
)
i = i + 1
acc = acc + (2*ampdrift*s0*om*i*h*np.cos(om*ti+phi0))/tau
i=i+1
return disp,acc
return disp, acc
###################################
# Generation and adaptation of Parameters of the Multi-Sine considering hardware constraints
def PR_MultiSine_adapt(
f1,
Nperiods,
Nsamples,
Nf=8,
fs_min=0,
fs_max=1e9,
frange=10,
log=True,
phases=None,
sample_inkr=1,
):
def PR_MultiSine_adapt(f1, Nperiods, Nsamples, Nf=8, fs_min=0,fs_max=1e9, frange=10, log=True, phases=None, sample_inkr=1):
"""
Returns an additive normalized Multisine time series. \n
f1 = start frequency (may be adapted) \n
Nperiods = number of periods of f1 (may be increased) \n
f1 = start frequency (may be adapted by the algorithm) \n
Nperiods = number of periods of f1 (may be increased by algorithm) \n
Nsamples = Minimum Number of samples \n
Nf = number of frequencies in multi frequency mix \n
fs_min = minimum sample rate of used device (default 0) \n
......@@ -506,103 +478,79 @@ def PR_MultiSine_adapt(
returns: freq,phase,fs,ti,multi \n
freq= array of frequencies \n
phase=used phases in deg or rad \n
T1 = (adapted) duration
fs=sample rate \n
ti=timestamps \n
multi=array of time series values \n
"""
if (
Nsamples // sample_inkr * sample_inkr != Nsamples
): # check multiplicity of sample_inkr
Nsamples = (
Nsamples // sample_inkr + 1
) * sample_inkr # round to next higher multiple
if Nsamples//sample_inkr*sample_inkr != Nsamples: # check multiplicity of sample_inkr
Nsamples = (Nsamples//sample_inkr +1)*sample_inkr # round to next higher multiple
T0 = Nperiods / f1 # given duration
fs0 = Nsamples / T0 # (implicitly) given sample rate
T0 = Nperiods/f1 # given duration
fs0 = Nsamples/T0 # (implicitly) given sample rate
if False:
print("0 Nperiods: " + str(Nperiods))
print("0 Nsamples: " + str(Nsamples))
print("0 fs: " + str(fs0))
print("0 T0: " + str(T0))
print("0 f1: " + str(f1))
print ("0 Nperiods: "+str(Nperiods))
print ("0 Nsamples: "+str(Nsamples))
print ("0 fs: "+str(fs0))
print ("0 T0: "+str(T0))
print ("0 f1: "+str(f1))
fs = fs0
if fs0 < fs_min: # sample rate too low, then set to minimum
if fs0 < fs_min: # sample rate too low, then set to minimum
fs = fs_min
print("sample rate increased")
elif fs0 > fs_max: # sample rate too high, set to max-allowed and
elif fs0 > fs_max: # sample rate too high, set to max-allowed and
fs = fs_max
Nperiods = sp.ceil(
Nperiods * fs0 / fs_max
) # increase number of periods to get at least Nsamples samples
T0 = Nperiods / f1
print("sample rate reduced, Nperiods=" + str(Nperiods))
Nsamples = T0 * fs
if (
Nsamples // sample_inkr * sample_inkr != Nsamples
): # check multiplicity of sample_inkr
Nsamples = (
Nsamples // sample_inkr + 1
) * sample_inkr # round to next higher multiple
T1 = Nsamples / fs # adapt exact duration
f1 = Nperiods / T1 # adapt f1 for complete cycles
Nperiods = np.ceil(Nperiods*fs0/fs_max) # increase number of periods to get at least Nsamples samples
T0 = Nperiods/f1
print("sample rate reduced, Nperiods="+str(Nperiods))
Nsamples = T0*fs
if Nsamples//sample_inkr*sample_inkr != Nsamples: # check multiplicity of sample_inkr
Nsamples = (Nsamples//sample_inkr +1)*sample_inkr # round to next higher multiple
T1 = Nsamples/fs # adapt exact duration
f1 = Nperiods/T1 # adapt f1 for complete cycles
if False:
print("Nperiods: " + str(Nperiods))
print("Nsamples: " + str(Nsamples))
print("fs: " + str(fs))
print("T1: " + str(T1))
print("f1: " + str(f1))
print ("Nperiods: "+str(Nperiods))
print ("Nsamples: "+str(Nsamples))
print ("fs: "+str(fs))
print ("T1: "+str(T1))
print ("f1: "+str(f1))
f_res = 1 / T1 # frequency resolution
f_res = 1/T1 # frequency resolution
# determine a series of frequencies (freq[])
if log:
fact = sp.power(frange, 1.0 / (Nf - 1)) # factor for logarithmic scale
freq = f1 * sp.power(fact, sp.arange(Nf))
else:
step = (frange - 1) * f1 / (Nf - 1)
freq = sp.arange(f1, frange * f1 + step, step)
fact = np.power(frange, 1.0/(Nf-1)) # factor for logarithmic scale
freq = f1*np.power(fact, np.arange(Nf))
else:
step = (frange-1)*f1/(Nf-1)
freq = np.arange(f1,frange*f1+step,step)
# auxiliary function to find the nearest available frequency
def find_nearest(
x, possible
): # match the theoretical freqs to the possible periodic freqs
idx = (sp.absolute(possible - x)).argmin()
def find_nearest(x,possible): # match the theoretical freqs to the possible periodic freqs
idx = (np.absolute(possible - x)).argmin()
return possible[idx]
fi_pos = sp.arange(f1, frange * f1 + f_res, f_res) # possible periodic frequencies
f_real = []
fi_pos = np.arange(f1,frange*f1+f_res, f_res) # possible periodic frequencies
f_real=[]
for f in freq:
f_real.append(find_nearest(f, fi_pos))
freq = sp.hstack(f_real)
f_real.append(find_nearest(f,fi_pos))
freq = np.hstack(f_real)
if True:
print("freq: " + str(freq))
print("freq: " +str(freq))
if phases is None: # generate random phases
phase = sp.randn(Nf) * 2 * sp.pi # random phase
else: # use given phases
phase = phases
phase = sp.randn(Nf)*2*np.pi # random phase
else: # use given phases
phase=phases
return freq, phase, T1, fs
return freq, phase,T1,fs
###################################
# Pseudo-Random-MultiSine for "quick calibration"
def PR_MultiSine(
f1,
Nperiods,
Nsamples,
Nf=8,
fs_min=0,
fs_max=1e9,
frange=10,
log=True,
phases=None,
deg=False,
sample_inkr=1,
):
def PR_MultiSine(f1, Nperiods, Nsamples, Nf=8, fs_min=0, fs_max=1e9, frange=10, log=True, phases=None,deg=False, sample_inkr=1):
"""
Returns an additive normalized Multisine time series. \n
f1 = start frequency (may be adapted) \n
......@@ -625,47 +573,39 @@ def PR_MultiSine(
multi=array of time series values \n
"""
freq, phase, T1, fs = PR_MultiSine_adapt(
f1,
Nperiods,
Nsamples,
Nf=Nf,
fs_min=fs_min,
fs_max=fs_max,
frange=frange,
log=log,
phases=phases,
sample_inkr=sample_inkr,
)
if deg: # rad -> deg
phase = phase * sp.pi / 180.0
freq,phase,T1,fs = PR_MultiSine_adapt(f1,Nperiods,Nsamples, Nf=Nf,fs_min=fs_min,fs_max=fs_max,frange=frange,log=log, phases=phases, sample_inkr=sample_inkr)
ti = sp.arange(T1 * fs, dtype=sp.float32) / fs
if deg : # rad -> deg
phase = phase*np.pi/180.0
ti = np.arange(T1*fs,dtype=np.float32)/fs
multi = sp.zeros(len(ti), dtype=sp.float64)
for f, p in zip(freq, phase):
multi = multi + sp.sin(2 * sp.pi * f * ti + p)
multi = np.zeros(len(ti), dtype=np.float64)
for f,p in zip(freq,phase):
multi = multi + np.sin(2*np.pi*f*ti + p)
multi = multi / sp.amax(sp.absolute(multi)) # normalize
multi = multi / np.amax(np.absolute(multi)) # normalize
if False:
import matplotlib.pyplot as mp
fig = mp.figure(1)
fig.clear()
pl1 = fig.add_subplot(211)
pl2 = fig.add_subplot(212)
pl1.plot(ti, multi, "-o")
pl2.plot(sp.hstack((ti, ti + ti[-1] + ti[1])), sp.hstack((multi, multi)), "-o")
pl1.plot(ti,multi,"-o")
pl2.plot(np.hstack((ti,ti+ti[-1]+ti[1])),np.hstack((multi,multi)),"-o")
mp.show()
return freq,phase,fs,multi # frequency series, sample rate, sample timestamps, waveform
return (
freq,
phase,
fs,
multi,
) # frequency series, sample rate, sample timestamps, waveform
# PR_MultiSine(1,10,1500,5,fs_max=101,sample_inkr=7)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment