Быстрое преобразование Фурье на Си

[linuxdude]

Мурик писал(а):

Сможете легко перенести программу скажем с PIC16 на STM32?

Ну я могу представить себе что вон тот Гёрцель писаный под AVR с гитхаба и на STM32 взлетит без особых изменений. Единственное что там придется поменять - как/откуда он сэмплы берет. Однако ж сам Гёрцель останется тем же самым. Правда в конкретно с тем кодом так делать не рекомендуется, автор очень в духе гитхаба лицензию указать вообще забыл и формально любой кто так сделает - как бы пират.

[YS]

О-о-о! Смотрите, коллеги, что я нашел!

Реализация целочисленного БПФ на процессорах с архитектурой ARM

[jcxz]

YS писал(а):

О-о-о! Смотрите, коллеги, что я нашел!
Реализация целочисленного БПФ на процессорах с архитектурой ARM

Что "О-о-о"? Что там такого интересного? К тому же - для устаревшего ядра...

[YS]

Ну, человек просил пример кода под ARM с целочисленной арифметикой - вот, ровно оно. Код будет выполняться и на M3, они совместимы.

А так по FFT много хороших публикаций есть. Вот, в частности.

[jcxz]

YS писал(а):

Ну, человек просил пример кода под ARM с целочисленной арифметикой - вот, ровно оно. Код будет выполняться и на M3, они совместимы.

Пример кода чего он просил? Видимо БПФ. А там что? Где там БПФ?
К тому же он просил для Cortex-M0, а не для ARM7.

YS писал(а):

А так по FFT много хороших публикаций есть.

Вот и почитайте хотя-бы сами, чтобы знать что такое БПФ.
И гуглом думаю все умеют пользоваться.....

[a5021]

Истеричка заявилась. Как обычно, ни одной полезной мысли, только вопли и визги.

[YS]

А там что?

А там - пример реализации ядра БПФ, элемента вычисления под названием "бабочка". Это, собственно, самый интересный этап БПФ.

[jcxz]

YS писал(а):

А там - пример реализации ядра БПФ, элемента вычисления под названием "бабочка". Это, собственно, самый интересный этап БПФ.

Ну так это - малая часть. Да и не особо оптимально написанная. А всё остальное, как пишет автор - на си.
Если хочется ассемблера, можно взять любой пример на си, скомпилить его для нужного ядра с максимальной оптимизацией, взять листинг и дальше допиливать оптимизировать его вручную.

[a5021]

а в конце концов -- в дурку, лечить таки шизофрению.

[YS]

Ну так это - малая часть.

Хехе, вот теперь я посоветую вам почитать хотя бы тот труд, ссылку на который я приводил. Производительность БПФ в основном определяется как раз реализацией "бабочки". Это основная часть.

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

Ну удачи, удачи...

[Repytw]

И вот интересно - среди кидания ссылок, требований подтянуть матчасть, обвинений друг друга в некомпетентности и других забавных проявлений человеческих недостатков кто-нибудь осмелится представить кусок кода на Си, т.е. подпрограммку, которая будет рассчитывать ДПФ целочисленного массива 12-битных данных, длинной хотя бы в 128 отсчётов.
Это не претензия.
Просто много говорится о ДПФ, обсуждаются возможности оптимизации, целочисленной реализации, бит-реверсной адресации какой-то абстрактной химеры, которой по-факту еще и нет.
Вот я, много ненужных сообщений назад, представил код, в котором есть такая подпрограммка - DFT.
Может попробует её совершенствовать, а не болтать без дела. Авось кому-нибудь и поможем.
Код компьютерный, писался для научных целей, но возможностей для оптимизации - вагон. Да еще и основная идея явно проглядывается.
Заодно посмотрим, кто в языке Си и ДПФ разбирается.
Ваши предложения, коллеги!

Спойлер

int DFT(FILE *SOURSE, FILE *SPECTR)
{
long k;
long n;
short temp;
double WN;
double wk;
double c;
double s;
double *XR;
double *XI;

COMPLEX *x;

x = (COMPLEX*)malloc(BUFFER*sizeof(COMPLEX));
XR = (double*)malloc(BUFFER*sizeof(double));
XI = (double*)malloc(BUFFER*sizeof(double));


rewind(SOURSE);
for(k=0; k<BUFFER;k++)
{
fread(&temp, 2 , 1 ,SOURSE);
x[k].real = (double)1.0*temp;
x[k].imag = 0;
}

WN = 2*PI/BUFFER;

for(k=0;k<BUFFER;++k)
{
XR[k]=0.0;
XI[k]=0.0;
wk = k*WN;
for(n=0;n<BUFFER;++n)
{
c = cos((double)n*wk);
s = sin((double)n*wk);
XR[k]=XR[k]+x[n].real*c+x[n].imag*s;
XI[k]=XI[k]-x[n].real*s+x[n].imag*c;
if(n==(BUFFER-10))
temp =0;
}
}
for(k=0;k<BUFFER;++k)
{
x[k].real=XR[k];
x[k].imag=XI[k];
x[k].modulus =sqrt(x[k].real*x[k].real+x[k].imag*x[k].imag);
x[k].angle = atan(x[k].imag/x[k].real);
}
for(k=0;k<BUFFER/2;k++)
{
fprintf(SPECTR, "%d\t%f\n", k, x[k].modulus/BUFFER);
}
free(x);
free(XR);
free(XI);

}

[dosikus]

Repytw, уже готовое но перенос с авр, со всеми ляпами/костылями/соплями http://we.easyelectronics.ru/reptile/8- ... 030f4.html

[VladislavS]

Repytw писал(а):

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

Если мне не изменяет память, одна из классический лаб по матану на 2-м курсе это запрограммировать ДПФ и БПФ, измерить разницу в скорости вычисления. Всё уже изобретено до нас.

[Repytw]

Если эта тема еще кому-то интересна.
Пример реализации 128-точечного ДПФ для 12-битных данных (по задумке от АЦП). Без float и double.
Пример опять-таки для ПК, но перенести его для приложений на ARM не составляет проблемы.

Спойлер

#include <stdio.h>

#define BUFFER 128 /*ДЛИНА БУФЕРА ДАННЫХ*/

int XR[BUFFER]; /*МАССИВ РЕЗУЛЬТАТОВ ДЕЙСТВИТЕЛЬНЫХ СОСТАВЛЯЮЩИХ ДПФ*/
int XI[BUFFER]; /*МАССИВ РЕЗУЛЬТАТОВ КОМПЛЕКСНЫХ СОСТАВЛЯЮЩИХ ДПФ*/

int RESULT[BUFFER]; /*МАССИВ РЕЗУЛЬТАТОВ ОПРЕДЕЛЕНИЯ АМПЛИТУДНОГО СПЕКТРА*/

/***МАССИВ 12-БИТНЫХ ВХОДНЫХ ДАННЫХ КОТОРЫЕ МЫ ЯКОБЫ ПОЛУЧАЕМ ОТ АЦП***/
int DATA[] = {
4096, 7502, 7880, 4895, 1200, 79, 2529, 6372,
8192, 6372, 2529, 79, 1200, 4895, 7880, 7502,
4096, 690, 312, 3297, 6992, 8113, 5664, 1820,
0, 1820, 5664, 8113, 6992, 3297, 312, 690,
4096, 7502, 7880, 4895, 1200, 79, 2529, 6372,
8192, 6372, 2529, 79, 1200, 4895, 7880, 7502,
4096, 690, 312, 3297, 6992, 8113, 5664, 1820,
0, 1820, 5664, 8113, 6992, 3297, 312, 690,
4096, 7502, 7880, 4895, 1200, 79, 2529, 6372,
8192, 6372, 2529, 79, 1200, 4895, 7880, 7502,
4096, 690, 312, 3297, 6992, 8113, 5664, 1820,
0, 1820, 5664, 8113, 6992, 3297, 312, 690,
4096, 7502, 7880, 4895, 1200, 79, 2529, 6372,
8192, 6372, 2529, 79, 1200, 4895, 7880, 7502,
4096, 690, 312, 3297, 6992, 8113, 5664, 1820,
0, 1820, 5664, 8113, 6992, 3297, 312, 690
};

/*******МАССИВ 12-БИТНЫХ КОСИНУСОВ ДЛЯ 128 ТОЧЕЧНОГО ДПФ*********/
/*****ЭТО - КОНСТАНТЫ, ПОЭТОМУ ИХ ЛУЧШЕ БЫ СОХРАНИТЬ ВО FLASH****/
signed int COS[] = {
2048, 2046, 2038, 2026, 2009, 1987, 1960, 1928,
1892, 1851, 1806, 1757, 1703, 1645, 1583, 1518,
1448, 1375, 1299, 1220, 1138, 1053, 965, 876,
784, 690, 595, 498, 400, 301, 201, 100,
0, -100, -201, -301, -400, -498, -595, -690,
-784, -876, -965, -1053, -1138, -1220, -1299, -1375,
-1448, -1518, -1583, -1645, -1703, -1757, -1806, -1851,
-1892, -1928, -1960, -1987, -2009, -2026, -2038, -2046,
-2048, -2046, -2038, -2026, -2009, -1987, -1960, -1928,
-1892, -1851, -1806, -1757, -1703, -1645, -1583, -1518,
-1448, -1375, -1299, -1220, -1138, -1053, -965, -876,
-784, -690, -595, -498, -400, -301, -201, -100,
0, 100, 201, 301, 400, 498, 595, 690,
784, 876, 965, 1053, 1138, 1220, 1299, 1375,
1448, 1518, 1583, 1645, 1703, 1757, 1806, 1851,
1892, 1928, 1960, 1987, 2009, 2026, 2038, 2046
};

/*******МАССИВ 12-БИТНЫХ СИНУСОВ ДЛЯ 128 ТОЧЕЧНОГО ДПФ***********/
/*****ЭТО - КОНСТАНТЫ, ПОЭТОМУ ИХ ЛУЧШЕ БЫ СОХРАНИТЬ ВО FLASH****/
signed int SIN[] = {
0, 100, 201, 301, 400, 498, 595, 690,
784, 876, 965, 1053, 1138, 1220, 1299, 1375,
1448, 1518, 1583, 1645, 1703, 1757, 1806, 1851,
1892, 1928, 1960, 1987, 2009, 2026, 2038, 2046,
2048, 2046, 2038, 2026, 2009, 1987, 1960, 1928,
1892, 1851, 1806, 1757, 1703, 1645, 1583, 1518,
1448, 1375, 1299, 1220, 1138, 1053, 965, 876,
784, 690, 595, 498, 400, 301, 201, 100,
0, -100, -201, -301, -400, -498, -595, -690,
-784, -876, -965, -1053, -1138, -1220, -1299, -1375,
-1448, -1518, -1583, -1645, -1703, -1757, -1806, -1851,
-1892, -1928, -1960, -1987, -2009, -2026, -2038, -2046,
-2048, -2046, -2038, -2026, -2009, -1987, -1960, -1928,
-1892, -1851, -1806, -1757, -1703, -1645, -1583, -1518,
-1448, -1375, -1299, -1220, -1138, -1053, -965, -876,
-784, -690, -595, -498, -400, -301, -201, -100
};

/***********ПОДПРОГРАММА РАСЧЕТА ДПФ**************/
/*АРГУМЕНТЫ: 1 - указатель на массив входных данных
2 - указатель на массив результатов расчёта*/
void DFT(int *INPUT, int *OUTPUT)
{
/*****ОБЪЯВЛЯЕМ КУЧКУ НУЖНЫХ ПЕРЕМЕННЫХ*****/
short n;
short k;
short temp;
signed int c;
signed int s;

/*****ВНЕШНИЙ ЦИКЛ ВЫЧИСЛЕНИЙ**************/
for(k=0; k<BUFFER/2; ++k)
{
XR[k] = 0;
XI[k] = 0;
for(n = 0; n<BUFFER; ++n)
{
temp = n*k;
temp %= BUFFER;
c = COS[temp];
s = SIN[temp];
XR[k] += INPUT[n]*c;
XI[k] -= INPUT[n]*s;
}
/*ОПРЕДЕЛЯЕМ ДИСКРЕТНЫЙ СПЕКТР АМПЛИТУД ИССЛЕДУЕМОГО СИГНАЛА*/
XR[k] >>= 16;
XI[k] >>= 16;
OUTPUT[k] = XR[k]*XR[k] + XI[k]*XI[k];
/*ЗДЕСЬ БЫ ПО УМУ КОРЕНЬ OUTPUT ПОСЧИТАТЬ, ТОЛЬКО ЭТО МНЕ НЕ
ПОД СИЛУ БЕЗ ЗНАЧИТЕЛЬНЫХ ЗАТРАТ РЕСУРСОВ */
}
}
/*************КОНЕЦ ПОДПРОГРАММЫ РАСЧЁТА ДПФ***************/

int main(void)
{
int k;

DFT(DATA, RESULT);

for(k=0;k<BUFFER/2;k++)
{
printf("%d\t%d\n", k, RESULT[k]);
}
_getch();
return 0;
}

[VladislavS]

Да вы издеваетесь? Кто ж ДПФ в лоб считает?

[Repytw]

Да вы издеваетесь?

Нет.

Кто ж ДПФ в лоб считает?

Те, кто хотят понять алгоритм.

Всё уже изобретено до нас.

Ждём Ваших усовершенствований)))

[dosikus]

Я это не писал...

[Repytw]

Я это не писал...

Это случайно вышло. Редактор косячит.

Всё уже изобретено до нас

[YS]

Те, кто хотят понять алгоритм.

Эээ, весь вопрос в этой теме (как я его понимаю) - в эффективной реализации быстрого преобразования Фурье.

Алгоритм обычного дискретного преобразования Фурье элементарен, там и понимать особо нечего - это просто расчет скалярных произведений входного сигнала на набор опорных частот, для каждой в двух вариантах, сдвинутых на 90° (синус и косинус, чтобы знать фазу).

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

Так что с обычным ДПФ все просто. А вот как реализуются ускоренные алгоритмы - это неочевидно и предмет для дискуссий.

[12val12]

del