Алгоритм Гёрцеля на fixed point

Если ваш вопрос не влез ни в одну из вышеперечисленных тем, вам сюда.
Ответить
Аватара пользователя
YS
Друг Кота
Сообщения: 7518
Зарегистрирован: Вс мар 29, 2009 22:09:05
Контактная информация:

Алгоритм Гёрцеля на fixed point

Сообщение YS »

Итак, товарищи, я задался идеей реализовать алгоритм Гёрцеля. С уклоном в эмбед, конечно. И потому с фиксированной точкой. На всякий случай - я в курсе, что есть готовые реализации, но мне интересно, почему не работает у меня.

Итак, для начала я реализовал этот алгоритм с плавающей точкой:

Код: Выделить всё

float Gz_GetFreqMagnitude(uint8_t block[BLOCK_SIZE],gz_tgt_freq_t *tgt_f)
{
    float Q0,Q1=0,Q2=0;
    uint8_t i;

    for (i=0; i<BLOCK_SIZE; i++)
    {
        Q0=(tgt_f->cf)*Q1-Q2+(block[i]/255.0);
        Q2=Q1;
        Q1=Q0;
    }

    return sqrt(Q1*Q1+Q2*Q2-Q1*Q2*(tgt_f->cf));
}
tgt_f - структура, содержащая служебную информацию и в том числе необходимый поворотный множитель.

Код: Выделить всё

typedef struct {

    float tgt_freq;
    float k;
    float w;
    float cos_w;
    float sin_w;
    float cf;

} gz_tgt_freq_t;
Она вычисляется другой функцией:

Код: Выделить всё

void Gz_ComputeFreqDescriptor(gz_tgt_freq_t *dsc)
{
    dsc->k=floor(0.5+((BLOCK_SIZE*(dsc->tgt_freq))/SAMPLING_RATE_HZ));

    dsc->w=((2*M_PI)/BLOCK_SIZE)*(dsc->k);

    dsc->cos_w=cos(dsc->w);

    dsc->sin_w=sin(dsc->w);

    dsc->cf=2*(dsc->cos_w);
}
Сначала заполняется поле tgt_freq, потом вызывается Gz_ComputeFreqDescriptor(), потом - Gz_GetFreqMagnitude() с вычисленными значениями в структуре.

Эта реализация с float работает. Чудеса начинаются при попытке перенести ее на фиксированную точку. Для этой благородной цели я написал файл, определяющий операции с фиксированной точкой:

Код: Выделить всё

typedef int32_t fixed;
#define FIX_BITS       10

#define fix_fill(wpart,fpart)       (((wpart) << FIX_BITS) | (((fpart) * (1 << FIX_BITS)) / 1000))

#define fix_add(num1,num2)          ((num1) + (num2))
#define fix_sub(num1,num2)          ((num1) - (num2))

#define fix_mul(num1,num2)          (((num1) * (num2)) >> FIX_BITS)
#define fix_div(num1,num2)          (((num1) << FIX_BITS) / (num2))

#define fix_whole(num)              (uint32_t)((num) >> FIX_BITS)
#define fix_frac(num)               (uint32_t)((((num) & ((1 << FIX_BITS)-1)) * 1000) / (1 << FIX_BITS))

#define fix_sqr(num)                (fix_mul((num),(num)))
И в отдельности он тоже работает!

Но, когда я пытаюсь совместить эти два файла, начинаются чудеса. Совмещаю я так:

Код: Выделить всё

uint32_t GzInt_GetFreqMagnitude(uint8_t block[BLOCK_SIZE],fixed *tgt_f)
{
    fixed Q0,Q1=0,Q2=0;
    uint8_t i;

    for (i=0; i<BLOCK_SIZE; i++)
    {
        //Q0=(tgt_f->cf)*Q1-Q2+(block[i]/255.0);
        Q0=fix_add(fix_sub(fix_mul((*tgt_f),Q1),Q2),fix_fill(block[i],0));
        Q2=Q1;
        Q1=Q0;
    }

    //return sqrt(Q1*Q1+Q2*Q2-Q1*Q2*(tgt_f->cf));
    return fix_whole(fix_sub(fix_add(fix_sqr(Q1),fix_sqr(Q2)),fix_mul(fix_mul(Q1,Q2),(*tgt_f))));
}
Коэффициент пересчитывается так:

Код: Выделить всё

