Lücken bei der Annäherung n. Bézier-Bernstein-Verfahren

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
gotridofmyphone
User
Beiträge: 33
Registriert: Mittwoch 15. März 2017, 08:54

Hallo,

mein Nicht-Echtzeit-Synthesizer (bzw. Klangdesigner) soll ermöglichen, Einschwing-, quasistationären und Ausklingphasen aller Teiltöne eines Klangs modellieren zu können, die noch mehr als die Amplitudenverhältnisse unter diesen Teiltönen ("Timbre") gemeinsam den Charakter eines Instruments bilden. Da ein natürlicher Klang viele solcher Teiltöne hat, und es eh unmöglich sein wird, ihn 100% nachzubilden, gilt es, den Spezifizierungsaufwand pro Teilton zu minimieren, um die maximale Detailtreue für den geringstmöglichen Aufwand herauszuholen. Nun soll der Benutzer nur ein zwei oder drei Koordinaten pro Teilton und Phase festlegen. Das System soll zwischen den Koordinaten eine Bezierkurve ziehen und diese "abtasten", so dass ich am Schluss alle y-Werte (Amplituden) zu allen x-Werten (0..Länge der Phase) kenne. Soweit zum Kontext meiner Frage.

Das klappt im Wesentlichen auch. Leider nicht für alle x-Werte. Manche bleiben unberücksichtigt. Kann mir jemand vielleicht bei diesem eigentlich nicht python-bezogenen, eher algorithmischen Problem helfen? Meine Vermutung ist, dass das mit der leidigen Gleitkommaberechnungsproblematik zu tun hat, aber an welcher Stelle tritt der Fehler konkret auf?

Erst mal habe ich mir beholfen, indem ich noch mal über die Liste marschiere und fehlende Werte interpoliere (Durchschnitt der beiden Nachbarwerte). Das ist natürlich nicht so schön, und sonderlich performant auch nicht.

Danke,
-- gotridofmyphone

Code: Alles auswählen

# Bezier-Basispolynome Grad 2:
# ----------------------------
# B_0,2 = (1-t)²; B_1,2 = 2*t*(1-t); B_2,2 = t²
# x(t) = x_0*B_0,2(t) + x_1*B_1,2(t) + x_2*B_2,2(t)
# 
# Bezier-Basispolynome Grad 3:
# ----------------------------
# B_0,3 = (1-t)³; B_1,3 = 3*(1-t)²*t; B_2,3 = 3*(1-t)t²; B_3,3 = t³
# x(t) = x_0*B_0,3(t) + x_1*B_1,3(t) + x_2*B_2,3(t) + x3*B_3,3(t)
# 
# B0,n wird nicht benötigt, weil der erste Koordinat immer (0,0) ist.
# 
# Beim Einschwingen:
# x_n ist die Einschwingzeit.
# y_n ist 1.
# Beim Ausklingen:
# x_n ist die Ausklingzeit.
# y_n ist 0.

import numpy as np
import matplotlib.pyplot as plt

def get_bezier_func(p1,p2,p3, tune_length):

    # scale x to tune length, y to 1
    x_max = p3[0] if p3 else p2[0]
    x_last = 0
    y_max = max(i[1] for i in (p1, p2, p3))

    for i in p1, p2, p3:
        if i is None: break

        if x_last > i[0]:
            raise Exception("Wrong order in coordinates")
        x_last = i[0]

        i[0] *= tune_length / x_max

        if not i[1] > 0:
            raise Exception("y must be greater than 0")
        i[1] /= y_max
        
    def B_12(x): return lambda t: x * 2 * t * (1-t)
    def B_22(x): return lambda t: x * t**2
    def B_13(x): return lambda t: x * 3 * (1-t)**2 * t
    def B_23(x): return lambda t: x * 3 * (1-t) * t**2
    def B_33(x): return lambda t: x * t**3
    
    if p3 is None:
       xf = [ B_12(p1[0]), B_22(p1[0]) ]
       yf = [ B_12(p1[1]), B_22(p1[1]) ]
    else:
       xf = [ B_13(p1[0]), B_23(p2[0]), B_33(p3[0]) ]
       yf = [ B_13(p1[1]), B_23(p2[1]), B_33(p3[1]) ]
    
    return lambda t: (sum( f(t) for f in xf ), sum( f(t) for f in yf ))

# Doch t wissen wir gar nicht.
# Das bedeutet, bis wir einen besseren Algorithmus finden,
# müssen t so lange abtasten, bis wir y zu allen int(x) kennen.

# Zu x=0 und x_n=x_max wissen wir y bereits vorab: Beim Einschwingen ist y=t=0 bzw. y=t=1,
# beim Abklingen verhält es sich umgekehrt, also t=0, y=1 bzw. t=1, y=0.

