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

В данной статье объясняется как удалось улучшить работу моего алгоритма с 6000 тактов процессорного времени до 284 тактов!!!

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

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

В листинге sqrtIQ.hpp обнаружена неточность, стоившая не соответствующих реальности показаний. Ошибка обнаружена на этапе

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

Вся загвоздка оказалась в строке 49 файла sqrtIQ.hpp:

quotient = radicand/firstEstimation;

Данная строка должна идти до строки 46. В итоге, с учетом всех замечаний, подмеченных в прошлых статьях и в этой, листинг должен

быть (изменения выделены жирным и прокомментированы):

1.// Version 2.0

2.#ifndef SQRTIQ_HPP

3.#define SQRTIQ_HPP

4.

5.#ifdef ITERATIONS

6.//for debug

7.int iterationCnt;

8.#endif

9.

10.template< class T >

11.void absIQ( T& i, T& q ) {

12.

13. i = ( i >= (T)0 ) ? i : -i;

14. q = ( q >= (T)0 ) ? q : -q;

15.}

16.

17.template< class T >

18.T sqrtIQ( T i, T q ) {

19.

20.#ifdef ITERATIONS

21.//for debug

22.iterationCnt = 1;

23.#endif

24.

25.//first estimation for Newton's algorithm. This estimation is found

26.//by means of "alpha max plus beta min" algorithm

27.T firstEstimation;

28.//current value of estimation

29.T nextEstimation;

30.T radicand;

31.T quotient;

32.//coefficients for "alpha max plus beta min":

33.//alpha = 1.0, beta = 0.5 therefore they don't need to be present in memory

34.

35.absIQ( i, q );//без этой строки 1е приближение м/б отрицательным - это ошибка

36.

37.//begin - "alpha max plus beta min" algorithm

38.if( i > q ) {

39. firstEstimation = i +  q*0.5;//везде замена /2

40.}

41.else {

42. firstEstimation = q +  i*0.5;

43.}

44.//end - "alpha max plus beta min" algorithm

45.//begin - Newton's algorithm

46. i *= i;

47. q *= q;

48. radicand = i + q;

49.

50. quotient = radicand/firstEstimation;

51. nextEstimation = ( firstEstimation + quotient )*0.5;

52. while( firstEstimation - nextEstimation != 0 ) {

53. firstEstimation = nextEstimation;

54.

55. quotient = radicand/firstEstimation;//поднял эту строку на две позиции вверх

56. nextEstimation = ( firstEstimation + quotient )*0.5;

57.

58. #ifdef ITERATIONS

59. //for debug

60. ++iterationCnt;

61. #endif

62.} //while

63.

64.return nextEstimation; //без разницы что возвращатть, т.к. они равны

65.

66.//end - Newton's algorithm

67.}

68.

69.#endif // SQRTIQ_HPP

70.

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

1й прогон:

<cmath> :___Result of sqrt( ( 0.0125 )^2 + ( 0.567 )^2 ) = 0.567138___Elapsed Time, in clocks: 5286
sqrtIQ:______Result of sqrtIQ( 0.0125, 0.567 ) = 0.567138__________Elapsed Time, in clocks: 900___Clock Ratio, stdLib/myFunc : 5.87333

2й прогон:

<cmath> :___Result of sqrt( ( 0.0125 )^2 + ( 0.567 )^2 ) = 0.567138___Elapsed Time, in clocks: 9515
sqrtIQ:______Result of sqrtIQ( 0.0125, 0.567 ) = 0.567138__________Elapsed Time, in clocks: 1916__Clock Ratio, stdLib/myFunc : 4.96607
и т.д..



Для выполнения на целевой платформе вместо цикла While() {} оставляем 3-4 итерации, т.е. из экспериментов этого достаточно для получения точного результата. Программа становится следующей:

1. /* Version 3.0 */
2. /**
3. @file sqrtIQ.hpp
4. @brief target platform - TMS320F28335
5. @date 25.12.2012
6. @note In this file template sqrtIQ(T i, T q) function is implemented that perform
7. calculation over i (in-phase) and q (quadrature) samples realizing
8. expression: sqrt( i^2 + q^2 )
9. */
11. #ifndef SQRTIQ_HPP
12. #define SQRTIQ_HPP
14. #include "abs.hpp"
16. namespace dsp {
18. /**
19. @function sqrtIQ() template function consists of two algorithm. The main
20. is Newton algorithm for computing the sqrt. For this algorithm
21. the first estimation must be known.We obtain the first estimation
22. by means of "alpha Max plus beta Min" algorithm ( sqrt
23.                             linearization algorithm ).
24. @param[in] i data sample  from in-phase channel of digital detector
25. @param[in] q data sample  from quadrature channel of digital detector
26. @return sqrt( i^2 + q^2 )
27. */
29. template< class T >
30. T sqrtIQ( T i, T q ) {
32. //first estimation for Newton's algorithm. This estimation is found
33. //by means of "alpha max plus beta min" algorithm
34. T firstEstimation;
35. //current value of estimation
36. T nextEstimation;
37. T radicand;
38. T quotient;
39. //coefficients for "alpha max plus beta min":
40. //alpha = 1.0, beta = 0.5 therefore they don't need to be present in memory
42. //without this func some samples equals negative values - it's
43. //impossible
44. dsp::absIQ( i, q );
46. //begin - "alpha max plus beta min" algorithm
47. if( i > q ) {
48. firstEstimation = i +  q*0.5;
49. }
50. else {
51. firstEstimation = q +  i*0.5;
52. }
53. //end - "alpha max plus beta min" algorithm
54. //begin - Newton's algorithm
55. i *= i;
56. q *= q;
57. radicand = i + q;
59. quotient = radicand/firstEstimation;
60. nextEstimation = ( firstEstimation + quotient )*0.5;
62. firstEstimation = nextEstimation;
63. quotient = radicand/firstEstimation;
64. nextEstimation = ( firstEstimation + quotient )*0.5;
66. firstEstimation = nextEstimation;
67. quotient = radicand/firstEstimation;
68. nextEstimation = ( firstEstimation + quotient )*0.5;
70. firstEstimation = nextEstimation;
71. quotient = radicand/firstEstimation;
72. nextEstimation = ( firstEstimation + quotient )*0.5;
74. return nextEstimation;
76. //end - Newton's algorithm
77. }
79. } //namespace dsp::
81. #endif // SQRTIQ_HPP

Результаты профилирования на целевой платформе TMS320F28335:
До последних внесенных изменений на выполнение функции уходило около 6000 тактов процессора. Сейчас уже лучше - 1112 тактов. С учетом того, что в алгоритме используется деление (что нехорошо) - 4 деления, каждое из которых занимает 227 тактов (см. здесь ), полученное профилировщиком время кажется правдоподобным.

Последний эксперимент.

Чтобы выжать максимум на TMS320, воспользуемся его библиотекой Fast RTS library, в которой представлены быстрые реализации многих стандартных функций, в том числе и деления. Посмотрим что это нам даст!
Как выжать максимум - читаем в статье увеличиваем производительность стандартных функций .
Настроили проект. Запускаем на отладку.

Полученные данные производительности в 5 раз лучше!! 

Получили результат 284 такта!

Вывод: получили неплохую производительность в 284 такта. Но есть более быстрые алгоритмы, которые уже реализованы в Fast RTS library (извлечение корня, инверсный корень квадратный, деление), поэтому проще, эффективней и экономичней использовать следующий вариант:

g_sqrtOut = std::sqrt( I*I + Q*Q );

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