Нравится? Делимся информацией!

понедельник, 14 января 2013 г.

Вычисление модуля вектора ( sqrt( I^2 + Q^2 ) ). Результаты




В статье производится тестирование алгоритма на производительность как на компьютере ,так и на целевой платформе (TMS320F28335).

Статья 3. Вычисление модуля вектора (работа над ошибками)



Часть 5.

Результаты профилирования на компьютере

Для получения результатов были написаны тестовые программы в Microsoft Visual C++ 2010 Express.
Для исследования времени выполнения алгоритма, за эталон была взята функция sqrt() из <cmath>. Счет времени производился при помощи rdtsc инструкции - для этих целей был создан класс CProfiler для профилирования кода.

  1. /* sqrt_profile.cpp */
  2. //объявление профилируемого алгоритма
  3. #include "sqrtIQ.hpp"
  4. #include "profiler.hpp"
  6. #include <iostream>
  7. #include <cmath>
  8. #include <fstream>
  10. int main()
  11. {
  12. std::fstream file;
  13. file.open("sqrt_profile_log.txt", std::ios::out | std::ios::trunc);
  15. float i = 0.0125f;
  16. float q = 0.567f;
  17. float temp = 0.0f;
  18. float resultLib = 0.f, resultMy = 0.f;
  19. DWORD clocksLib = 0, clocksMy = 0;
  21. CProfiler profile;
  23. profile.prefix();
  24. //профилируемый код
  25. temp = i*i + q*q;
  26. resultLib = std::sqrt( temp );
  27. // конец профилируемого кода
  28. clocksLib = profile.postfix();
  30. profile.clear();
  32. profile.prefix();
  33. resultMy = sqrtIQ( i, q );
  34. clocksMy = profile.postfix();
  36. file << "<cmath> : Result of sqrt( ( "<< i <<" )^2 + ( "<< q <<" )^2 ) = ";
  37. file << resultLib << "\t Elapsed Time, in clocks: " << clocksLib << std::endl;
  38. file << "sqrtIQ: Result of sqrtIQ( "<< i <<", "<< q <<" ) = " << resultMy ;
  39. file << "\t\t\t Elapsed Time, in clocks: " << clocksMy ;
  40. file << "\t Clock Ratio, stdLib/myFunc : " << (float) clocksLib / clocksMy <<std::endl;
  42. std::cout << " Profiling is complete. Refer to log file in project folder" << std::endl;
  44. file.close();
  45. std::cin.get();
  46. return 0;
  47. }

Результат в файле sqrt_profile_log.txt :

<cmath> : Result of sqrt( ( 0.0125 )^2 + ( 0.567 )^2 ) = 0.567138 Elapsed Time, in clocks: 5820
sqrtIQ: Result of sqrtIQ( 0.0125, 0.567 ) = 0.567138 Elapsed Time, in clocks: 1704 Clock Ratio, stdLib/myFunc : 3.41549


Интересно было сравнить абсолютную и относительную ( в % ) точность измерения, а так же количество заходов алгоритма в цикл while() (sqrtIQ.hpp, 43 строка). Все измерения производятся в файле sqrt_iterations.cpp:

  1. /* sqrt_iterations.cpp */
  2. #include <iostream>
  3. #include <cmath>
  4. #include <fstream>
  6. #define ITERATIONS
  7. #include "sqrtIQ.hpp"
  10. int main() {
  12. std::fstream file;
  13. //ios::out - открыть для записи с начала файл
  14. //ios::trunc - очистить файл если он существует
  15. file.open( "sqrt_iterations_log.txt", std::ios::out | std::ios::trunc );
  17. float i, q, temp;
  18. //variables for sqrt's result
  19. float sqrtLib, sqrtMy;
  20. float error, errorPer;
  22. file << " " << std::endl;
  23. file << "/**********************************************************************************************/" << std::endl;
  24. file << "Опыт №1: Определяем количество итераций вычисления значения корня квадратного по методу Ньютона." << std::endl;
  25. file << "Это количество прохода цикла while( firstEstimation - nextEstimation != 0 ) { <..> } " << std::endl;
  26. file << "error = std::fabs( sqrtLib - sqrtMy);" << std::endl;
  27. file << "error, % вычисляется как: 100 * error/ sqrtLib; " << std::endl;
  28. file << " " << std::endl;
  30. for( int k = 0, g = 100; k < 100; ++k, --g) {
  32. i = ( float )k/10;
  33. q = ( float )g/10;
  34. temp = i*i + q*q;
  36. sqrtLib = sqrt( temp );
  37. sqrtMy = sqrtIQ(i, q);
  38. error = std::fabs( sqrtLib - sqrtMy);
  39. errorPer = 100 * error/ sqrtLib;
  40. file << "<cmath> sqrt("<< temp <<"): " << sqrtLib << "\t my sqrt("<< i <<","<< q <<"): " << sqrtMy;
  41. file << "\t perform in " << iterationCnt << " iterations" ;
  42. file << "\t error: " << error << "\t error, %: " << errorPer <<std::endl;
  44. } // for
  46. file.close();
  47. std::cin.get();
  48. return 0;
  49. }

