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

added 'periods' parameter to sequential fits, added XV_MultiSine()

parent a0e40b1f
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 358 additions and 227 deletions
SineTools.py
+ 358
− 227
View file @ 037ca90b
......@@ -4,15 +4,16 @@ Created on Fri Jun 14 20:36:01 2013
SineTools.py
auxiliary functions related to sine-approximation
@author: bruns01
@author: Th. Bruns, B.Seeger
"""
import scipy as sp
import numpy as np
from scipy import linalg as la
import matplotlib.pyplot as mp
import warnings
def sampletimes(Fs, T): #
def sampletimes(Fs, T): #
"""
generate a t_i vector with \n
sample rate Fs \n
......@@ -23,7 +24,7 @@ def sampletimes(Fs, T): #
# 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,18 +39,25 @@ 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 * np.random.randn(len(ti))
if absnoise != 0:
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
n = n + absnoise * np.random.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
......@@ -65,20 +73,30 @@ def fm_counter_sine(fm, f, x, phi, ti, offset=0, noise=0, absnoise=0, drift=0, a
"""
Tau = ti[-1] - ti[0]
n = 0
if noise != 0:
n = x*noise*sp.randn(len(ti))
if noise != 0:
n = x * noise * np.random.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)*np.sin(2*np.pi*f * ti - phi) + n +d*ti + offset)
s = np.floor(s+fm*ti)
n = n + absnoise * np.random.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
......@@ -88,15 +106,16 @@ 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*np.pi*f0
a = np.array([np.cos(w0*t),np.sin(w0*t),np.ones(t.size)])
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
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
......@@ -106,20 +125,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*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
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*np.pi*f0
return abc[0]*np.cos(w0*t)+abc[1]*np.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):
......@@ -127,7 +146,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 np.absolute(abc[1]+1j*abc[0])
return np.absolute(abc[1] + 1j * abc[0])
def phase(abc, deg=False):
......@@ -136,9 +155,10 @@ 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 np.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
......@@ -147,7 +167,8 @@ 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
......@@ -157,32 +178,42 @@ 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, periods=1):
"""
period-wise sine-fit at a known frequency\n
y vector of sample values \n
t vector of sample times\n
f0 known frequency\n
periods number of full periods for every fit (default=1)\n
\n
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(np.floor(t.size/N)) ## number of sections or periods
if y.size < t.size:
raise ValueError("Dimension mismatch in input data y<t")
if t.size < y.size:
warnings.warn(
"Dimension mismatch in input data y>t. fiting only for t.size values",
RuntimeWarning,
)
Tau = 1.0 / f0 *periods
dt = t[1] - t[0]
N = int(Tau // dt) ## samples per section
M = int(t.size // N) ## number of sections or periods
abc = np.zeros((M, 3))
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
......@@ -192,34 +223,39 @@ 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*np.pi*f0
abcd = threeparsinefit(y, t, f0)
w = 2 * np.pi * f0
err = 1
i = 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)
])
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]
abcd = (la.lstsq(D.transpose(), y))[0]
dw = abcd[3]
w = w+0.9*dw
i+=1
err = np.absolute(dw/w)
w = w + 0.9 * dw
i += 1
err = np.absolute(dw / w)
assert i < nmax, "iteration error"
return np.hstack((abcd[0:3],w/(2*np.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*np.pi*abcf[3]
return abcf[0]*np.cos(w0*t)+abcf[1]*np.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)
......@@ -244,52 +280,62 @@ 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, periods=1, debug_plot=False):
"""
period-wise sine-fit at a known frequency\n
sliced or period-wise sine-fit at a known frequency\n
y vector of sample values \n
t vector of sample times\n
f0 estimate of excitation frequency\n
periods integer number of periods used in each slice for fitting
nmax maximum of iteration to improve f0 \n
debug_plot Flag for plotting the sequential fit for dubugging \n
\n
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(np.floor(Tau/dt)) ## samples per section
M = int(np.floor(t.size/N)) ## number of sections or periods
if y.size < t.size:
raise ValueError("Dimension mismatch in input data y<t")
if t.size < y.size:
warnings.warn(
"Dimension mismatch in input data y>t. fiting only for t.size values",
RuntimeWarning,
)
Tau = 1.0 / f0 *periods
dt = t[1] - t[0]
N = int(Tau // dt) ## samples per section
M = int(t.size // N) ## number of sections or periods
abcd = np.zeros((M, 4))
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): # f0 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
......@@ -300,36 +346,36 @@ def multi_threeparsinefit(y,t,f0): # f0 vector of frequencies
for y = sum_i (a_i*sin(2*pi*f0_i*t) + b_i*cos(2*pi*f0_i*t) + c_i
"""
w0 = 2 * np.pi * f0
D = np.ones((len(t),1)) # for the bias
D = np.ones((len(t), 1)) # for the bias
# set up design matrix
for w in w0[::-1]:
D = np.hstack(( np.cos(w*t)[:,np.newaxis],
np.sin(w*t)[:,np.newaxis],
D))
D = np.hstack((np.cos(w * t)[:, np.newaxis], np.sin(w * t)[:, np.newaxis], D))
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[:-1][0::2] + 1j*abc[:-1][1::2] # make complex without Bias (last value)
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][0::2] + 1j*abc[:-1][1::2] # make complex without Bias (last value)
return np.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
......@@ -339,12 +385,13 @@ 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*np.cos(2*np.pi*fi*t)+b*np.sin(2*np.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):
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
......@@ -357,19 +404,20 @@ def multi_waveform_mp(f,m,p,t, bias=0, deg=True):
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
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
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
......@@ -381,21 +429,22 @@ 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(np.floor(Tau/dt)) ## samples per section
M = int(np.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 = np.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):
"""
......@@ -403,18 +452,21 @@ 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 :
if diff:
d = np.diff(y)
d = d - np.mean(d)
y = np.hstack((0,np.cumsum(d)))
y = np.hstack((0, np.cumsum(d)))
else:
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
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=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):
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
......@@ -426,41 +478,55 @@ 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*np.pi*f
om2 = om**2
tau = ti[-1]-ti[0]
disp = np.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*np.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 = np.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)*np.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*np.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*np.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*np.cos(om*ti+phi0))/tau
i=i+1
return disp,acc
acc = (
acc
+ (2 * ampdrift * s0 * om * i * h * np.cos(om * ti + phi0)) / tau
)
i = i + 1
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 by the algorithm) \n
......@@ -481,76 +547,99 @@ def PR_MultiSine_adapt(f1, Nperiods, Nsamples, Nf=8, fs_min=0,fs_max=1e9, frange
T1 = (adapted) duration
fs=sample rate \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 = 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
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 = 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)
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
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 = np.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))
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*np.pi # random phase
else: # use given phases
phase=phases
phase = np.random.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
......@@ -573,39 +662,81 @@ def PR_MultiSine(f1, Nperiods, Nsamples, Nf=8, fs_min=0, fs_max=1e9, frange=10,
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)
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*np.pi/180.0
ti = np.arange(T1*fs,dtype=np.float32)/fs
if deg: # rad -> deg
phase = phase * np.pi / 180.0
ti = np.arange(T1 * fs, dtype=np.float32) / fs
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)
for f, p in zip(freq, phase):
multi = multi + np.sin(2 * np.pi * f * ti + p)
multi = multi / np.amax(np.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(np.hstack((ti,ti+ti[-1]+ti[1])),np.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
###################################
# MultiSine for "quick calibration" based on displacement and velocity without
# explicit time-information (used for Linearmotor-control)
def XV_MultiSine(freq, phase, x0, tau, N=1000, deg=True):
"""
returns two arrays of positions and velocities of length N which describe
a multi-sine motion with the frequencies given in f and displacement amplitudes
given in x0.
----
freq numpy array of frequencies\n
phase numpy array of initial phases\n
x0 numpy array of displacement amplitudes\n
tau double giving the duration of the motion
N integer length of result array, len(x) and len(v)\n
deg boolean True:phase values in degree, False:phase values in rad
"""
assert (
len(freq) == len(x0) == len(phase)
), "XV_MultiSine: Unequal length of frequency and amplitude arrays!"
if deg: # rad -> deg
phase = phase * np.pi / 180.0
ti = np.arange(0, tau, tau / N)
multi_x = np.zeros(len(ti), dtype=np.float64)
multi_v = np.zeros(len(ti), dtype=np.float64)
for f, x, p in zip(freq, x0, phase):
om = 2 * np.pi * f
multi_x = multi_x + x * np.sin(om * ti + p)
multi_v = multi_v + x / om * np.cos(om * ti + p)
return multi_x, multi_v
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment