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

Есть результаты физического эксперимента, которые с некоторой ошибкой могут быть представлены в виде функции
y = A*sin(B*x + x0) + y0
Всего пар x;y около 1000, причем в эксперимент попадает несколько(8-10) периодов функции.
Иными словами точки насажены довольно плотно.

Вопрос - каким образом можно определить A, B, x0, y0 так, что бы точность определения параметров росла с ростом числа известных точек, ну и желательно что бы минимизировалась ошибка(но уже необязательно)?

Это называется "нелинейный МНК". Используй Левенберга-Марквардта для минимизации суммарной ошибки, т.е. функции вида F=(y(x0)-y0)^2 + ... + (y(xn)-yn)^2

Реализация метода - либо пакет levmar (c/c++), либо из ALGLIB - http://alglib.sources.ru/optimization/levenbergmarquardt.php

Добавлено 11.10.09, 17:04
Добавление 1: да, минимизация ведется по A,B, xo, yo.

Добавление 2: в такой постановке задачи у тебя минимизируется ошибка аппроксимации (т.е. ошибка вычисления функции y), а не ошибка нахождения неизвестных параметров. Однако на практике можно исходить из того, что малая ошибка аппроксимации обозначает малую ошибку в параметрах функции.

(спасибо, пока думаю)

Добавлено 11.10.09, 19:13
Тут нужна некоторая помощь.

У нас выборка x(i), y(i)
Как я понимаю, в качестве f(i) берется (Y - A*sin(B*x(i) + X0) + Y0)^2
Где Y,A,X0,Y0 - переменные от которых зависит функция.
Я взял реализацию на C# и применял ее вот так:

        Dim s As New alglib.minlm.lmstate
        Dim sol() As Double = {0.25, 1591.549431, 0, -150} ' {0.1, 1000, 0, -120} '  'Примерно подобранные начальные значения
        'ReDim sol(results.Count)
        alglib.minlm.minlmfj(4, results.Count, sol, 0.0000001, 0.0000001, 10000, s)' Тут числа взяты наугад
 
        While alglib.minlm.minlmiteration(s)
            Dim A As Double = s.x(0)
            Dim B As Double = s.x(1)
            Dim X0 As Double = s.x(2)
            Dim Y0 As Double = s.x(3)
            If s.needf Then
                Dim _F As Double = 0
                For Each R As PointF In results
                    _F += (R.Y - (A * Math.Sin(B * R.X + X0) + Y0)) ^ 2
                Next
                s.f = _F
                Debug.Print(_F.ToString) 'ТУТ вывожу значения вычисляемой функции
            End If
            If s.needfij Then
                For i As Integer = 0 To results.Count - 1
                    Dim r As PointF = results(i)
                    s.fi(i) = r.Y - (A * Math.Sin(B * r.X + X0) + Y0)
                    s.j(i, 0) = Math.Sin(B * r.X + X0) REM dF/dA
                    s.j(i, 1) = A * B * Math.Cos(B * r.X + X0) REM dF/dB
                    s.j(i, 2) = A * Math.Cos(B * r.X + X0) REM dF/dX0
                    s.j(i, 3) = 1 REM dF/dY0
                Next
            End If
            If s.needfg Or s.needfgh Then
                Throw New Exception()
            End If
        End While
        Dim rep As alglib.minlm.lmreport
        Dim X(4) As Double
        alglib.minlm.minlmresults(s, X, rep)
        Dim z As String = "{" & X(0) & ", " & X(1) & ", " & X(2) & ", " & X(3) & "}"
        MsgBox(z)

Использую вот такой код.


Теперь пробую прогнать ее на одной и той же выборке, при разных значениях стартовых параметров
При стартовых параметрах {0.25, 1591.549431, 0, -150} (это довольно близко к правде) алгоритм дает ответ
{0.250000043588598, 1591.54946567717, 2.17883057160585E-08, -150.00006958343}
Что похоже на правду.
При этом последовательные значения функции вот такие:

Цитата

175285943.807378
175279093.429274
174998443.565119
155484743.497379
47898495.6142789
43850845.497696

Если же взять другие начальные значения, например {0.1, 1000, 0, -120} что уже довольно далеко от истины, получается ответ
{0.0999999859957626, 1000.00000574627, 5.74626982824354E-09, -120.000062631127}
А это совсем неверно
При этом вычисленные значения функции следующие:

Цитата

142013745.650736
142008195.060041
141780793.005068
125971771.737748
38809425.2111686
35530193.4296744

Вообще какие-то громадные числа получаются. Причина завершения в обоих случаях "relative step is no more than EpsX"

Ну, насколько я помню этот метод, он ищет локальный минимум, а если их несколько (а в реальных многопараметрических задачах - обычно очень много) - в какой он попадет, как раз зависит от начальных условий.

Гм. Ну а как тогда все таки решать исходную задачу? Нужно найти, получается, глобальный минимум, причем в нем f(x) = 0
И посмотрите плз, правильно ли посчитан якобиан

А нельзя ли для нахождения параметров синуса применить спектральный анализ?

Что?)

Э-э-э... может я не понял суть задачи .. в более понятных мне терминах радиотехнических, она звучит так: необходимо определить амплитуду, частоту, начальную фазу и постоянную составляющую синусоидального сигнала, заданного дискретными значениями. Только значения y должны быть "взяты" через равные промежутки времени, а в Вашей формуле вместо времени фигурирует некий x, и о нем ничего неизвестно. Если растояния между соседними значениями x одинаковы, можно вычислить спектр сигнала.

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

Цитата

А нельзя ли для нахождения параметров синуса применить спектральный анализ?

