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

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 ?

TvWidget napisał:

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.

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.

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ą).

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ą.

TvWidget napisał:

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'.
Na końcu programu dla sprawdzenia policzyłem iloczyn macierzy - to wymaga odwracania macierzy:(jest symetryczna 2x2, więc łatwo - poprzestawiać, zmienić znak 'cs' i podzielić przez wyznacznik)

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:

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.

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.

khoam napisał:

Ten soft to ma pracować na µP, a nie wielordzeniowym Intel czy AMD

Zawsze można przetłumaczyć na C, np. set s [expr { sin($tpi*$t) }]; => s=sin(tpi*t);. Ale trzeba mieć funkcje trygonometryczne i pierwiastek kwadratowy.

Należy używać atan2(), a nie arcsin() - bo w okolicach π/2 arcsin() da duży błąd. Do tego, może dać zły znak 'b'.

Znalazłem w sieci artykuł na podobny temat:
http://pe.org.pl/articles/2016/5/26.pdf
Może będzie przydatny lub nawiążesz kontakt z autorami.

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.

TvWidget napisał:

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.

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ń.

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°.

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));
}
}

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).

TvWidget napisał:

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.