[serjj1333 0]

И, если вы прочитали всю приведенную ссылку, в ней, опять же, подтверждают мои слова о разрывности дифференцируемого вами арктангенса

И прочитал и в своё время сделал в железке, не поверите :rolleyes:

Его там исключают потому что далее делается дифференцирование и можно вместо d/dt(arg(y(t))) рассчитывать конечное выражение и собрать эквивалентую схему. Никто не запрещает сделать схему "в лоб" с использованием арктангенса + unwrap и всё будет работать (проверено на практике, работает, alas!).

Да и где ж разрыв?

 

Поэтому, в реальной обработке, арктангенс старательно исключают.

Приехали. Преобразование координат и демодуляцию PSK-N созвездий наверное тоже без арктангеса делают, ага.

 

Расскажите мне, как определить фазу сигнала по одному отсчету АЦП? Ну ладно, учитывая нарисованные вами производные - пусть будет целых три последовательных отсчета с интервалом 1 мксек: -1В; 0В; +1В.

Ну, хотя бы, любимую "мгновенную" частоту определите?

А мне не нужно работать на уровне АЦП. И ТС'у никто не предлагал работать на уровне АЦП и по одному отчёту. Есть сигнал после гетеродина и ФНЧ, и он комплексный. Для него понятие фазы определено строго математически и применимо к каждому символу. А раз определено понятие фазы, то значит можно оценить мгновенную частоту. Точность такой оценки будет зависеть от сигнал-шума. Процесс демодуляции ЧМ сигналов привёл как простой и распространённый пример.

 

Что ещё Вы хотите тут опровергать?

Изменено 14 сентября, 2015 пользователем serjj

[mcheb 0]

Прошу прощения за нескромный вопрос. А в каком месте atan(x) разрывается?

Наверное, подразумевалась неоднозначность.

[blackfin 16]

Как я уже говорил ранее, понятие "мгновенная частота" в метрологии отсутствует. Оно есть только в головах пишущих формулы теоретиков.

Можно вполне строго доказать невозможность измерения "мгновенной частоты" на практике.

 

Для этого нужно воспользоваться формулой для вычисления частоты по трем значениям гармонической функции S(t) из соседней темы:

 

f = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

 

Это абсолютно точная формула, позволяющая получить точное значение частоты "f" при условии,

 

что нам известны точные значения гармонической функции S(t) в моменты времени t₁,t₂ и t₃:

 

S₁ == S(t₁),

S₂ == S(t₂),

S₃ == S(t₃),

 

где:

 

Δt == t₃ - t₂ = t₂ - t₁,

 

и при этом:

 

f*(t₃ - t₁) < 1.

 

Чтобы строго доказать невозможность измерения "мгновенной частоты" на практике, необходимо вычислить полный дифференциал функции трех переменных:

 

f(S₁,S₂,S₃) = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

 

А именно:

 

df == (δf/δS₁)*dS₁ + (δf/δS₂)*dS₂ + (δf/δS₃)*dS₃,

 

где частные производные функции "f" равны:

 

δf/δS₁ = -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[1/(2*S₂)],

 