Результаты в файле sqrt_iterations_log.txt (приведу неполный файл):
/**********************************************************************************************/
Опыт №1: Определяем количество итераций вычисления значения корня квадратного по методу Ньютона.
Это количество прохода цикла while( firstEstimation - nextEstimation != 0 ) { <..> }
error = std::fabs( sqrtLib - sqrtMy);
error, % вычисляется как: 100 * error/ sqrtLib;

<cmath> sqrt(100): 10 my sqrt(0,10): 10 perform in 1 iterations error: 0 error, %: 0 MIN
<cmath> sqrt(98.02): 9.9005 my sqrt(0.1,9.9): 9.90051 perform in 26 iterations error: 6.67572e-006 error, %: 6.74281e-005
<cmath> sqrt(96.08): 9.80204 my sqrt(0.2,9.8): 9.80205 perform in 26 iterations error: 1.14441e-005 error, %: 0.000116752
<cmath> sqrt(94.18): 9.70464 my sqrt(0.3,9.7): 9.70465 perform in 26 iterations error: 1.71661e-005 error, %: 0.000176886
<cmath> sqrt(92.32): 9.60833 my sqrt(0.4,9.6): 9.60833 perform in 34 iterations error: 1.90735e-006 error, %: 1.9851e-005
<cmath> sqrt(90.5): 9.51315 my sqrt(0.5,9.5): 9.51318 perform in 26 iterations error: 2.86102e-005 error, %: 0.000300744
<cmath> sqrt(88.72): 9.41913 my sqrt(0.6,9.4): 9.41916 perform in 26 iterations error: 3.33786e-005 error, %: 0.00035437
<cmath> sqrt(86.98): 9.32631 my sqrt(0.7,9.3): 9.32635 perform in 26 iterations error: 3.8147e-005 error, %: 0.000409025
<cmath> sqrt(85.28): 9.23472 my sqrt(0.8,9.2): 9.23476 perform in 26 iterations error: 4.29153e-005 error, %: 0.000464717
<cmath> sqrt(83.62): 9.1444 my sqrt(0.9,9.1): 9.14445 perform in 26 iterations error: 4.76837e-005 error, %: 0.000521453
<cmath> sqrt(82): 9.05539 my sqrt(1,9): 9.05544 perform in 26 iterations error: 5.14984e-005 error, %: 0.000568705
<cmath> sqrt(80.42): 8.96772 my sqrt(1.1,8.9): 8.96778 perform in 26 iterations error: 5.62668e-005 error, %: 0.000627437
<cmath> sqrt(78.88): 8.88144 my sqrt(1.2,8.8): 8.8815 perform in 26 iterations error: 6.00815e-005 error, %: 0.000676483
<cmath> sqrt(77.38): 8.79659 my sqrt(1.3,8.7): 8.79665 perform in 26 iterations error: 6.38962e-005 error, %: 0.000726374
<cmath> sqrt(75.92): 8.71321 my sqrt(1.4,8.6): 8.71328 perform in 26 iterations error: 6.67572e-005 error, %: 0.000766161 MAX
<cmath> sqrt(74.5): 8.63134 my sqrt(1.5,8.5): 8.63134 perform in 40 iterations error: 0 error, %: 0 MIN
<cmath> sqrt(73.12): 8.55102 my sqrt(1.6,8.4): 8.55102 perform in 39 iterations error: 1.90735e-006 error, %: 2.23055e-005
<cmath> sqrt(71.78): 8.47231 my sqrt(1.7,8.3): 8.47231 perform in 39 iterations error: 9.53674e-007 error, %: 1.12564e-005
<cmath> sqrt(70.48): 8.39524 my sqrt(1.8,8.2): 8.39524 perform in 39 iterations error: 9.53674e-007 error, %: 1.13597e-005
<cmath> sqrt(69.22): 8.31986 my sqrt(1.9,8.1): 8.31985 perform in 39 iterations error: 1.90735e-006 error, %: 2.29253e-005
<cmath> sqrt(68): 8.24621 my sqrt(2,8): 8.24621 perform in 39 iterations error: 9.53674e-007 error, %: 1.1565e-005
<cmath> sqrt(66.82): 8.17435 my sqrt(2.1,7.9): 8.17435 perform in 39 iterations error: 9.53674e-007 error, %: 1.16667e-005
<cmath> sqrt(65.68): 8.10432 my sqrt(2.2,7.8): 8.10432 perform in 41 iterations error: 9.53674e-007 error, %: 1.17675e-005
<cmath> sqrt(64.58): 8.03617 my sqrt(2.3,7.7): 8.03617 perform in 39 iterations error: 9.53674e-007 error, %: 1.18673e-005
<cmath> sqrt(63.52): 7.96994 my sqrt(2.4,7.6): 7.96994 perform in 39 iterations error: 1.43051e-006 error, %: 1.79488e-005
<cmath> sqrt(62.5): 7.90569 my sqrt(2.5,7.5): 7.90569 perform in 39 iterations error: 9.53674e-007 error, %: 1.20631e-005
<cmath> sqrt(61.52): 7.84347 my sqrt(2.6,7.4): 7.84347 perform in 39 iterations error: 9.53674e-007 error, %: 1.21588e-005
<cmath> sqrt(60.58): 7.78332 my sqrt(2.7,7.3): 7.78331 perform in 39 iterations error: 9.53674e-007 error, %: 1.22528e-005
<cmath> sqrt(59.68): 7.72528 my sqrt(2.8,7.2): 7.72528 perform in 39 iterations error: 1.43051e-006 error, %: 1.85173e-005
<cmath> sqrt(58.82): 7.66942 my sqrt(2.9,7.1): 7.66942 perform in 39 iterations error: 1.43051e-006 error, %: 1.86521e-005
<cmath> sqrt(58): 7.61577 my sqrt(3,7): 7.61577 perform in 39 iterations error: 9.53674e-007 error, %: 1.25224e-005
<cmath> sqrt(57.22): 7.56439 my sqrt(3.1,6.9): 7.56439 perform in 39 iterations error: 9.53674e-007 error, %: 1.26074e-005
<cmath> sqrt(56.48): 7.51532 my sqrt(3.2,6.8): 7.51532 perform in 39 iterations error: 9.53674e-007 error, %: 1.26897e-005
<cmath> sqrt(55.78): 7.4686 my sqrt(3.3,6.7): 7.4686 perform in 39 iterations error: 1.43051e-006 error, %: 1.91537e-005
<cmath> sqrt(55.12): 7.42428 my sqrt(3.4,6.6): 7.42428 perform in 42 iterations error: 4.76837e-007 error, %: 6.42267e-006
<cmath> sqrt(54.5): 7.38241 my sqrt(3.5,6.5): 7.38241 perform in 42 iterations error: 4.76837e-007 error, %: 6.4591e-006
<cmath> sqrt(53.92): 7.34302 my sqrt(3.6,6.4): 7.34302 perform in 39 iterations error: 9.53674e-007 error, %: 1.29875e-005
<cmath> sqrt(53.38): 7.30616 my sqrt(3.7,6.3): 7.30616 perform in 39 iterations error: 9.53674e-007 error, %: 1.3053e-005
<cmath> sqrt(52.88): 7.27186 my sqrt(3.8,6.2): 7.27186 perform in 42 iterations error: 4.76837e-007 error, %: 6.55729e-006
<cmath> sqrt(52.42): 7.24017 my sqrt(3.9,6.1): 7.24016 perform in 39 iterations error: 9.53674e-007 error, %: 1.3172e-005
< не полный файл sqrt_iterations_log.txt  >

