Как определить параметры синуса?

[shadeofgray]

Цитата ANDLL @ 12.10.09, 10:29

А то результат все равно довольно неправильный(

Похоже, проблема в свойствах функции - много локальных минимумов, метод попадает в них. Я так понимаю, ему очень сложно менять частоту синуса, т.к. при этом меняется знак функции и проходится большое количество минимумов, в которых мы и застреваем.

Есть например такой вариант - использовать спектральный метод, предложенный Prince, для определения начального приближения. Затем, если качество не устроит - Левенберг-Марквардт для окончательной подгонки.

[Prince]

Цитата

К слову сказать, точки и правда лежат на равномерной сетке

Тогда,

Цитата

Во-вторых, спектральные методы позволяют выбирать множитель "B" только из конечного диапазона частот, а Левенберг-Марквардт позволяет непрерывно выбирать оптимальный параметр. В-третьих, спектральный анализ совершенно не годится для ситуации, когда нельзя однозначно определить значение функции в точке (есть несколько точек с совпадающими абсциссами, но разными значениями функции) - ситуация, часто встречающаяся в практике.

Второй пункт может оказаться несущественным, ведь неизвестна пока необходимая точность, кроме того, погрешность измерений и искажения самой синусоиды позволят ли однозначно выбрать "оптимальный" параметр? Я не знаком абсолютно с рассматриваемым методом, так что вполне серьезно интересуюсь. Если, скажем, в выборке кроме основного "тона" присуствуют гармоники? Можно ли подобрать оптимальные параметры?
Последнее замечание я не понимаю. Как может быть несколько значений y для одной абсциссы? Каждой точке х соответсвует одна и только одна точка y.

Сообщение отредактировано: Prince - 12.10.09, 10:47

[ANDLL]

А ссылку можно на спектральный метод?
http://www.google.ru/search?q=%D1%81%D0%BF%D0%B5%D0%BA%D1%82%D1%80%D0%B0%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9+%D0%BC%D0%B5%D1%82%D0%BE%D0%B4&ie=utf-8&oe=utf-8&aq=t&rls=org.mozilla:ru:official&client=firefox показывает мало интересного

Добавлено 12.10.09, 10:46
shadeofgray
А есть алгоритмы, которые решают задачу не поиска минимума функции, а поиска нуля(или близкого к нему) функций такого вида?

[Pavia]

А я не согласен преобразование Фурье спокойно работает с неравномерной сеткой. Только быстрое нельзя организовать.
По поводу разных значений y при одинаковом x. Оно усреднится.

[shadeofgray]

Цитата Prince @ 12.10.09, 10:43

Последнее замечание я не понимаю. Как может быть несколько значений y для одной абсциссы? Каждой точке х соответсвует одна и только одна точка y.

Например, измерения в одной точке проводились несколько раз, и каждый раз из-за шума получались разные значения.

Цитата ANDLL @ 12.10.09, 10:44

shadeofgray
А есть алгоритмы, которые решают задачу не поиска минимума функции, а поиска нуля(или близкого к нему) функций такого вида?

Есть алгоритмы решения систем нелинейных уравнений, но я не уверен, что их можно применять в данном случае. И что имеет смысл применять. Существование минимума гарантировано, существование ноля - нет.

Мне кажется, лучше всего вариант с двухступенчатым алгоритмом - сначала взять БПФ, и выбрать гармонику с максимальной амплитудой (и ещё константу). Затем довести до оптимума при помощи ЛМ.

[ANDLL]

Так, а для тех, кто в танке можно?

[Prince]

Цитата

Например, измерения в одной точке проводились несколько раз, и каждый раз из-за шума получались разные значения.

А-а, ещё неизвестно, каким образом получаются пары значений икс игрек? То есть могут быть единичные измерения, а может многократные, с последующим "усреднением".
В любом случае, имхо, как раз спектральный анализ позволит исключить хотя бы часть шума, вне гармоник Фурье, на которые приходится полезный сигнал, даже при "небольшом" соотношении сигнал/шум, в то время как вот тот другой метод, который вы обсуждали, скорее всего, загнется.

Добавлено 12.10.09, 11:45
ANDLL, литературу можно брать отсюда:
http://dsp-book.narod.ru/books.html

[Pavia]

ANDLL
Так как сетка по x равномерная то можно использовать БПФ (FFT). То берешь свой массив Y. Заносишь в реальную часть. Делаешь преобразование.
http://alglib.sources.ru/fasttransforms/fft.php

После преобразование будет массив комплексный чисел. Амплитуда вычисляется так
A[i]=sqrt(Re[i]*Re[i]+Im[i]*Im[i])
A[0] даст смещение y0
Ищешь максимум. Начиная со 2 или 3 до n/2
Где максим A будет соответствовать амплитуде твоей синусоиды.

В этой точке ArcTan(Re/Im) даст начальную фазу x0 только надо привести в соответствии с x.
i - номер максимального отсчета даст частоту B только надо привести в соответствии с x.

Сообщение отредактировано: Pavia - 12.10.09, 12:07

[ANDLL]

Положим у меня такой массив точек:

Цитата

142.4831; 142.8679; 143.2807; 143.7201; 144.1843; 144.6715; 145.1797; 145.7068; 146.251; 146.81; 147.3815; 147.9633; 148.5532; 149.1489; 149.7478; 150.3478; 150.9464; 151.5413; 152.13; 152.7104; 153.2802; 153.8367; 154.3784; 154.9027; 155.4076; 155.8913; 156.3515; 156.787; 157.1954; 157.5756; 157.9259; 158.2448; 158.5312; 158.7838; 159.002; 159.1845; 159.3308; 159.4403; 159.5125; 159.5471; 159.5442; 159.5035; 159.4252; 159.3099; 159.1578; 158.9695; 158.7459; 158.4876; 158.196; 157.872; 157.5169; 157.132; 156.7193; 156.2798; 155.8157; 155.3285; 154.8203; 154.2932; 153.7491; 153.1901; 152.6185; 152.0367; 151.4468; 150.8511; 150.2522; 149.6521; 149.0536; 148.4587; 147.87; 147.2895; 146.7199; 146.1633; 145.6217; 145.0973; 144.5924; 144.1088; 143.6485; 143.213; 142.8046; 142.4244; 142.0742; 141.7552; 141.4688; 141.2162; 140.9981; 140.8155; 140.6693; 140.5597; 140.4876; 140.4528; 140.4559; 140.4966; 140.5748; 140.6901; 140.8422; 141.0305; 141.2541; 141.5123; 141.804; 142.128; 142.4831; 142.8679; 143.2807; 143.7202; 144.1843; 144.6715; 145.1797; 145.7068; 146.251; 146.81; 147.3815; 147.9633; 148.5532; 149.1488; 149.7479; 150.3479; 150.9465; 151.5413; 152.13; 152.7104; 153.2802; 153.8367; 154.3783; 154.9027; 155.4076; 155.8913; 156.3516; 156.787; 157.1954; 157.5755; 157.9259; 158.2448; 158.5312; 158.7838; 159.002; 159.1845; 159.3308; 159.4403; 159.5125; 159.5472; 159.5442; 159.5035; 159.4252; 159.3099; 159.1578; 158.9695; 158.7459; 158.4877; 158.196; 157.872; 157.5169; 157.1321; 156.7193; 156.2797; 155.8157; 155.3285; 154.8203; 154.2932; 153.7491; 153.1901; 152.6185; 152.0367; 151.4468; 150.8511; 150.2522; 149.6521; 149.0535; 148.4587; 147.87; 147.2895; 146.7199; 146.1633; 145.6216; 145.0973; 144.5924; 144.1088; 143.6485; 143.213; 142.8046; 142.4244; 142.0741; 141.7552; 141.4688; 141.2162; 140.9981; 140.8155; 140.6693; 140.5597; 140.4876; 140.4529; 140.4559; 140.4966; 140.5747; 140.6901; 140.8422; 141.0305; 141.2541; 141.5123; 141.804; 142.128; 142.4831; 142.8679; 143.2807; 143.7202; 144.1843; 144.6715; 145.1797; 145.7068; 146.2509; 146.8099; 147.3815; 147.9633; 148.5532; 149.1488; 149.7479; 150.3478; 150.9465; 151.5413; 152.13; 152.7105; 153.2802; 153.8368; 154.3784; 154.9027; 155.4076; 155.8913; 156.3515; 156.787; 157.1954; 157.5756; 157.9259; 158.2448; 158.5312; 158.7838; 159.002; 159.1845; 159.3308; 159.4403; 159.5125; 159.5472; 159.5442; 159.5035; 159.4252; 159.3099; 159.1578; 158.9694; 158.7459; 158.4876; 158.196; 157.872; 157.5169; 157.1321; 156.7193; 156.2797; 155.8157; 155.3285; 154.8203; 154.2932; 153.749; 153.1901; 152.6185; 152.0367; 151.4468; 150.8511; 150.2522; 149.6521; 149.0536; 148.4586; 147.87; 147.2896; 146.7199; 146.1633; 145.6216; 145.0973; 144.5924; 144.1087; 143.6484; 143.213; 142.8046; 142.4244; 142.0741; 141.7552; 141.4688; 141.2161; 140.9981; 140.8155; 140.6693; 140.5597; 140.4876; 140.4528; 140.4559; 140.4966; 140.5747; 140.6901; 140.8422; 141.0305; 141.2541; 141.5123; 141.804; 142.128; 142.4831; 142.8679; 143.2807; 143.7202; 144.1843; 144.6715; 145.1797; 145.7068; 146.251; 146.81; 147.3815; 147.9633; 148.5532; 149.1489; 149.7479; 150.3479; 150.9465; 151.5413; 152.13; 152.7104; 153.2802; 153.8367; 154.3784; 154.9027; 155.4076; 155.8913; 156.3515; 156.787; 157.1954; 157.5756; 157.9259; 158.2448; 158.5312; 158.7838; 159.0019; 159.1845; 159.3307; 159.4403; 159.5125; 159.5471; 159.5442; 159.5035; 159.4252; 159.3099; 159.1578; 158.9695; 158.7459; 158.4877; 158.196; 157.872; 157.5169; 157.1321; 156.7193; 156.2798; 155.8157; 155.3285; 154.8203; 154.2932; 153.7491; 153.1901; 152.6185; 152.0367; 151.4468; 150.8511; 150.2522; 149.6521; 149.0535; 148.4586; 147.87; 147.2895; 146.7198; 146.1632; 145.6216; 145.0973; 144.5924; 144.1088; 143.6484; 143.213; 142.8046; 142.4244; 142.0741; 141.7553; 141.4688; 141.2162; 140.9981; 140.8155; 140.6693; 140.5597; 140.4876; 140.4528; 140.4559; 140.4966; 140.5747; 140.69; 140.8422; 141.0305; 141.2541; 141.5123; 141.804; 142.128; 142.4831; 142.8679; 143.2808; 143.7202; 144.1843; 144.6715; 145.1797; 145.7068; 146.251; 146.8099; 147.3815; 147.9633; 148.5532; 149.1489; 149.7478; 150.3478; 150.9464; 151.5413; 152.13; 152.7104; 153.2801; 153.8367; 154.3784; 154.9027; 155.4076; 155.8913; 156.3516; 156.787; 157.1954; 157.5756; 157.9259; 158.2448; 158.5312; 158.7838; 159.002; 159.1845; 159.3308; 159.4403; 159.5125; 159.5472; 159.5441; 159.5035; 159.4252; 159.3099; 159.1578; 158.9695; 158.7459; 158.4877; 158.1961; 157.872; 157.5169; 157.132; 156.7193; 156.2797; 155.8157; 155.3285; 154.8203; 154.2933; 153.7491; 153.1902; 152.6185; 152.0367; 151.4468; 150.8511; 150.2522; 149.6521; 149.0535; 148.4587; 147.87; 147.2896; 146.7199; 146.1632; 145.6216; 145.0973; 144.5924; 144.1088; 143.6485; 143.213; 142.8046; 142.4245; 142.0741; 141.7553; 141.4688; 141.2162; 140.9981; 140.8155;

Где x(0) берется при y = 50, x(1) при y = 100, x(2) при y = 150, ну вы поняли
Я делаю над ним fftc1d(из алголибовской реализации)
Потом определяю что максимум a(i) достигается при i = 5
В этой точке a(i) = 2217.1317596732765
Что даже по данным видно - совсем не амплитуда.
Почему?
Код показывать?

[Pavia]

Код можно и показать. Лучше проектом.

[ANDLL]

Вот
Запускать проект ConsoleApplication1. Там нет одного проекта, но и без него будет запускаться

Добавлено 12.10.09, 16:09
Вот просто код:

    Const SOURCEFILE As String = "......."
 
    Sub Main()
        Dim rz() As String = SOURCEFILE.Split(CChar(";"))
        Dim p(rz.Length) As AP.Complex
        For i As Integer = 0 To rz.Length - 1
            p(i) = Double.Parse(rz(i))
        Next
 
        alglib.fft.fftc1d(p, rz.Length)
 
        Dim a(rz.Length) As Double
        For i As Integer = 0 To rz.Length - 1
            a(i) = Math.Sqrt(p(i).x * p(i).x + p(i).y * p(i).y)
        Next
 
        Dim maxI As Integer = 3
        For i As Integer = 3 To (rz.Length - 1) \ 2
            If a(i) > a(maxI) Then
                maxI = i
            End If
        Next
 
    End Sub

Прикреплённый файлWindowsApplication1.zip (63.71 Кбайт, скачиваний: 111)

[Pavia]

shadeofgray
К тебе притензия.

Надо делить на N при прямом преобразовании, а не при обратном. Чтобы размерности совпадали.

ANDLL
Предлагаю просто поделить на N как Re так Im части p.

Тогда A[0] равно y0=150 что правильно.
A[5] является максимум тоже верно так как у нас 5 полных синусоид.
A[5]=4.56 так как у нас есть и отрицательная частота то в два раза меньше.
А еще я так понимаю ошибка порядка 1/5=0.2
B:=2*pi/n*MaxI / 50; 50 - твой шаг
И фазу сам определишь через arctan2.

Сообщение отредактировано: Pavia - 12.10.09, 17:20

[ANDLL]

Цитата Pavia @ 12.10.09, 17:12

B:=2*pi/n*MaxI / 50; 50 - твой шаг

Pavia
Только B = 1590, а по формуле B:=1 / (pi/n*MaxI / 50) получаются, в зависимости от выборки числа в разбросе от 1100 до 1900.

[Pavia]

1/(pi/486*4/50) =1933,7325
1/(pi/486*5*/50) =1546,986
1/(pi/486*6/50) =1289,1550
1/(pi/486*7/50)=1104,990033
Точность растет чем больше частота.

Так что надо уточнять. Либо используем предложенный метод многомерной оптимизации выше. Лучше его в качестве начальных значений бери что получили после спектрального анализа.


Либо использовать прямое ПФ.

/была ошибка по правил/

function FT(w:Real; N:Integer; a:PComplex):Complex;Overload;
var i:Integer;
begin
Result.Re:=0;
Result.Im:=0;
for i:=0 to N-1 do
 begin
 Result.Re:=Result.Re+A[i].Re*Cos(-2*Pi*i*w/N)-A[i].Im*Sin(-2*Pi*i*w/N);
 Result.Im:=Result.Im+A[i].Re*Sin(-2*Pi*i*w/N)+A[i].Im*Cos(-2*Pi*i*w/N);
 end;
Result.Re:=Result.Re/N;
Result.Im:=Result.Im/N;
end;

Подставляешь w чем точнее тем больше будет результат FT. И дальше метод градиентного спуска или любой др. Правда дольше работать будет.

Сообщение отредактировано: Pavia - 12.10.09, 20:15

[ANDLL]

Цитата Pavia @ 12.10.09, 18:51

Точность растет чем больше частота.

Я и говорю, получается что ничего подобного. При увеличении числа точек точность вовсе не растет, а плавает