Elektroda.pl
Elektroda.pl
X

Search our partners

Find the latest content on electronic components. Datasheets.com
Elektroda.pl
Please add exception to AdBlock for elektroda.pl.
If you watch the ads, you support portal and users.

Jak obliczyć: 'b'? Chodzi o aproksymację trygonometryczną.

TvWidget 12 Aug 2021 22:13 510 21
  • #1
    TvWidget
    Level 36  
    Załóżmy, że napięcie zmienia się okresowo w czasie można go opisać funkcją U(t)=a*sin(t+b) gdzie a i b to stałe wartości. Mam zbiór próbek pomiarowych tego napięcia (są to pary liczb U i t). Jak na ich podstawie obliczyć b ?
  • #2
    khoam
    Level 41  
    Skoro b to stała wartość, to po co masz obliczać b?
  • #3
    krzysiek_krm
    Level 40  
    TvWidget wrote:
    Załóżmy, że napięcie zmienia się okresowo w czasie można go opisać funkcją U(t)=a*sin(t+b) gdzie a i b to stałe wartości. Mam zbiór próbek pomiarowych tego napięcia (są to pary liczb U i t). Jak na ich podstawie obliczyć b ?

    Ja bym kombinował z analizą miejsc zerowych.
  • #4
    khoam
    Level 41  
    b = arcsin(U/a) - t
    :)
  • #5
    piotrek0207
    Level 20  
    Jeszcze trzeba pamiętać, że argumentem sinusa nie jest tylko t, a omega*t - pulsacja*t. Musimy znać częstotliwość przebiegu. B to będzie jakieś przesunięcie fazowe.
  • #6
    Maciej Gonet
    VBA, Excel specialist
    TvWidget, dlaczego nie dałeś przykładowych danych?
    Czy jesteś pewien, że w Twoim wzorze nie brakuje częstości ω?
    Generalnie jest to problem regresji nieliniowej i np. w Excelu można go rozwiązać za pomocą Solvera.
    Przykład w załączniku. Użyłem danych symulowanych, do których wprowadziłem losowe zaburzenie wartości U. Formuła z komórki B2 została użyta w kolumnie B, a następnie zamieniona na wartości.
    Do tego dopasowałem gładką krzywą w kol. C. Sumę kwadratów odchyleń zminimalizowałem Solverem, zmieniając wartości w komórkach I1:I2. Gdyby jednak była potrzebna też ω, to trzeba dodać trzeci parametr w analogiczny sposób.
    Na wykresie widać punkty użyte do symulacji (niebieskie) i dopasowaną linię (pomarańczową).
  • #7
    TvWidget
    Level 36  
    Nie potrzebuję policzyć tego w Excelu tylko zaimplementować algorytm w uP.
    Częstotliwość ω to jedynie parametr. To czy we wzorze jest ωt czy samo t nie ma znaczenia.
    Zdaje się, że kiedyś rozwiązanie identycznego problemu znalazłem w książce "Metody numeryczne". Algorytm był dość prosty. Wartości próbek były jakoś mnożone przez cos(t) i sin(t). Coś było sumowane i uzyskiwało się wartości a i b.
    W transformacie Fouriera sygnał można przybliżyć sumą funkcji okresowych. W tym wypadku sprawa jest prostsza. Chodzi o przybliżenie tylko jedną funkcją okresową.
  • #8
    _jta_
    Electronics specialist
    TvWidget wrote:
    Wartości próbek były jakoś mnożone przez cos(t) i sin(t).

    Dokładnie. Potem to się sumuje (oddzielnie z mnożenia przez cos(t) i sin(t)), dostajesz dwie liczby - nazwijmy je Uc i Us. Dalsze dwie wyliczasz sumując kwadraty cos(t) i sin(t) dla tych samych 't' - one będą potrzebne do znormalizowania pierwszych dwu, nazwijmy je Nc i Ns. Na koniec liczysz a=sqrt((Uc/Nc)^2+(Us/Ns)^2), b=atan2(Us/Ns,Uc/Nc). Myślę, że to zadziała.

    Okazuje się, że to jest bardziej skomplikowane, i trzeba wyliczyć jeszcze jedną sumę: cos(t)*sin(t). Teraz sumy bez U nazwałem 'cc', 'cs', 'ss', i zamiast 'U' mam 'v'.
    Code: tcl
    Log in, to see the code

    Na końcu programu dla sprawdzenia policzyłem iloczyn macierzy - to wymaga odwracania macierzy:
    Code: less
    Log in, to see the code
    (jest symetryczna 2x2, więc łatwo - poprzestawiać, zmienić znak 'cs' i podzielić przez wyznacznik)
  • #9
    khoam
    Level 41  
    @jta Ten soft to ma pracować na µP, a nie wielordzeniowym Intel czy AMD :)
  • #10
    Maciej Gonet
    VBA, Excel specialist
    Tej metody, którą podał _jta_, nie znałem. Daje ona dobre wyniki aproksymacji, choć nieco inne niż przy metodzie najmniejszych kwadratów. Widocznie tam zastosowano inne kryterium.
    Metodę najmniejszych kwadratów można też zastosować bezpośrednio, bez pomocy Solvera. Opiera się to na wykorzystaniu wzoru na sinus sumy:
    Code: text
    Log in, to see the code

    W Excelu trzeba utworzyć kolumny sin(t) i cos(t) i potraktować je jak zmienne niezależne w analizie regresji liniowej. Otrzymuje się wyniki takie same jak przy użyciu Solvera.
    Jest też jeszcze metoda przybliżona, w której amplitudę a szacuje się na podstawie różnicy między wartością maksymalną a minimalną, a wartość b na podstawie punktu przejścia funkcji przez 0. Ta metoda wymaga, żeby przebieg obejmował maksimum i minimum sinusoidy oraz żeby nie był mocno zakłócony.
    Trzeba też pamiętać, że b jest określone z dokładnością do stałej, bo funkcja sinus jest okresowa. Z reguły wybiera się przedział [-pi, +pi] albo [0, 2*pi]. Nie podałeś, który to ma być wariant.
    W pliku przykłady w Excelu. To musisz sobie przetłumaczyć na swój język.
  • #11
    StaryVirus_e_Wiarus
    Level 20  
    Cześć
    Kompilator powinien mieć możliwość i udostępniać bibliotekę <math.h>. W niej jest zapewne udostępniona funkcja "arcsin(X)". Wzór z postu #4 miałby tu zastosowanie.
  • #13
    khoam
    Level 41  
    @staryvirusewiarus Nawet nie trzeba używać funkcji bilbiotecznej arcsin(). Wystarczy aproksymacja:

    Jak obliczyć: 'b'? Chodzi o aproksymację trygonometryczną.
  • #16
    _jta_
    Electronics specialist
    W tym artykule pomiary napięcia wykonuje się w równych odstępach czasu ∆T - czy takie ma być założenie? Mój sposób działa bez takiego wymogu, odstępy czasu mogą być prawie dowolne (byle nie wielokrotności połowy okresu, bo wtedy wyznacznik macierzy wyjdzie 0, a trzeba wykonać dzielenie przez ten wyznacznik; a najlepiej, jakby były po ćwierć okresu, wtedy można prościej). Ale za to w tym artykule wyznacza się fazę i amplitudę zakłóceń, aby odjąć zakłócenia od sygnału użytecznego - tego z kolei nie daje mój sposób.
  • #17
    rb401
    Level 38  
    TvWidget wrote:
    W transformacie Fouriera sygnał można przybliżyć sumą funkcji okresowych. W tym wypadku sprawa jest prostsza. Chodzi o przybliżenie tylko jedną funkcją okresową.


    No to w sam raz do tego celu pasuje Ci algorytm Goertzela, taki właśnie "jednoprążkowy FFT". Tym bardziej że jest bardzo mało obciążający obliczeniowo, nie wymaga konieczności zbierania próbek w tablicy, ani innych operacji na tablicach, w sam raz nawet na słabe uC.

    Tu ładnie opisany algorytm i praktyczne wskazówki o dobraniu wielkości okna:

    https://www.embedded.com/the-goertzel-algorithm/

    a tu źródła do przykładu z tego artykułu:

    https://github.com/pramasoul/jrand/blob/master/goertzel.c

    Z ciekawości dodałem tam tylko instrukcję "angle = atan2(imag, real);" i miałem policzoną fazę.
    Z testów mi wyszło że regulując fazę w funkcji generacji próbek, z algorytmu dostaję w wyniku precyzyjnie (do trzeciego miejsca po przecinku dla kąta w radianach) fazę która zadaję.
    Jedyny problem jaki widzę to konieczność "skalibrowania" zera fazy względem początku bloku próbek (lub innego arbitralnie wybranego punktu bloku próbek).
    W tym przykładzie jeśli blok próbek zaczynał się od próbki sin(0), jak oryginalnie to robi funkcja Generate(), w wyniku pozyskanym z algorytmu pomierzona faza nie wychodziła zerowa. Akurat dla danych w tym przykładzie gdy dodałem do generacji sinusa stały offset równy 1.9544rad (nie wnikam dlaczego akurat tyle) to już wartości fazy zadane w generacji próbek i pomierzone przez algorytm były liczbowo identyczne.
  • #18
    StaryVirus_e_Wiarus
    Level 20  
    Można też wartości kątów i sinusa (arcsin) ująć w tabelce dwuwymiarowej i prostą instrukcją przeszukiwać taką tabelę. A najbliższe wartości aproksymować lub uśrednić do najbliższej właściwego wartości. Trzeba to wykonać z założonym błędem obliczeń.
  • #19
    _jta_
    Electronics specialist
    Wciąż pozostaje otwarta sprawa: jakie są odstępy czasów pomiędzy próbkami? Jeśli można mieć stałe i dostosować je do potrzeby wyliczania, to jest prościej.

    Przykładowe wyniki dla mojego programu (częstotliwość jest 1, częstość 2pi): zadane a=12 b=0.123 t1=0.1 t2=0.4; wyliczone a=12.0 b=0.12300000000000001.

    Błąd na poziomie 10^(-17) - ale tak wychodzi, jak do wyliczania podaje się dokładne wygenerowane wartości, a różnica czasów jest bliska T/4.
    Przykładowy wynik dla "gorszych" parametrów: a=12 b=0.123 t1=0.06 t2=0.10 t3=0.55 => a=12.000000000000027 b=0.12300000000000019.

    Jest jeszcze kwestia, czy stablicowanie sin() i cos(), a potem interpolacja (albo użycie wzorów na sin() i cos() sumy kątów, dla małych kątów można przyjąć sin(x)=x, cos(x)=1-x^2/2) daje oszczędność czasu procesora przy wyliczaniu, czy nie. Dla kąta 30° znamy ze szkoły sin() i cos(); dla kąta 18° można je wyliczyć z konstrukcji pięciokąta foremnego; wykorzystując wzory na sin() i cos() sumy i różnicy kątów, oraz połowy kąta można je stablicować np. co 3°.
  • #20
    TvWidget
    Level 36  
    Chyba udało mi się znaleźć algorytm. Wydaje się jednak dziwnie prosty.
    Załóżmy, że mamy 20 próbek danych. W poniższym kodzie generuje te 20 próbek testowo zmieniając fazę o 0.1.
    Różnica pomiędzy przyjętą wartością fazy a obliczoną nie przekracza 0.01.
    for (fi=0;fi<3;fi=fi+0.1){
    a=0;
    b=0;
    for (i=1;i<20;i++){
    y=Math.sin(i+fi); //testowe próbki

    a=a+y*Math.cos(i)
    b=b+y*Math.sin(i)
    }
    console.log(fi-Math.atan2(a,b));
    }
    }
  • #21
    _jta_
    Electronics specialist
    To jest rzeczywiście za proste. Działa w miarę dobrze, jeśli sumy kwadratów sin(i) i cos(i) są zbliżone, a suma ich iloczynów jest bliska 0.

    Jak zrobisz próbkowanie sygnału co π/4, i weźmiesz parzystą liczbę próbek, to wynik będzie dokładny (na ile dokładne będą obliczenia).
  • #22
    rb401
    Level 38  
    TvWidget wrote:
    Chyba udało mi się znaleźć algorytm. Wydaje się jednak dziwnie prosty.


    Ok, działa, bo w zasadzie musi. Liczysz po prostu korelacje próbek z zadanym sinusem i kosinusem, tak z definicji.
    Ale jednak masz jeszcze dość ciężki obliczeniowo algorytm. Na każdą próbkę liczysz sinus i kosinus i masz dwa mnożenia.
    Co prawda sinus można wstępnie stablicować.
    Ale jeśli to ma chodzić na mniej wydajnym sprzęcie to jednak Twój algorytm mocno przegrywa z Goertzelem, który to samo wylicza ale na każdą próbkę jest potrzebne tylko jedno mnożenie a żadnej tablicy nie potrzeba. Wystarczą tylko trzy stałe jako kompletna nastawa pomiaru. No i też Goertzel jest opisany w pełni, sprawdzony przez lata i mnóstwo o nim informacji, i przykładów implementacji.