def get_amplitudes (approx, length):

    results = [None]*length
    results[0] = 0

    def scan_bezier(start, pos, max):
    
        x, y = approx( start + pos ); x = int(x)
        results[ x ] = y
        half_pos = pos / 2

        print "start: {} | pos: {} | max: {} => ({}, {})".format(
            start, pos, max, x, y
        )
 
        if results[ x-1 ] is None:
            scan_bezier( start, half_pos, x-1 )
        else: return
    
        if x < max:
            scan_bezier( start+pos, half_pos, max )

        return
    
    scan_bezier( 0, .5, length )

    for i, value in enumerate(results):
        if value is None:
            results[i] = (results[i-1] + results[i+1]) / 2
            print "Empty index {} interpolated to {}".format(
                i, results[i]
            )

    return results

length = 530 # 12ms bei 44,1KHz Abtastrate: typische Einschwingzeit eines Klaviertons
example = get_bezier_func([1, 2.5], [2.5, 1], [5, 3], length)

results = get_amplitudes(example, length)
plt.plot(np.arange(length), results )
plt.show()
BlackJack

@gotridofmyphone: Ohne mich jetzt tatsächlich mit dem Algorithmus beschäftigt zu haben: Die ``print``-Anweisungen bedeuten das ist für Python 2, und Du teilst da an einigen Stellen durch eine ganze Zahl: Bist Du sicher das der erste Operand entweder immer glatt teilbar oder eine Gleitkommazahl ist?

Falls nicht: Wird es besser wenn Du als ersten Import ein ``from __future__ import division`` machst?
gotridofmyphone
User
Beiträge: 33
Registriert: Mittwoch 15. März 2017, 08:54

@BlackJack: Ja, gute Idee ... hat leider keine Besserung gebracht, trotzdem danke.
gotridofmyphone
User
Beiträge: 33
Registriert: Mittwoch 15. März 2017, 08:54

Ich habe das Problem selbst gelöst bekommen, siehe Zeile 66. Außerdem wird der Binominalkoeffizient jetzt on-demand berechnet und gecacht. Für jedwede Codekritik bin ich gern zu haben, für Optimierungsvorschläge natürlich auch.

Code: Alles auswählen

    _b_cache = {}
    @staticmethod
    def _get_bezier_func(*coords):

        def b(n,k):

            v = Shape._b_cache.get( (n,k) )
            if v: return v

            if k == 0: v = 1

            elif 2*k > n:
                v = b(n,n-k)

            else:
                v = (n+1-k) / k * b(n,k-1)

            Shape._b_cache[ (n,k) ] = v
            return v

        def B(x,n,k): return lambda t: b(n,k) * x * t**k * (1-t)**(n-k)

        xf, yf, cnum = [], [], len(coords)-1
        for i, c in enumerate(coords):
            xf.append( B(c[0], cnum, i) )
            yf.append( B(c[1], cnum, i) )

        func = lambda t: (
            sum( f(t) for f in xf ),
            sum( f(t) for f in yf )
        )

        return func

    def render (self, unit_length=1, y_scale=1, adj_length=False ):
        "Get samples according to the bezier curve"

        coords = copy.deepcopy(self.coords)

        if adj_length and isinstance(adj_length, bool):
            unit_length *= y_scale
            coords = self.new_coords(y_scale=y_scale)
            length = int( unit_length * self.length + .5 )
        else:
            length = adj_length or self.length
            coords = self.new_coords(adj_length, y_scale)
            length = int( unit_length * length + .5 )

        if not length: return []

        approx = Shape._get_bezier_func(*coords)

        results = [None]*(length+1)
        results[0] = coords[0][1]
        results[-1] = coords[-1][1]

        def scan_bezier(start, pos, max):

            t = start + pos
            x, y = approx( t ); x = int(x * length)
            results[ x ] = y
            half_pos = pos / 2

            if results[ x-1 ] is None:
                scan_bezier( start, half_pos, x-1 )
            elif x == length or results[ x+1 ] is not None: # HIER war der Fehler: else reicht nicht
                return

            if x < max:
                scan_bezier( start+pos, half_pos, max )

            return

        scan_bezier( 0, .5, length )

        return results
Leider ist das Abtasten des Bezierpfades sehr langsam, in Anbetracht der Abtastrate von mehreren zehntausend Hertz pro Sekunde. So ein naturalistischer Ton eines Instruments hat zudem dutzende Teiltöne, jeder mit seiner eigenen ADSR-Charakteristik (Attack, Decay, Sustain, Release heißen die Schwingphasen komplett). Sobald es nicht mehr tolerabel für mich wird, versuche ich, das mittels Boost.Python über C++ laufen zu lassen, oder habt ihr eine bessere Idee? Und wenn auch das nichts hilft, muss ich halt doch auf lineares ADSR-Modelling umschwenken. Pro Ton ne Minute zu warten ... bei allem was recht ist. Cachen wäre noch zu überlegen, damit fertig gerenderte Töne ggf. mehrmals verwendet werden können.
Antworten