fixed Gz_ConvertCoefficient(gz_tgt_freq_t *float_coeff)
{
    //printf("Freq: %.2f,coeff: %.5f\n",float_coeff->tgt_freq,float_coeff->cf);

    double f_part,i_part;
    fixed fixed_coeff;

    f_part=modf(float_coeff->cf,&i_part);

    fixed_coeff=fix_fill((int32_t)i_part,(int32_t)(f_part*1000));


    //printf("IF: %i;%i\n",(int32_t)i_part,(int32_t)(f_part*1000));
    //printf("IG: %i,%i\n",fix_whole(fixed_coeff),fix_frac(fixed_coeff));

    return fixed_coeff;
}
Виден отладочный вывод - я проверял, все переводится 100% правильно. Но сам алгоритм рассыпается. Версия с плавающей точкой при тестах дает четкий максимум на корректной частоте. С фиксированной - какой-то левый случайный набор значений.

Отчего так? Неужели алгоритму Гёрцеля нужна такая точность, что он рассыпается от ограниченной точности фиксированной точки?

Я уже всю голову сломал, ничего не помогает. И главное - по отдельности все работает!
Разница между теорией и практикой на практике гораздо больше, чем в теории.
Реклама
Аватара пользователя
igor-x
Мудрый кот
Сообщения: 1817
Зарегистрирован: Пн ноя 29, 2010 15:58:43

Re: Алгоритм Гёрцеля на fixed point

Сообщение igor-x »

если не секрет, на какой контроллер планируете ставить этот алгоритм?
я очень давно тоже пытался этот алгоритм реализовать на АВР, но у меня недостаток теории по этому вопросу . поэтому проект повис ...
Реклама
Аватара пользователя
Kavka
Мудрый кот
Сообщения: 1810
Зарегистрирован: Чт июн 10, 2010 08:55:35
Откуда: Сибирские Афины

Re: Алгоритм Гёрцеля на fixed point

Сообщение Kavka »

YS, вероятно, проблема именно в точности. Какое минимальное и максимальное число у тебя в фиксированной точке? Есть места где значение может уходить за эти пределы? Дробную часть увеличивать не пробовал? Хотя бы для эксперимента.
В качестве бредовой идеи: сведи варианты в одну прогу и сравни значения по шагам - станет ясно где проблема и в чём именно. Возможно придётся наставить проверок в места, чтобы избегать выхода значений за пределы допустимого диапазона.

Вот что попалось с ходу - http://www.avrfreaks.net/wiki/index.php ... :AVR_Float
float : 32 bit IEEE 754
double : identical to float (32 bit IEEE 754)
Type__________ Size____ Range (+/-)__________________ Exponent_ Mantissa
float__________ 32 bits +-1.18E-38 to +-3.39E+38___________ 8 bits 23 bits
Когда уже ничего не помогает - прочтите, наконец, инструкцию.
Лучший оптимизатор находится у вас между ушей. (Майкл Абраш, программист Quake и QuakeII)
Избыток информации ведёт к оскудению души - Леонтьев А. (сказано в 1965 г.)
Аватара пользователя
YS
Друг Кота
Сообщения: 7518
Зарегистрирован: Вс мар 29, 2009 22:09:05
Контактная информация:

Re: Алгоритм Гёрцеля на fixed point

Сообщение YS »

на какой контроллер планируете ставить этот алгоритм?
И на AVR в том числе. Хочу получить кросс-платформенный алгоритм с фиксированной точкой. Но главная цель - прокачать скилл вычислительной математики. :)
о у меня недостаток теории по этому вопросу
Я делал по мануалу отсюда. Очень толковое описание без углубления в дебри математики.

А вообще преобразование Гёрцеля - частный случай ДПФ и выводится из него.
Какое минимальное и максимальное число у тебя в фиксированной точке?
У меня формат 22.10, так что максимум/минимум уходят за два миллиона (с учетом знаковости). На дробную часть десять бит, соответственно погрешность представления что-то около 0.001.
Дробную часть увеличивать не пробовал?
Пробовал. Результат непредсказуемо меняется.
сведи варианты в одну прогу и сравни значения по шагам
Похоже, так и сделаю...

Еще думаю написать универсальную реализацию - задефайнить арифметические действия, и подставлять либо float, либо фиксированную точку.

UPD:

Заработало. Похоже, все-таки это было переполнение.

В расчете на каждой итерации заменил

Код: Выделить всё

Q0=fix_add(fix_sub(fix_mul((*tgt_f),Q1),Q2),fix_fill(block[i],0));
на

Код: Выделить всё

Q0=fix_add(fix_sub(fix_mul((*tgt_f),Q1),Q2),fix_div(fix_fill(block[i],0),fix_fill(255,0)));
Т.е., привожу значение сэмпла к диапазону [0;1]. При точности вещественной части 12 бит результат почти не отличается.

Причешу, протестирую, оформлю и, если время будет, запилю статью...
Разница между теорией и практикой на практике гораздо больше, чем в теории.
Реклама
Эиком - электронные компоненты и радиодетали
Ответить

Вернуться в «Разные вопросы по МК»