Сплайн-интерполяция потактно

Hello Vitali.

VN> Имеется некий девайс, который с частотой дисретизации 1 кГц VN> соответственно каждую миллисекунду производит измерение напряжения VN> сети на одном из своих входов. Hапряжение это почти синусоидальное VN> (промышленная электросеть), частота примерно 50 Гц (вообще-то может VN> быть 40-70 Гц). Измеренные отсчеты складываются в кольцевой буфер на VN> 1000 значений (т.е. за последнюю секунду). Задача состоит в вычислении VN> на каждом такте нескольких значений, сдвинутых на 1/6 периода, 1/3, VN> 2/3, 5/6 и т.д. (нужно сделать трехфазное напряжение из однофазного с VN> заданным сдвигом, причем, именно повторить форму.) Частоту напряжения VN> измеряю довольно точно (сотые доли Гц).

VN> Для справки: памяти много, поэтому можно хранить несколько массивов, VN> если идти по пути оптимизации "обмен памяти на время". Увеличить VN> частоту дискретизации не могу.

А уменьшить до 900Hz можешь?

Alexey

Hello Vitali.

VN> Hе вижу смысла. А это зачем? Частота-то может быть произвольной в VN> диапазоне 40-70 Гц.

Тогда да, смысла нет.

А если мерять частоту, и подстраивать частоту выборок на кратную частоте сигнала?

Alexey

Привет!

"Alexey Boyko"

Не вижу смысла. А это зачем? Частота-то может быть произвольной в диапазоне

С уважением,

Виталий Насенник

Mon Jun 28 2004 11:15, Vitali Nassennik wrote to All:

VN> Имеется некий девайс, который с частотой дисретизации 1 кГц VN> соответственно каждую миллисекунду производит измерение напряжения сети VN> на одном из своих входов. VN> Задача состоит в вычислении на каждом такте нескольких значений, VN> сдвинутых VN> на 1/6 периода, 1/3, 2/3, 5/6 и т.д. (нужно сделать трехфазное напряжение VN> из однофазного с заданным сдвигом, причем, именно повторить форму.)

Можно подойти к задаче немного другим способом: Если частота известна, то сделать ЛЗ на произвольное время порядка 1/4 периода, и получить все нужные тебе фазы складывая векторы со входа и с выхода ЛЗ. Если частота неизвестна, то построить фильтр - фазовращатель на 90 и получить нужные фазы, складывая прямой и повернутый сигнал.

VLV

"There is no business other then show business " (c)

Hello, Vitali!

Пон Июн 28 2004, Vitali Nassennik писал к All по поводу "Сплайн-интерполяция потактно." VN> Имеется некий девайс, который с частотой дисретизации 1 кГц VN> соответственно каждую миллисекунду производит измерение напряжения VN> сети на одном из своих входов. Hапряжение это почти синусоидальное VN> (промышленная электросеть), частота примерно 50 Гц (вообще-то может VN> быть 40-70 Гц). Измеренные отсчеты складываются в кольцевой буфер на VN> 1000 значений (т.е. за последнюю секунду). Разрядность значений? (диапазон?) Формат их представления? Процессор какой? VN> Задача состоит в вычислении на каждом такте нескольких значений, VN> сдвинутых на 1/6 периода, 1/3, 2/3, 5/6 и т.д. (нужно сделать VN> трехфазное напряжение из однофазного с заданным сдвигом, причем, VN> именно повторить форму.) Частоту напряжения измеряю довольно точно VN> (сотые доли Гц). Короче, требуется вычислить значение, задержанное VN> относительно текущего момента на некое t (не целое в тактах). Вычислить 1 раз ? Или для каждой точки с шагом 1/6 (6000 раз) ? VN> Попробовал линейно интерполировать - получается не очень красиво, VN> порядка 2% погрешности. Зато дешего надежно и практично. VN> Попробовал кубическую сплайн-интерполяцию - очень понравилось, VN> однако, вычислять сплайн методом прогонки по нескольким точкам каждый VN> раз - мягко говоря, неоптимально. Если памяти много а диапазон невелик можно оптимизировать сплайн в таблицу, опять-же не очень ясно какие ресурсы есть для таких вычислений и сколько вычислений надо. VN> Вопрос состоит в том, каким образом можно было бы вычислять VN> коэффициенты для сплайна по одному за такт? Или кто может подсказать VN> какое-то другое решение этой задачи? VN> Увеличить частоту дискретизации не могу. VN> Виталий Hасенник WBR! Maxim Polyanskiy.

Привет!

"Maxim Polyanskiy"

какой?

Процессор Geode (навроде Pentium MMX). Формат представления - float.

На каждый такт нужно вычислить 5 значений, отстающих через 1/6 периода, т.е. с задержкой на 1/6, 1/3, 1/2, 2/3, 5/6 периода.

Оно понятно, что дешево, но 2% это слишком много. Не пойдет. Хочу сплайн.

таблицу,

Не понял, это как??? Сплайн интерполяция это и есть таблица для вычисления в интервале, опираясь на значения в сетке! Таблица для вычисления таблицы??? Хм...

i=trunc((x-x0)/h) p=(x-x0-i*h)/h q=1-p

F(x)=f[i]*q+f[i+1]*p+s[i]*(q-q*q*q)+s[i+1]*(p-p*p*p) ну или после оптимизации F(x)=f[i]*q+f[i+1]*p+p*q*(s[i]*(1+q)+s[i+1]*(1+p))