δf/δS₂ = +[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[(S₁+S₃)/(2*S₂²)].

 

δf/δS₃ = -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[1/(2*S₂)],

 

Соответственно, полный дифференциал равен:

 

df == -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[(dS₁+dS₃-(S₁+S₃)*dS₂/S₂)/(2*S₂)].

 

Теперь мы можем найти погрешность измерения частоты "σ_F" связанную с погрешностью измерения сигнала "σ_S":

 

σ_F = σ_S*[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*sqrt[2+(S₁+S₃)²/S₂²]/[2*abs(S₂)],

 

или:

 

σ_F = σ_S*[1/[pi*(t₃ - t₁)]]*[1/sqrt[(2*S₂)²-(S₁ + S₃)²]]*sqrt[2+(S₁+S₃)²/S₂²].

 

Знаменатель в этой формуле примечателен тем, что при стремлении t₃ к t₁ он стремится к нулю как const*(t₃-t₁)²

 

Как следствие, при стремлении t₃ к t₁ для любой погрешности "σ_S" измерения амплитуды сигнала S(t), погрешность измерения "мгновенной частоты" "σ_F" устремляется к бесконечности как σ_S/(t₃-t₁)²..

[rudy_b 1]

Нет ничего проще: На основании именно этих полученных от АЦП отсчетов можно смело сказать, что на его входе присутствует синусоидальный сигнал с частотой 250 kHz, c амплитудой 1В. Центральному отсчету соответствует фаза 0. Мы считаем, что амплитуда сигнала нормирована под полный "раскрыв" АЦП +/-1В.

 

Смелый вы человек. Но с чего вы взяли, что это синус? Боженька на ушко шепнул? А может это линейно нарастающее напряжение? Или пила, трапеция и т.д. с произвольным периодом?

Или вообще непериодический сигнал?

 

Но, даже если предположить чистый синус, то с чего вы решили, что амплитуда сигнала равна 1В?

Вот набор пар амплитуда/частота который дает именно эти отсчеты. Ваш вариант - это только один из них.

Амплитуда ----------- Частота --------- Напряжение при t=1 мксек

1.000000E+00 -- 2.500000E+05 1 -- Ваш вариант

1.500000E+00 -- 1.161398E+05 1

2.250000E+00 -- 7.329944E+04 1

3.375000E+00 -- 4.787579E+04 1

5.062500E+00 -- 3.164613E+04 1

7.593750E+00 -- 2.101973E+04 1

1.139063E+01 -- 1.399046E+04 1

1.708594E+01 -- 9.320293E+03 1

2.562891E+01 -- 6.211555E+03 1

3.844336E+01 -- 4.140452E+03 1

5.766504E+01 -- 2.760129E+03 1

8.649756E+01 -- 1.840034E+03 1

и далее ...

 

А если взять частоты более 250 кГц (не совсем корректно, но чего только в жизни не бывает) - то снова получим неограниченный набор пар частота/амплитуда, котрые дадут те же отсчеты.

 

Если вы собираетесь заниматься гаданием на кофейной гуще, то вы не туда попали, эта тема связана с метрологией.

 

А мне не нужно работать на уровне АЦП. И ТС'у никто не предлагал работать на уровне АЦП и по одному отчёту. Есть сигнал после гетеродина и ФНЧ, и он комплексный.

Так как же быть с вашим определением "мгновенной" частоты? Опять увиливаете от ответа? Или уже убедились, что это чушь? Тогда так и скажите и закончим на этом эту бесплодную дискуссию.

 

Прошу прощения за нескромный вопрос. А в каком месте atan(x) разрывается?

Действительно не совсем корректно выразился. Разрывность появляется при попытке определить именно фазовый угол (что, при грамотной обработке, практически не требуется) по соотношению Re и Im компонент применяя функцию арктангенса. Нужна, как минимум, atan2, но и с ней получаемый фазовый угол будет не линейно нарастающим (если есть постоянная разность частот), а скакать при переходе через +/-пи.

 

Наверное, подразумевалась неоднозначность.

Это как минимум.

[serjj1333 0]

f = [1/[pi*(t3 - t1)]]*arccos[(S1 + S3)/(2*S2)].

...

Как следствие, при стремлении t3 к t1 для любой погрешности "σS" измерения амплитуды сигнала S(t), погрешность измерения "мгновенной частоты" "σF" устремляется к бесконечности как σS/(t3-t1)2..

А тыкните носом в раздел анализа и теории оценок, где написано, что надо так считать. Просто интересно :rolleyes:

В моём понимании f = [1/[pi*(t3 - t1)]]*arccos[(S1 + S3)/(2*S2)] при t3-t1->0 примет вид acos(1)/0, что даёт ноль в числителе и ноль в знаменателе, нормальное поведение для производной.

Кстати ваш вывод будет справедлив только для гармонического сигнала, но всё равно любопытно.

 

С другой стороны для идеального (нет шума) действительного гармонического сигнала можно написать такой простой скрипт:

clear all;
Ns = 1e4;
fs = 100e5;
Ts = 1/fs;
t = 0:Ts:(Ns-1)*Ts;
f0 = 1333.756;
x = sin(2*pi*f0*t);
x1 = sqrt(1 - x.^2);
ph = angle(x1 + 1i*x);
index = randi([2 Ns]);
fcalc = abs(ph(index)-ph(index-1))/(Ts * 2 * pi);
str = sprintf('real frequency = %E -- estimated frequency = %E', f0, fcalc);
disp(str);

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

Увеличивая частоту дискретизации, мы сокращаем время между двумя точками, с другой стороны разница напряжений в этих точках тоже уменьшается, т.о. мы приближаемся к определению мгновенной частоты. Если шума нет, то мгновенная частота будет определяться с точностью до fs, если добавить шум, то погрешность, которую он внесёт, превысит погрешность от fs (при fs >> искомой частоты). Погрешность такого эстиматора уже будет оцениваться с помощью CRB.

 

Так как же быть с вашим определением "мгновенной" частоты? Опять увиливаете от ответа? Или уже убедились, что это чушь? Тогда так и скажите и закончим на этом эту бесплодную дискуссию.

Увиливаете тут только вы. Если вы хотите, чтобы ваши агрументы были приняты, опишите это математически или, если удобнее скриптом. Мы не на уроке литературы. А за одним посмотрите, чем отличается неоднозначность обратных тригонометрических функций от понятия разрыва и непрерывности. Это довольно познавательно. :rolleyes:

Изменено 15 сентября, 2015 пользователем serjj

[FatRobot 0]

Ваши рассуждения довольно любопытны для человека, который декларирует, что хоть сколько-то занимается метрологическим оборудованием. Сначала вы ставите условие и просите получить по данным из условия оценку частоты. После того, как оценка частоты получена, вы дополняете условия, и оценка оказывается не верной и/или не может быть получена. Знаменитый принцип ведения дискуссии: "так да не так".

 

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

 

А второе: любой оценщик, как использующий 3 отсчета, так и 333 отсчета, можно обмануть и придумать бесконечное множество "а если", которые заставят его показывать не то, что вы подразумевали. Это не повод опускать руки. Нужно ограничить область применения доп. условиями, предоставить больше априорной информации о входном воздействии к началу измерений.

 

Выводы простые в общем-то:

Можно ли оценивать частоту сигнала по 3м отсчетам? Можно, если имеется другая необходимая априорная информация или предположения о сигнале.

Можно ли "обмануть" любой, наперед заданный оценщик? Можно, если выйти из базиса гипотез об оцениваемом сигнале.

 

А существование понятия "мгновенная частота".. Угловое положение вращающегося диска можно измерить. Время тоже.

 

Смелый вы человек. Но с чего вы взяли, что это синус? Боженька на ушко шепнул? А может это линейно нарастающее напряжение? Или пила, трапеция и т.д. с произвольным периодом?

Или вообще непериодический сигнал?

 

Но, даже если предположить чистый синус, то с чего вы решили, что амплитуда сигнала равна 1В?

Вот набор пар амплитуда/частота который дает именно эти отсчеты. Ваш вариант - это только один из них.

 

[...]

 

А если взять частоты более 250 кГц (не совсем корректно, но чего только в жизни не бывает) - то снова получим неограниченный набор пар частота/амплитуда, котрые дадут те же отсчеты.

 

Если вы собираетесь заниматься гаданием на кофейной гуще, то вы не туда попали, эта тема связана с метрологией.

[blackfin 16]

А тыкните носом в раздел анализа и теории оценок, где написано, что надо так считать. Просто интересно :rolleyes:

Читайте про преобразования случайной величины, Якобиан и проч.

 

Или: Г.Корн Т.Корн, стр. 561, Гл. 18.5-4, "Функции от случайных величин. Замена переменных."

 

Ну и имейте ввиду, что случайная величина "f" есть однозначная функция трех независимых случайных величин: f = f(S1,S2,S3)..

[thermit 1]

С другой стороны для идеального (нет шума) действительного гармонического сигнала можно написать такой простой скрипт:

Код

clear all;

Ns = 1e4;

fs = 100e5;

Ts = 1/fs;

t = 0:Ts:(Ns-1)*Ts;

f0 = 1333.756;

x = sin(2*pi*f0*t);

x1 = sqrt(1 - x.^2);

ph = angle(x1 + 1i*x);

index = randi([2 Ns]);

fcalc = abs(ph(index)-ph(index-1))/(Ts * 2 * pi);

str = sprintf('real frequency = %E -- estimated frequency = %E', f0, fcalc);

disp(str);

 

Для неопределенной амплитуды скрипт, естественно мертвый. 3 точки нужно минимум.

[serjj1333 0]

Для неопределенной амплитуды скрипт, естественно мертвый. 3 точки нужно минимум.

Да, точно. Ну и на практике никто так не делает (x1 = sqrt(1 - x.^2)), если нужно получить аналитический сигнал, то ставят Гильберта. Там никакой зависимости от амплитуды + класс входных сигналов не ограничен только гармоническими. Ну или гетеродин + ФНЧ - если диапазон измерения априори задан.

Изменено 15 сентября, 2015 пользователем serjj

[Pridnya 0]

Можно вполне строго доказать невозможность измерения "мгновенной частоты" на практике.

 

Для этого нужно воспользоваться формулой для вычисления частоты по трем значениям гармонической функции S(t) из соседней темы:

 

f = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

 

Это абсолютно точная формула, позволяющая получить точное значение частоты "f" при условии,

 

что нам известны точные значения гармонической функции S(t) в моменты времени t₁,t₂ и t₃:

 

S₁ == S(t₁),

S₂ == S(t₂),

S₃ == S(t₃),

 

где:

 

Δt == t₃ - t₂ = t₂ - t₁,

 

и при этом:

 

f*(t₃ - t₁) < 1.

 

Единственная формула, которую сразу можно посчитать в C++ (IDE NetBeans, компилятор MinGW ). Два результата: для 16-ти разрядного АЦП и для идеального АЦП (см. результат в конце C++ кода). Для хорошего АЦП метод точный, т.е. если фильтрануть основную гармонику, то будет хорошо.

 

Другие способы (MUSIC, MLE - Метод максимального правдоподобия...), по крайней мере, на первый взгляд очень сложные. Какие-то для матлаба, какие-то в общем виде многоэтажные формулы. Кто такие формулы пишет? Шифротекст для потомков. Неужели один я такой тупой, что не понимаю.

 

#include <cstdlib>
#include <iostream>   // std::cout
#include <math.h>

using namespace std;

#define SAMPLES 64
//short Data[SAMPLES] = {0}; // АЦП 16-ти разрядный.
float Data[SAMPLES] = {0}; // идеальный АЦП, но его же нет (сигнал 16-ти разрядного АЦП как-то нужно улучшить???).

int main(int argc, char** argv) {

// Создаем сигнал частотой 50 Гц (первая гармоника): косинус, амплитуда 1000, фаза 0.
// таблица содержит 64 выборки на интервале 0.02 сек (период промышленной частоты 50 Гц).    
for(int i=0;i<SAMPLES;i++)
    {Data[i]=1000*cos(2*M_PI*i*1/SAMPLES+2*M_PI*0/360.0F);} // 1000*sin(1wt+0).

// Три равноотсоящих точки выборки (первые три точки).
float S1 = Data[0];
cout << "S1 = " << S1 << endl;  // 1000
float S2 = Data[1];
cout << "S2 = " << S2 << endl;  // 995
float S3 = Data[2];
cout << "S3 = " << S3 << endl;  // 980
float dT = 1.0F/3200.0F;
cout << "dT = " << dT << endl;  // 0.0003125 сек (период между выборками).
float t1 = 0; 
cout << "t1 = " << t1 << endl;  // 0 сек.
float t3 = 2*dT;
cout << "t3 = " << t3 << endl;  // 0.000625 сек.

//f = [1/[pi*(t3 - t1)]]*acos[(S1 + S3)/(2*S2)] // Оригинальная формула.
float f = (1/(M_PI*(t3 - t1)))*acos((S1 + S3)/(2*S2)); // Заменил квадратные скобки на круглые.

cout << "f = " << f << endl;    // 51.078 Гц ??? для массива short Data[64]
cout << "f = " << f << endl;    // 50.0001 Гц для массива float Data[64]

return 0;

Изменено 15 сентября, 2015 пользователем Pridnya

[blackfin 16]

Единственная формула, которую сразу можно посчитать в C++.

"Сложные проблемы всегда имеют простые, легкие для понимания неправильные решения."

[serjj1333 0]

Другие способы (MUSIC, MLE - Метод максимального правдоподобия...), по крайней мере, на первый взгляд очень сложные. Какие-то для матлаба, какие-то в общем виде многоэтажные формулы. Кто такие формулы пишет? Шифротекст для потомков. Неужели один я такой тупой, что не понимаю.

Вы никак не улучшите ваш сигнал оцифрованный 16-битным АЦП. При таком соотношении символьной и искомой частот (3200/50) аргумент функции acos -> 1, а значит сама функция - к 0, т.е. некоторому малому числу. Далее вы делите на 2pi*dt, другое малое число. С вычислительной точки зрения мало устойчивая задача. Любое незначительное искажение сигнала внесёт значительное искажение в оценку. Пропустив сигнал через АЦП вы добавили шум квантования, который внёс ошибку в ваш эстиматор. А теперь представте - в реальности при оцифровке стараются сделать так, чтобы оцифрованный аддитивный шум (тепловой шум, шум канала, кароче говоря неустронимый физический шум) был много больше шума квантования - условие необходимое для сохранения SNR на том уровне, который был до оцифровки. И посчитайте, какую ошибку внесёт он. Будете неприятно удивлены, т.к. даже незначительный шум даст полностью бесмыссленные значения частот. И даже узкополосный полосовой фильтр на 50 Гц вряд ли существенно улучшит результат.

Так что вам придётся применять какой либо из "сложных" методов.

Изменено 15 сентября, 2015 пользователем serjj

[blackfin 16]

При таком соотношении символьной и искомой частот (3200/50) аргумент функции acos -> 1, а значит сама функция - к 0, т.е. некоторому малому числу.

Далее вы делите на 2pi*dt, другое малое число. С вычислительной точки зрения мало устойчивая задача.

Так Вы, оказывается, не поняли основной идеи предложенной формулы.. :rolleyes:

 

Я же там где-то раньше указывал, что для получения точного значения частоты должно выполняться условие: abs(S₁ + S₃) < abs(S₂).

 

Это требование фактически означает, что для вычисления "f" годятся не любые моменты времени t₁,t₂ и t₃, а только те,

 

в которых значения S₁ и S₃ находятся "вблизи" нуля и при этом значение S₂ "автоматически" оказывается вблизи максимума функции S(t).

 

Это примерно соответствует вычислению частоты "вручную", то есть когда мы просто измеряем моменты времени в которые гармонический сигнал S(t) пересекает ось абсцисс,

 

потом находим соответствующий этим моментам времени период гармоники "T" и уже затем вычисляем её частоту f = 1/T. Просто в предложенном способе не нужно точно вычислять моменты времени, когда S(t) = 0.

 

PS. А про шумы полностью согласен. Формула с arccos применима только для сферических функций в вакууме..

[Pridnya 0]

...

PS. А про шумы полностью согласен. Формула с arccos применима только для сферических функций в вакууме..

 

Формула сама интересная и, если сигнал будет чистый, то должна дать хороший результат.

Попробую сигнал с выхода АЦП обработать НЧ фильтром с КИХ и полосой 75Гц -3 дБ, 100 Гц -40 дБ, 64 коэффициента float.

А там можно разные способы попробовать. Что-то похожее советовал мне petrov.

Pridnya

 

Выделяйте комплексным полосовым КИХ фильтром интересующий диапазон, АЧХ в полосе пропускания должна быть как можно более равномерная, гармоники(в том числе на отрицательных частотах) должны быть задавлены как можно сильнее, затем вычисляете несколько бинов ДПФ, далее по статье считаете частоту:

 

http://electronix.ru/forum/index.php?s=&am...t&p=1141831

 

Советую модельку в симулинке сделать и погонять с различными параметрами.

Только вот не пойму, откуда бины возьмутся после фильтрации сигнала фильтром с частотой среза 75 Гц, вроде там должен быть один бин, т.е. 50 Гц. Или я что-то не понимаю. Или он имеет ввиду посчитать 1024 точечное ДПФ от сигнала на интервале 0,02 секунды и по тем бинам посчитать?

 

[blackfin 16]

А там можно разные способы попробовать. Что-то похожее советовал мне petrov.

Эту задачу на форуме решали уже десятки раз..

 

См. например:

Но вот скажем сигнал: синус частотой 40...60Гц. Требуемая точность =>0,01Гц.

Частота дискретизации 6400 (для примера).

 

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

Как я понимаю чем выборка больше тем лучше?

Увеличивая интервал обработки можно получить как-угодную точность каким-угодно методом. В этих задачах оптимальность означает - как быстро растёт точность при увеличении интервала.

Оптимальность это не только научные понты для ученых. Мы не можем реально значительно увеличивать интервал измерения, только в пределах стабильности самих параметров сигнала (амплитуды, частоты)

 

Сравните с оцениванием спектра и обратной квадратичной интерполяцией. Если интервал частот узкий, то сразу можно начинать с того, что называют ML-extension:

1. Умножить на центральную комплексную экспоненту (снести сигнал в 0) и отфильтровать ФНЧ. Соответственно прорежение.

2. Взять 5-7 сумм типа ДПФ в интервале [-7.5, 7.5]. Оценка энергии

3. Найти максимум энергии и по трём точкам построить параболу. Аргумент максимума - частота

 

Точность будет примерно такая-же.

И там ещё в окрестностях есть что почитать.. ;)