Графическая интерпретация полученных результатов по поводу неточности вычислений функции sqrtIQ() представлена ниже. Первый график (первые два рисунка) показывает зависимость ошибки от отношения синфазной и квадратурной составляющих цифрового сигнала. Синфазная составляющая - это вектор от 0 до 10 с шагом 0,1; квадратурная составляющая - это вектор от 10 до 0 с шагом 0,1. На втором графике равномерно по оси абсцисс просто выводятся результаты вычисленной ошибки.





Вывод: полученный алгоритм быстрее библиотечной функции, примерно, в 3 раза, за что мы расплачиваемся ухудшением точности, но это ухудшение незначительное для наших целей (максимум погрешности 0,00076 % ). Если же сравнивать с аппроксимацией 4 порядка (видимо, по  ряду Тейлора) по этой ссылке , то точность sqrtIQ() в 10 раз лучше.

Примечание 1. Если в алгоритме, что записан в файле sqrtIQ.hpp в строках 42, 46 заменить операцию деления на 2 умножением на 0,5, то по скорости получится еще больший выигрыш (до замены соотношение было около 3.3):

<cmath> : Result of sqrt( ( 0.0125 )^2 + ( 0.567 )^2 ) = 0.567138 Elapsed Time, in clocks: 5792
sqrtIQ: Result of sqrtIQ( 0.0125, 0.567 ) = 0.567138 Elapsed Time, in clocks: 1176 Clock Ratio, stdLib/myFunc : 4.92517

