Аппроксимация квадратного корня из суммы двух квадратов на микроконтроллере

Я работаю над реализацией алгоритма FFT в сборке на 8-битном микроконтроллере (HCS08) для развлечения. Как только алгоритм будет завершен, у меня будет массив 8-битных пар вещественных и мнимых чисел, и я хочу найти величину каждого из этих значений. То есть, если x сложный, я хочу найти

|x| = sqrt(Re{x}^2 + Im{x}^2)

Теперь мне доступны 16-битный регистр и 8-битный регистр. Я думал просто возвести их в квадрат, сложить и извлечь квадратный корень из результата, но это создает проблему: максимально возможное значение суммы квадратов двух 8-битных чисел составляет ~ 130k, что больше, чем максимальное значение, которое может содержать 16-битный регистр (65,5 КБ).

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

Подводя итог: скажем, у меня есть вектор с двумя 8-битными компонентами, и я хочу найти длину вектора. Как я могу аппроксимировать это без фактического вычисления квадратов и квадратных корней?

Спасибо!


microcontroller signal-processing assembly numerical-methods
person user599599    schedule 03.04.2011    source источник
comment
Алгоритм CORDIC (en.wikipedia.org/wiki/CORDIC) можно использовать для ротации вектор <x,y> в какой-то новый вектор <x1,0> (или эквивалентно <0,y1>. x1 (или y1) дает величину исходного вектора, а CORDIC можно реализовать без умножения. Я никогда не делал этого сам, и понятия не имею, насколько это сложно .   -  person mtrw    schedule 03.04.2011
comment
Это случайно не для аудио? Собираетесь ли вы впоследствии вычислять log10, чтобы получить значение в дБ?   -  person Paul R    schedule 03.04.2011
comment
Зависит от цели: если вам нужна длина, другого способа вычисления нет, но когда вам действительно нужна норма (обычно это длина), вы можете использовать другую норму вместо нормы L2 по умолчанию, например манхэттенское расстояние (= |реальное|+|изображение|).   -  person flolo    schedule 03.04.2011
comment
@Paul R: Да, это для аудиопроекта, над которым я работаю. Аппаратное обеспечение, с которым я взаимодействую, ожидает линейное напряжение и преобразует его в логарифмическую шкалу.   -  person user599599    schedule 04.04.2011
comment
@ user599599: Хорошо, в таком случае вы, вероятно, можете избавиться от sqrt - см. Ответ ниже.   -  person Paul R    schedule 04.04.2011

Ответы (6)


arrow_upward
3
arrow_downward

Если сумма больше 65535, разделите ее на 4 (сдвиньте вправо на 2 бита), извлеките квадратный корень и умножьте на 2. Вы потеряете один бит точности, и, естественно, не гарантируется, что результат уместится в 8 биты.

person Mark Ransom    schedule 03.04.2011
comment
Спасибо за ответ. Меня беспокоит только то, что если сумма больше 65535, она переполнится, и я не смогу узнать об этом. (У меня есть только 16-разрядный регистр, поэтому добавление двух 16-разрядных чисел может привести к непредсказуемым результатам.) Я думаю, что могу добиться того же, сначала разделив Re{x} и Im{x} на 2, а затем умножив конечный результат. ответ на 2; это звучит эквивалентно тому, что вы предложили? - person user599599; 03.04.2011
comment
Вы уже приняли этот ответ, так что я думаю, вы поняли: разделите входные данные на 4 и умножьте выходные данные на 2. - person Mark Ransom; 04.04.2011

arrow_upward
6
arrow_downward

Существует веб-страница с описанием быстрого оценщика величины. Основная идея состоит в том, чтобы получить метод наименьших квадратов (или другое высокое качество), подходящее к уравнению:

Mag ~= Alpha * max(|I|, |Q|) + Beta * min(|I|, |Q|)

для коэффициентов Альфа и Бета. Несколько пар коэффициентов перечислены со среднеквадратичными ошибками, максимальными ошибками и т. д., включая коэффициенты, подходящие для целочисленных АЛУ.

person hotpaw2    schedule 03.04.2011
comment
Похоже, для этого приложения подойдет один из вариантов 61/64. - person Mark Ransom; 04.04.2011

arrow_upward
0
arrow_downward

Ну, вы можете написать x в полярной форме:

x = r[cos(w) + i sin(w)]

где w = arctan(Im(x)/Re(x)), значит

|x| = r = Re(x)/cos(w)

Здесь нет больших чисел, но, возможно, вы потеряете точность тригонометрических функций (то есть, если у вас есть доступ к тригонометрическим функциям :-/)

person Eelvex    schedule 03.04.2011
comment
Хм, интересная идея. К сожалению, у меня нет доступа к триггерным функциям, и в любом случае микроконтроллер не поддерживает операции с плавающей запятой, поэтому я в значительной степени ограничен базовыми целочисленными операциями. Однако я планирую иметь таблицу поиска триггеров, поэтому я буду иметь это в виду. - person user599599; 03.04.2011

arrow_upward
0
arrow_downward

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

|x| ~ max(|Re{x}|,|Im{x}|) + min(|Re{x}|,|Im{x})/2;

Это приведет к переоценке |x| где-то между 0 и 12%.

person swestrup    schedule 03.04.2011

arrow_upward
0
arrow_downward

Если вы впоследствии собираетесь преобразовывать величину в дБ, то вы вообще отказываетесь от операции sqrt. т.е. если ваш расчет:

magnitude = sqrt(re*re+im*im); // calculate magnitude of complex FFT output value
magnitude_dB = 20*log10(magnitude); // convert magnitude to dB

вы можете переписать это как:

magnitude_sq = re*re+im*im; // calculate squared magnitude of complex FFT output value
magnitude_dB = 10*log10(magnitude_sq);  // convert squared magnitude to dB
person Paul R    schedule 04.04.2011
comment
Хороший вопрос, но моя проблема в том, что log10 также требует больших вычислительных ресурсов. У меня все еще была бы проблема поиска ближайшего целого числа или использования таблицы поиска. - person user599599; 10.04.2011
comment
@ user599599: да, у вас все еще есть log, но раньше у вас было sqrt + log, а теперь у вас просто log. - person Paul R; 10.04.2011

arrow_upward
0
arrow_downward

Возможно, вы ограничены всего двумя регистрами, но вы можете посмотреть этот код по адресу http://www.realitypixels.com/turk/opensource/index.html Квадратный корень с фиксированной точкой Тригонометрические функции фиксированной точки с использованием CORDIC

person Reality Pixels    schedule 08.05.2016