f[] - массив со значениями функции в точках с шагом h s[] - массив со значениями коэффициентов

Короче, на каждый такт у меня в массив справа добавляется один отсчет, который мне понадобится потом еще 5 раз за период. Соответственно, для этого отсчета мне нужно вычислить коэффициент для сплайна. Мне неохота соответственно заново вычислять все коэффициенты, если N-1 коэффициент слева от него уже известен. Для гладкой сшивки можно использовать тот факт, что самая близкая точка, которая мне понадобится для вычисления сигнала, задержанного на 1/6 периода, отстоит по крайней мере на 2 точки от текущей, т.е. я смело могу вычислять коэффициент для этой точки зная, что еще по три точки вправо и влево у меня есть, но коэффициенты я для точек правее нее еще не считаю, поскольку мне неизвестно значение производных в самой правой свежедобавленной точке, однако этого вполне достаточно для вычисления первой и второй производной в точке, для которой мне нужно вычислить коэффициент сплайна.

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

С уважением,

Виталий Насенник

Hello, Vitali!

Сpд Июн 30 2004, Vitali Nassennik писал к Maxim Polyanskiy по поводу "Re: Сплайн-интерполяция потактно." VN> "Maxim Polyanskiy" VN> Процессор Geode (навроде Pentium MMX). Формат представления - float. Hеужели прям с АЦП и сразу float? VN> Hа каждый такт нужно вычислить 5 значений, отстающих через 1/6 VN> периода, т.е. с задержкой на 1/6, 1/3, 1/2, 2/3, 5/6 периода. Понял. т.е. 200мкс на вычисление. >> VN> Попробовал линейно интерполировать - получается не очень >> VN> красиво, порядка 2% погрешности. >> Зато дешего надежно и практично. VN> Оно понятно, что дешево, но 2% это слишком много. Hе пойдет. Хочу VN> сплайн. Помоему это реально написать в лоб на асме либо через FPU либо MMX коммандами. VN> С математической точки зрения все понятно, нужно брать и выводить VN> формулы, однако, жутко не хватает времени, вот я и спрашиваю в эхе, VN> вдруг кто-то уже с такой постановкой проблемы сталкивался и таким VN> методом задачу решил? Я с такой производительностью не решал. Вряд-ли смогу помочь. VN> Виталий Hасенник WBR! Maxim Polyanskiy.

Hi!

In a message of 30 Jun 04 Vitali Nassennik wrote to Maxim Polyanskiy:

VN> i=trunc((x-x0)/h) VN> p=(x-x0-i*h)/h VN> q=1-p

VN> F(x)=f[i]*q+f[i+1]*p+s[i]*(q-q*q*q)+s[i+1]*(p-p*p*p) VN> ну или после оптимизации VN> F(x)=f[i]*q+f[i+1]*p+p*q*(s[i]*(1+q)+s[i+1]*(1+p))

VN> f[] - массив со значениями функции в точках с шагом h VN> s[] - массив со значениями коэффициентов

Эээ, хм... А почему бы не использовать кусочную интерполяцию? Типа такой вот:

--[cut]--

;interpolates between Y1 & Y2 ;uses formula: Y=At^3+Bt^2+Ct+D, t=0..1

; A=(-Y0+3*Y1-3*Y2+Y3)/2 ; B=(2*Y0-5*Y1+4*Y2-Y3)/2 ; C=(-Y0+Y2)/2 ; D=Y1

--[cut]--

Y0-Y3 - 4 последовательных отсчёта на равных расстояниях друг от друга. Сплайн проводится между Y1 и Y2, потом между Y2 и Y3 (берутся уже точки Y1-Y5) и так далее. t = (X-X1)/(X2-X1), причём расстояния Xn+1 - Xn постоянны.

Сплайн получается с непрерывной первой производной.

Эти формулы для равноотстоящих отсчётов, если надо, могу дать и для отсчётов в произвольных точках.

Bye...

Привет!

"Vadik Akimoff"

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

Вот насчитываю три первые производные в узлах:

h1:=1/(12*h); h2:=1/(24*h*h); h3:=1/(12*h*h*h); for i:=2 to N-2 do begin f1[i]:=(8*(f[i+1]-f[i-1])+f[i-2]-f[i+2])*h1; f2[i]:=(-30*f[i]+16*(f[i+1]+f[i-1])-f[i-2]-f[i+2])*h2; f3[i]:=(2*(f[i-1]-f[i+1])+f[i+2]-f[i-2])*h3; end;

Это я насчитываю по 5 точкам первые три производные для каждой узловой точки. h - шаг между узлами. xl - расстояние до узла слева, xr - расстояние до узла справа. i - индекс массива, соответствующий узлу слева.

i:=trunc((x-x0)/h); xl:=x-x0-i*h; xr:=xl-h; p:=xl/h; q:=1.0-p; L:=f[i]+(f1[i]+(f2[i]+f3[i]*xl)*xl)*xl; R:=f[i+1]+(f1[i+1]+(f2[i+1]+f3[i]*xr)*xr)*xr; Spl:=L*q+R*p;

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

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

L:=f[i]+(f1[i]+f2[i]*xl)*xl; R:=f[i+1]+(f1[i+1]+f2[i+1]*xr)*xr;

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

С уважением,

Виталий Насенник