<cmath> : Result of sqrt( ( 0.0125 )^2 + ( 0.567 )^2 ) = 0.567138 Elapsed Time, in clocks: 4915
sqrtIQ: Result of sqrtIQ( 0.0125, 0.567 ) = 0.567138 Elapsed Time, in clocks: 1172 Clock Ratio, stdLib/myFunc : 4.19369



Примечание 2. Графики строились в среде МатЛаб при помощи шаблонного скрипта, который требовалось лишь немного подправить. Сам скрипт:

  1. %This script is used for quick plotting the experimental data from file
  3. %specify the name of file with data that must be plotted
  4. file = 'sqrt_iterations_error.dat';
  6. fid = fopen( file, 'r');
  7. if fid == -1
  8.   error('file is not opened');
  9. end
  11. %read the data
  12. data = fscanf( fid,'%f', inf );
  14. %DELETE - example specified
  15. %
  16. %form additional data for plotting data
  17. % in-phase and quadrature samples that used for calculating sqrt()
  18. I = 0:0.1:10;
  19. Q = sort( I, 'descend');    %contains 10 9.9 9.8 ...
  21. %preallocate 1D column vector
  22. quotient = zeros( 1, 100, 'single');
  23. for i=1:100
  24.   quotient(i) = I(i) / Q(i);
  25. end
  26. display( quotient );
  27. %
  28. %DELETE
  30. %print into screen
  31. fprintf(1, '%.11fn', data );
  33. plot( quotient, data );
  36. %MODIFY - example specified
  37. %the first graph
  38. title('Error of sqrtIQ() function (with respect to <cmath> sqrt())');
  39. xlabel('Ratio of the in-phase to quadrature: I/Q');
  40. ylabel('Error, %');
  42. %the second graph
  43. %plot( data );
  44. %title('Error of sqrtIQ() function (with respect to <cmath> sqrt())');
  45. %xlabel('Uniform scale');
  46. %
  47. %MODIFY
  49. fclose ( fid );