Да, совершенно верно, я о том же подумал.
Левенберг-Марквардт в данном случае - это чересчур.

ANDLL, неправильно считается Якобиан.

Цитата ANDLL @ 11.10.09, 18:49

s.j(i, 0) = Math.Sin(B * r.X + X0) REM dF/dA
s.j(i, 1) = A * B * Math.Cos(B * r.X + X0) REM dF/dB
s.j(i, 2) = A * Math.Cos(B * r.X + X0) REM dF/dX0
s.j(i, 3) = 1 REM dF/dY0

Во-первых, в каждой строчке потерян минус. Во-вторых, в строке "s.j(i,1)=" есть лишний множитель B.

Гм.
Только наверное
s.j(i, 1) = - A * r.X * Math.Cos(B * r.X + X0) REM dF/dB

Цитата 0797293 @ 12.10.09, 10:21

Цитата
А нельзя ли для нахождения параметров синуса применить спектральный анализ?

Да, совершенно верно, я о том же подумал.
Левенберг-Марквардт в данном случае - это чересчур.

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


P.S. Вообще, наверное надо добавить в алгоритм какой-то код для проверки точности вычисления Якобиана/градиента/Гессиана.

Добавлено 12.10.09, 10:28

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

Гм.
Только наверное
s.j(i, 1) = - A * r.X * Math.Cos(B * r.X + X0) REM dF/dB

А, ну да, конечно же вы правы я неясно выразился.

Заменил код.
Вот так правильный якобиан?

                    s.j(i, 0) = -Math.Sin(B * r.X + X0) REM dF/dA
                    s.j(i, 1) = -A * r.X * Math.Cos(B * r.X + X0) REM dF/dB
                    s.j(i, 2) = -A * Math.Cos(B * r.X + X0) REM dF/dX0
                    s.j(i, 3) = -1 REM dF/dY0

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

Цитата

Во-первых, методы, основанные на БПФ, позволяют работать только с точками на равномерной сетке

Вот в этом и была неясность. Раз сетка неравномерная, со спектром я промахнулся.

К слову сказать, точки и правда лежат на равномерной сетке(ну или по крайней мере можно выкинув небольшое число точек ее получить) , чето я сразу не подумал что это важно.

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

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

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

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

Цитата

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

Тогда,

Цитата

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

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

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

А ссылку можно на спектральный метод?
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
А есть алгоритмы, которые решают задачу не поиска минимума функции, а поиска нуля(или близкого к нему) функций такого вида?

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

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

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

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

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

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

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

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

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

Цитата

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

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

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

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

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

Цитата

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
Что даже по данным видно - совсем не амплитуда.
Почему?
Код показывать?

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

Вот
Запускать проект 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)

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

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

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

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

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

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

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

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

ANDLL
Если мы увеличиваем частоту дискретизации нам это ничего не даст. Но если мы соберем больше по времени точек то точность повысится.
Я увеличил в 2 раза число точек у меня ошибка у пала с 0.5 до 0.125 Не как не пойму куда точность девается.

Люди, вы меня запутали. Объясните, что такое за формула:
B:=1 / (pi/n*MaxI / 50)
Я пытаюсь следить за ходом событий, и не почень получается.
Вот Pavia приводит формулу для В:
B:=2*pi/n*MaxI / 50?
50 - это шаг сетки по Х? Максимальная круговая частота в этом случае pi/50. Так? И значения В должны быть значительно меньше 1.
И тут же следует другая формула. И В за 1000.

Цитата

Если мы увеличиваем частоту дискретизации нам это ничего не даст. Но если мы соберем больше по времени точек то точность повысится.

Если максимум расположен в районе 5-й гармоники, может, уменьшить частоту дискретизации? В этом случае точность тоже возрастет?

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

Нашел где точность теряется. Надо убирать постоянную составляющую. А то она дает SinC на спекторе.
Убирать как m:=mean(p); sub(p,m)

Сразу точность возросла до максимума. Осалось понять что у меня с фазой.

Добавлено 12.10.09, 21:38
Есть подозрение где еще теряется точность. Если сигнал обрублен не целое число гармоник укладывается. Надо искать переход через ноль.

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

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

Это неоднозначный вопрос - половина людей делит на N при прямом, вторая половина - при обратном преобразовании. Как ни делай, кто-то будет недоволен Я решил сделать так, как реализовано в MATLAB - чтобы минимизировать количество непоняток у пользователей.

Аппроксимируется с детерминацией 0,999988927 и остаточной дисперсией 0,0005146545
методом множественной регрессии с линеаризацией синуса....

На графике синяя линия - результат аппроксимации...(в прикрепленном файле)

y=Asin(BX+X0)+Y0 имеет гармонику низшего порядка...а имеет ли она отношение к физической сути описываемого процесса...
Прикреплённый файлКнига1.xlsx (72.61 Кбайт, скачиваний: 152)

repz
А можно все таки уточнить, откуда именно появилась синяя линия?
Интересует все таки метод получения, а не результат

Цитата ANDLL @ 13.10.09, 11:19

repz
А можно все таки уточнить, откуда именно появилась синяя линия?
Интересует все таки метод получения, а не результат

Здесь принципиально определиться что есть по "ИКС"
Если это аргумент функции, который является некоторой физической величиной, и шкала не регулярна - это одна тема (хотя по представленным данным она вроде как регулярна)

Ежели это шкала категории (номер события, шкала времени) - простая линейная регрессия (заявленной функции) с предиктором sin(BX+X0), где B и X0 параметры категории шкалы (просто задаются)

Спасибо всем кто участвовал, общими усилиями вопрос решился.