| |
|
  |
 Измерение частоты основной гармоники (50 Гц) с точностью 0.01 Гц |
|
|
|
|
 Nov 10 2015, 10:57
|
Частый гость
Группа: Свой
Сообщений: 142
Регистрация: 11-01-11
Из: Орел
Пользователь №: 62 159
|
Цитата(blackfin @ Nov 10 2015, 12:11)  А вы перед тем, как считать БПФ, сначала умножьте свою основную гармонику (50,00000 Гц 49,9092 Гц) на весовой множитель: g(t) = exp(-t 2/T 2), где: T = 0.1 [сек]. А потом сравните ошибку метода "с параболической интерполяцией" с ошибкой своего метода "по двум бинам".. PS. Формула для g(t) предполагает, что сигнал отличен от нуля на интервале времени: -Tсигн/2 ... +Tсигн/2. Извините, но мне совсем не понятно, как и зачем умножать. Единственное, что я понял, так это то, что вы предлагаете умножать на затухающую экспоненту (посмотрел g(t) для 32-х точек, через каждые 0,000625 секунд). Т.е. это не окно Гаусса...половина окна Гаусса что ли? Даже не понимаю, откуда возник такой способ. Это авторская методика? Где почитать про это вообще? Всю выборку что ли умножить на затухающую экспоненту? Код float t=0.0F;
float T=0.1F;
int SAMPLES = 32;
float g[SAMPLES];
for(int i=0;i<SAMPLES;i++)
{
g[i]= expf(-(t*t/T*T));
cout<<"g["<<i<<"]="<<g[i]<<endl;
t+=0.000625F;
}
/*
g[0]=1
g[1]=1
g[2]=1
g[3]=0.999999
g[4]=0.999998
g[5]=0.999998
g[6]=0.999996
g[7]=0.999995
g[8]=0.999994
g[9]=0.999992
g[10]=0.99999
g[11]=0.999988
g[12]=0.999986
g[13]=0.999983
g[14]=0.999981
g[15]=0.999978
g[16]=0.999975
g[17]=0.999972
g[18]=0.999968
g[19]=0.999965
g[20]=0.999961
g[21]=0.999957
g[22]=0.999953
g[23]=0.999948
g[24]=0.999944
g[25]=0.999939
g[26]=0.999934
g[27]=0.999929
g[28]=0.999923
g[29]=0.999918
g[30]=0.999912
g[31]=0.999906
*/ Цитата(blackfin @ Nov 10 2015, 12:38) Вы по ссылке то из поста #134 ходили? По ссылке ходил. Вообще ничего теперь не понимаю. Что это и как этим пользоваться.
Сообщение отредактировал Pridnya - Nov 10 2015, 11:04
|
|
|
|
|
|
|
|
Nov 10 2015, 11:31
|
Гуру
Группа: Свой
Сообщений: 3 106
Регистрация: 18-04-05
Пользователь №: 4 261
|
Цитата(Serg76 @ Nov 10 2015, 14:11) Окно необходимо только в случае, если в спектре кроме интересующей гармоники присутствуют еще какие либо мешающие воздействия, при наличии только "чистой" синусоиды лучшее окно - это прямоугольное, т. е. никакого. Вы забыли про отрицательные частоты в спектре вещественного сигнала: f - = -50 [Гц], f + = +50 [Гц].
|
|
|
|
|
|
|
|
Nov 10 2015, 11:41
|
Частый гость
Группа: Свой
Сообщений: 142
Регистрация: 11-01-11
Из: Орел
Пользователь №: 62 159
|
Цитата(blackfin @ Nov 10 2015, 14:31) Вы забыли про отрицательные частоты в спектре вещественного сигнала: f- = -50 [Гц], f+ = +50 [Гц]. А можно от абстрактного уровня перейти поближе к физическому. Если спектр состоит из одних амплитуд и отражает картину происходящего (максимум соответствует заданной частоте). Как мне на выходе БПФ ощутить наличие отрицательных частот и как это связано с предложенной вами формулой? И вообще, вписывается ли такая методика в какую-либо нацчную концепцию. А то я вижу, что вы называете "окном Гаусса" совсем другую функцию Цитата На практике, для окна Гаусса заданного формулой:
g(t) = exp(-t2/T2), , не похожую на функцию Гаусса. Выше половины частоты дискретизации все фильтруется антиалисинговым фильтром. В спктре присутствует только одна частота, например 52 Гц, отношение сигнал/шум пусть будет 40 дБ (это чтобы никто на влияние шумов не указывал). Сигнал частотой 52 Гц я задал с генератора сигналов. В спектре вижу максимум около 52-х Гц. Дальше считаю.
Сообщение отредактировал Pridnya - Nov 10 2015, 11:45
|
|
|
|
|
|
|
|
Nov 10 2015, 11:42
|
Гуру
Группа: Свой
Сообщений: 3 106
Регистрация: 18-04-05
Пользователь №: 4 261
|
Формулу поправьте: Код int SAMPLES = 32;
float g[SAMPLES];
for(int i=0;i<SAMPLES;i++)
{
g[i]= expf(-((0.16/0.1)*(0.16/0.1)*(i-SAMPLES/2+0.5)*(i-SAMPLES/2+0.5)/((SAMPLES/2)*(SAMPLES/2))));
cout<<"g["<<i<<"]="<<g[i]<<endl;
} См. Gaussian_window.
|
|
|
|
|
|
|
|
Nov 10 2015, 11:55
|
Частый гость
Группа: Свой
Сообщений: 142
Регистрация: 11-01-11
Из: Орел
Пользователь №: 62 159
|
Цитата(blackfin @ Nov 10 2015, 14:42) Формулу поправьте: Код int SAMPLES = 32;
float g[SAMPLES];
for(int i=0;i<SAMPLES;i++)
{
g[i]= expf(-((0.16/0.1)*(0.16/0.1)*(i-SAMPLES/2+0.5)*(i-SAMPLES/2+0.5)/(SAMPLES*SAMPLES)));
cout<<"g["<<i<<"]="<<g[i]<<endl;
} Формулу поправил. Теперь она чем-то на верхнюю часть функции Гаусса похожа ( дуга какая-то, на рисунке для 32-х точек, для примера). Теперь мне нужно всю выборку, все мои 1024 точек выборки умножить на такое окно?
Сообщение отредактировал Pridnya - Nov 10 2015, 11:58
Эскизы прикрепленных изображений
|
|
|
|
|
|
|
|
Nov 10 2015, 12:01
|
Гуру
Группа: Свой
Сообщений: 3 106
Регистрация: 18-04-05
Пользователь №: 4 261
|
Цитата(Pridnya @ Nov 10 2015, 14:55) Формулу поправил. Теперь она чем-то на верхнюю часть функции Гаусса похожа ( дуга какая-то, на рисунке для 32-х точек, для примера). Что-то не так поправили. На границе окна функция g(t) должна спадать до: g(0.16[сек]) = exp(-1.6*1.6) = exp(-2.56) = 0,0773. То есть, почти до нуля.. Цитата(Pridnya @ Nov 10 2015, 14:55) Теперь мне нужно всю выборку, все мои 1024 точек выборки умножить на такое окно? Ессно..
|
|
|
|
|
|
|
|
Nov 10 2015, 12:18
|
Частый гость
Группа: Свой
Сообщений: 142
Регистрация: 11-01-11
Из: Орел
Пользователь №: 62 159
|
Цитата(blackfin @ Nov 10 2015, 15:01) Что-то не так поправили.
На границе окна функция g(t) должна спадать до: g(0.16[сек]) = exp(-1.6*1.6) = exp(-2.56) = 0,0773.
То есть, почти до нуля..
Ессно.. Я не правил, с скопировал из вашего поста "как есть" Цитата g[i]= expf(-((i-SAMPLES/2+0.5)*(i-SAMPLES/2+0.5)/(SAMPLES*SAMPLES))); , затем вставил в свой код Посчитал коэффициенты этого окна, умножил на 1000 и через Open Office Calc и таблицу Столбец А - индекс, столбец B - значение, посторил диаграмму XY. Это чтобы иметь представление, какая там кривая за формулой скрывается. Теперь поправил (везде float-ы): g[i]= expf(-((0.16F/0.1F)*(0.16F/0.1F)*(i-SAMPLES/2.0F+0.5F)*(i-SAMPLES/2.0F+0.5F)/(SAMPLES*SAMPLES))); результат изменися, но все равно не похож на функцию Гаусса. PS: Что правил, что нет, а результат изменился не значительно. По двум бинам точнее в разы. Но после этого окна формула по двум бинам даёт большую ошибку, чем без окна вообще. Это окно плохо влияет на оценку частоты. Оно как-то пропорционально уменьшило все бины. Код fin = 52.00; // Гц
// Бины без окна.
x1 = 16; y1 = 456.366F;
x2 = 17; y2 = 793.953F;
x3 = 18; y3 = 205.923F;
// xmax = 16.9324
// fmax = 52.9136
// deltaf = -0.91362 Гц.
xmax = x2 + 0.25*(y3 - y1)/(2*y2 - y3 - y1); // = 16.0 + 0.25*(29.052 - 28.9915)/(2*999.485 - 29.052 - 28.9915) = 16.0 + 0.25*0.0605/1940.9265 = 16.00000779267.
fmax = xmax*3200.0/1024.0;
deltaf = fin - fmax;
cout <<"fin = " << fin << endl;
cout <<"xmax = " << xmax << endl;
cout <<"fmax = " << fmax << endl;
cout <<"deltaf = " << deltaf << endl;
// Бины после окна.
x1=16; y1=439.371F;
x2=17; y2=680.025F;
x3=18; y3=84.8844F;
//xmax = 16.894
//fmax = 52.7936
//deltaf = -0.793648
xmax = x2 + 0.25*(y3 - y1)/(2*y2 - y3 - y1); // = 16.0 + 0.25*(29.052 - 28.9915)/(2*999.485 - 29.052 - 28.9915) = 16.0 + 0.25*0.0605/1940.9265 = 16.00000779267.
fmax = xmax*3200.0/1024.0;
deltaf = fin - fmax;
cout <<"fin = " << fin << endl;
cout <<"xmax = " << xmax << endl;
cout <<"fmax = " << fmax << endl;
cout <<"deltaf = " << deltaf << endl;
// *************************************************************************** //
// Способ по двум бинам (хорошо считает, если левый бин не равен правому))
xmax = (16*y1 + 17*y2)/(y1+y2);
cout <<"fmax = (16*y1 + 17*y2)/(y1+y2)" << endl;
fmax = xmax*3200.0/1024.0;
deltaf = fin - fmax;
cout <<"fin = " << fin << endl;
cout <<"xmax = " << xmax << endl;
cout <<"fmax = " << fmax << endl;
cout <<"deltaf = " << deltaf << endl;
//fmax = (16*y1 + 17*y2)/(y1+y2)
//xmax = 16.6075
//fmax = 51.8984
//deltaf = 0.101585
Сообщение отредактировал Pridnya - Nov 10 2015, 12:26
Эскизы прикрепленных изображений
|
|
|
|
|
|
|
|
Nov 10 2015, 12:23
|
Знающий
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730
|
Код g[i]= expf(-((0.16F/0.1F)*(0.16F/0.1F)*(i-SAMPLES/2.0F+0.5F)*(i-SAMPLES/2.0F+0.5F)/(SAMPLES*SAMPLES))); SAMPLES -> SAMPLES/2
|
|
|
|
|
|
|
|
Nov 10 2015, 12:31
|
Частый гость
Группа: Свой
Сообщений: 142
Регистрация: 11-01-11
Из: Орел
Пользователь №: 62 159
|
Цитата(thermit @ Nov 10 2015, 15:23) Код g[i]= expf(-((0.16F/0.1F)*(0.16F/0.1F)*(i-SAMPLES/2.0F+0.5F)*(i-SAMPLES/2.0F+0.5F)/(SAMPLES*SAMPLES))); SAMPLES -> SAMPLES/2 Спасибо большое, что нашли ошибку в формуле blackfin-a. Если знаете где ошибка, то лучше пишите формулу правильно, а не размножайте неправильную формулу, а то так тяжело бывает понять эти намеки. К этому посту мы родили окно: g[i]= expf(-((0.16F/0.1F)*(0.16F/0.1F)*(i-SAMPLES/2.0F+0.5F)*(i-SAMPLES/2.0F+0.5F)/((SAMPLES/2.0F)*(SAMPLES/2.0F)))); PS: И даже после умножения выборки на это правильное окно результат улучшился незначительно. Цитата // БПФ без окна и параболическая аппроксимация.
fin = 52
xmax = 16.9324
fmax = 52.9136
deltaf = -0.91362
// БПФ с окном и параболическая аппроксимация.
fin = 52
xmax = 16.8536
fmax = 52.6674
deltaf = -0.667435
// БПФ и аппроксимация по двум бинам.
// В случае использования окна формула по двум бинам начинает ухудшать оценку, было 0.015625 Гц, стало 0.25 Гц..
fmax = (16*y1 + 17*y2)/(y1+y2)
fin = 52
xmax = 16.5592
fmax = 51.7476
deltaf = 0.252434 На данный момент при неравенстве левого и правого бина лучшую оценку дает расчет по двум бинам, причем умножать выборку на окно не нужно, это только ухудшает результат.
Сообщение отредактировал Pridnya - Nov 10 2015, 12:44
Эскизы прикрепленных изображений
|
|
|
|
|
|
|
|
Nov 10 2015, 12:41
|
Гуру
Группа: Свой
Сообщений: 3 106
Регистрация: 18-04-05
Пользователь №: 4 261
|
Цитата(Pridnya @ Nov 10 2015, 15:31) PS: И даже после умножения выборки на это правильное окно результат улучшился незначительно. Можно попробовать уменьшить "ширину" окна до минимально возможного: g[i]= expf(-((0.16F/ 0.02F)*(0.16F/ 0.02F)*(i-SAMPLES/2.0F+0.5F)*(i-SAMPLES/2.0F+0.5F)/((SAMPLES/2.0F)*(SAMPLES/2.0F))));
|
|
|
|
|
|
|
|
Nov 10 2015, 12:52
|
Частый гость
Группа: Свой
Сообщений: 142
Регистрация: 11-01-11
Из: Орел
Пользователь №: 62 159
|
Цитата(blackfin @ Nov 10 2015, 15:41) Можно попробовать уменьшить "ширину" окна до минимально возможного:
g[i]= expf(-((0.16F/0.02F)*(0.16F/0.02F)*(i-SAMPLES/2.0F+0.5F)*(i-SAMPLES/2.0F+0.5F)/((SAMPLES/2.0F)*(SAMPLES/2.0F)))); Окно стало совсем узким. Результат параболической аппроксимации еще чуть улучшился. Результат по двум бинам ухудшился. Цитата int main(int argc, char** argv) {
float fin, x1, x2, x3, y1, y2, y3, xmax, fmax, deltaf;
// *************************************************************************** //
// Параболическая аппроксимация по вашей формуле.
fin = 52.00; // Гц
// Бины без окна.
x1 = 16; y1 = 456.366F;
x2 = 17; y2 = 793.953F;
x3 = 18; y3 = 205.923F;
xmax = x2 + 0.25*(y3 - y1)/(2*y2 - y3 - y1);
fmax = xmax*3200.0/1024.0;
deltaf = fin - fmax;
cout <<"fin = " << fin << endl;
cout <<"xmax = " << xmax << endl;
cout <<"fmax = " << fmax << endl;
cout <<"deltaf = " << deltaf << endl;
/*
fin = 52
xmax = 16.9324
fmax = 52.9136
deltaf = -0.91362
*/
// Бины после окна "весовое окно узкое".
x1=16; y1=109.043F;
x2=17; y2=110.226F;
x3=18; y3=103.154F;
xmax = x2 + 0.25*(y3 - y1)/(2*y2 - y3 - y1);
fmax = xmax*3200.0/1024.0;
deltaf = fin - fmax;
cout <<"fin = " << fin << endl;
cout <<"xmax = " << xmax << endl;
cout <<"fmax = " << fmax << endl;
cout <<"deltaf = " << deltaf << endl;
/*fin = 52
xmax = 16.8217
fmax = 52.5677
deltaf = -0.567665
*/
// *************************************************************************** //
// Способ по двум бинам (хорошо считает, если левый бин не равен правому))
xmax = (16*y1 + 17*y2)/(y1+y2);
cout <<"fmax = (16*y1 + 17*y2)/(y1+y2)" << endl;
fmax = xmax*3200.0/1024.0;
deltaf = fin - fmax;
cout <<"fin = " << fin << endl;
cout <<"xmax = " << xmax << endl;
cout <<"fmax = " << fmax << endl;
cout <<"deltaf = " << deltaf << endl;
/*fmax = (16*y1 + 17*y2)/(y1+y2)
fin = 52
xmax = 16.5027
fmax = 51.5709
deltaf = 0.429073
*/
return 0;
} Напрвавление движения ясно. Два бина без окна, только нужно разобраться возле точек N*3.125. Цитата // *************************************************************************** //
Fin - заданная частота
Fcalc - расчетная частота
Err - разница
Yleft - левый бин (относительно максимума)
Yleft - правый бин (относительно максимума)
FFT 1024, без окна
// *************************************************************************** //
Fin=45 Fcalc=44.9805 Err=0.0194887 Yleft=224.695 Ymax=764.906 Yright=496.823
В этом промежутке не более Err 0,03 Гц
Fin=46.7697 Fcalc=46.7731 Err=-0.00334748 Yleft=33.6973 Ymax=999.206 Yright=33.5755
Fin=46.7797 Fcalc=46.9675 Err=-0.187824 Yleft=30.3978 Ymax=999.448 Yright=30.4946
Fin=46.7897 Fcalc=46.9583 Err=-0.168642 Yleft=27.1192 Ymax=999.654 Yright=27.3905
Fin=46.7997 Fcalc=46.949 Err=-0.149342 Yleft=23.8616 Ymax=999.826 Yright=24.2634
Fin=46.8097 Fcalc=46.9396 Err=-0.129922 Yleft=20.6249 Ymax=999.962 Yright=21.1137
Fin=46.8197 Fcalc=46.9301 Err=-0.110381 Yleft=17.4092 Ymax=1000.06 Yright=17.9415
Fin=46.8297 Fcalc=46.9204 Err=-0.0907161 Yleft=14.2145 Ymax=1000.13 Yright=14.747
Fin=46.8397 Fcalc=46.9106 Err=-0.0709257 Yleft=11.0407 Ymax=1000.16 Yright=11.5307
Fin=46.8497 Fcalc=46.9007 Err=-0.0510076 Yleft=7.88788 Ymax=1000.16 Yright=8.29256
Fin=46.8597 Fcalc=46.8906 Err=-0.0309596 Yleft=4.75595 Ymax=1000.12 Yright=5.03299
В этом промежутке не более Err 0,03 Гц
Fin=49.8992 Fcalc=49.9022 Err=-0.00300645 Yleft=32.2884 Ymax=999.253 Yright=32.1374
Fin=49.9092 Fcalc=50.0883 Err=-0.179093 Yleft=28.9915 Ymax=999.485 Yright=29.052
Fin=49.9192 Fcalc=50.079 Err=-0.159874 Yleft=25.7155 Ymax=999.682 Yright=25.9436
Fin=49.9292 Fcalc=50.0697 Err=-0.140537 Yleft=22.4604 Ymax=999.844 Yright=22.8125
Fin=49.9392 Fcalc=50.0602 Err=-0.12108 Yleft=19.2262 Ymax=999.971 Yright=19.6588
Fin=49.9492 Fcalc=50.0507 Err=-0.101501 Yleft=16.013 Ymax=1000.06 Yright=16.4828
Fin=49.9592 Fcalc=50.041 Err=-0.0817981 Yleft=12.8206 Ymax=1000.12 Yright=13.2847
Fin=49.9692 Fcalc=50.0311 Err=-0.0619689 Yleft=9.64908 Ymax=1000.14 Yright=10.0648
Fin=49.9792 Fcalc=50.0212 Err=-0.0420115 Yleft=6.4985 Ymax=1000.13 Yright=6.82334
В этом промежутке не более Err 0,03 Гц
Fin=53.0287 Fcalc=53.0314 Err=-0.00270608 Yleft=30.8695 Ymax=999.305 Yright=30.7067
Fin=53.0387 Fcalc=53.209 Err=-0.17037 Yleft=27.5761 Ymax=999.527 Yright=27.6163
Fin=53.0486 Fcalc=53.1998 Err=-0.151113 Yleft=24.3035 Ymax=999.714 Yright=24.5031
Fin=53.0586 Fcalc=53.1904 Err=-0.131737 Yleft=21.0517 Ymax=999.867 Yright=21.3672
Fin=53.0686 Fcalc=53.1809 Err=-0.112241 Yleft=17.8207 Ymax=999.984 Yright=18.209
Fin=53.0786 Fcalc=53.1713 Err=-0.092622 Yleft=14.6106 Ymax=1000.07 Yright=15.0285
Fin=53.0886 Fcalc=53.1615 Err=-0.0728786 Yleft=11.4214 Ymax=1000.11 Yright=11.8262
Fin=53.0986 Fcalc=53.1516 Err=-0.0530086 Yleft=8.25301 Ymax=1000.13 Yright=8.60214
Fin=53.1086 Fcalc=53.1416 Err=-0.0330097 Yleft=5.10552 Ymax=1000.11 Yright=5.35662
Далее до Fin 55 Гц Err не более Err 0,03Гц.
// *************************************************************************** //
Сообщение отредактировал Pridnya - Nov 10 2015, 12:58
Эскизы прикрепленных изображений
|
|
|
|
|
|
|
|
3 чел. читают эту тему (гостей: 3, скрытых пользователей: 0)
Пользователей: 0
|
|
|