Часть 6.
Результаты работы на целевой платформе TMS320F28335.
Проверка работы алгоритма производилась в программе, которая является частью цифрового детектора. В этой программе была реализована первая ступень фильтрации-децимации цифрового детектора как в синфазном, так и квадратурном канале, и как логическое завершение - находился модуль вектора по I-Q отсчетам при помощи шаблонной функции sqrtIQ():


  1. #include "fircoeff.h" //from matlab
  2. #include "filter.hpp"
  3. #include "adcmodel.hpp"
  4. #include "multiplication.hpp"
  5. #include "sqrtIQ.hpp"
  7. #define WDCR        *((volatile int *)0x7029)    /* WD Control reg */
  8. #define DISABLE_WD  0x0068
  10. void Disable_WD(void)
  11. {
  12.   asm(" eallow");
  13.   WDCR |= DISABLE_WD;
  14.   asm(" edis");
  15. }
  17. //this is the samples amount in the *.dat file
  18. const int kLength = 5208;
  19. float inPhaseSamples[ kLength ];
  20. float quadratureSamples[ kLength ];
  22. static void dataIO();
  24. float g_sampleI = 0, g_sampleQ = 0;
  25. float g_firOutI = 0, g_firOutQ = 0;
  26. float g_sqrtOut = 0;
  28. int main() {
  30. Disable_WD();
  32. //function dataIO() to catch the breakPoint and read
  33. //all massive of samples into inPhaseSamples buffer
  34. dataIO();
  35. //load samples into quadratureSamples
  36. dataIO();
  38. CAdcModel adcI( inPhaseSamples, kLength );
  39. CAdcModel adcQ( quadratureSamples, kLength);
  40. //filter for in-phase channel
  41. dsp::CFirFilter firFilterI( FIR_coeff, BL, 10 );
  42. //filter for quadrature channel
  43. dsp::CFirFilter firFilterQ( FIR_coeff, BL, 10 );
  45. bool isReadyFirI, isReadyFirQ;
  47. while(1) {
  48. //reading ADC
  49. adcI.getSample( g_sampleI );
  50. adcQ.getSample( g_sampleQ );
  52. isReadyFirI = firFilterI.filtrate( g_sampleI, g_firOutI );
  53. isReadyFirQ = firFilterQ.filtrate( g_sampleQ, g_firOutQ );
  55. if( ( isReadyFirI == true ) &&
  56. ( isReadyFirQ  == true ) ) {
  58. //here "g_out"/"g_firOutI/Q " contains proof FIR-decimator-output value.
  60. g_sqrtOut = dsp::sqrtIQ( g_firOutI, g_firOutQ );
  61. //Function to view the output of sqrt()
  62. dataIO();
  63. }
  65. }// while
  67. return 0;
  68. }
  70. /**
  71. @function dataIO
  72. The dataIO function acts as a placeholder.It is a convenient place to connect
  73. a Probe Point (or BreakPoint in the new CCS) that injects data from a PC file.
  74. @detailed - for details see spru301, 4.3 Adding a Probe Point for File I/O
  75. */
  76. static void dataIO() {
  77. /* do data I/O */
  78. return;
  79. }

Результаты работы алгоритма, оказались неожиданными:


Что то не так! Должна была быть на графике низкочастотная огибающая:


Первая идея - что-то не так с алгоритмом sqrtIQ(). Действительно. Оказалось, что отрицательные значения на выходе sqrtIQ(), появляются тогда, когда значения ВЧ-заполнения входного сигнала становятся отрицательными (когда ВЧ несущая меньше нуля). В этом случае алгоритм “альфа максимум + бета минимум” выдает отрицательное начальное приближение, и уже начиная с этого отрицательного значения, запускается расчет корня по алгоритму Ньютона, который, естественно, сходится к отрицательному числу.
Для преодоления проблемы написана шаблонная функция вычисления модуля (вычисляет по отдельности модуль для I и Q):

  1. template< class T >
  2. void absIQ( T& i, T& q ) {
  4. i = ( i >= (T)0 ) ? i : -i;
  5. q = ( q >= (T)0 ) ? q : -q;
  7. }

absIQ() была добавлена в начало функции sqrtIQ() (до начала "alpha max plus beta min"). Результат не заставил долго себя ждать:

Автор: на
Ярлыки: , , , , , ,

Комментариев нет:

Отправить комментарий

Подписаться на: Комментарии к сообщению (Atom)