ivan219 0 16 ноября, 2008 Опубликовано 16 ноября, 2008 · Жалоба Уковонибудь есть алгоритм вычесления БПФ по основанию 4. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Связист 0 16 ноября, 2008 Опубликовано 16 ноября, 2008 · Жалоба Уковонибудь есть алгоритм вычесления БПФ по основанию 4. Эти алгоритмы (как и куча других быстрых алгоритмов преобразований Фурье) описаны у Блейхута в книге "Быстрые алгоритмы цифровой обработки сигналов", Москва, Мир, 1989. Очень распространённая книга. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
ivan219 0 24 ноября, 2008 Опубликовано 24 ноября, 2008 · Жалоба А в исходниках есть??? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Связист 0 29 ноября, 2008 Опубликовано 29 ноября, 2008 · Жалоба А в исходниках есть??? Есть. У Texas Instruments. Искать надо в Code Composer Studio исходники библиотеки dsplib. Они там кажется уже оптимизированы на ассемблере. И ещё возможно в матлабе есть. Но точно не скажу где именно. Кажется там где поддержка процессоров TI в симулинке. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
shasik 0 30 ноября, 2008 Опубликовано 30 ноября, 2008 · Жалоба А в исходниках есть??? Эх, да простят меня модераторы ибо уже раз в четвертый рекламирую этот источник: Алгоритмы для программеров Там есть куча алгоритмов, в том числе и FFT4. Ищите там книгу (она есть в нескольких форматах, лично я скачивал pdf) и архив с исходниками к книге. Ну, а если не найдете то вот оно здесь: Algorithms_for_Programmers.rar sources.rar ЗЫЖ Можно было бы, конечно, вырезать только нужно, но вдруг вас что-нибудь еще заинтересует, поэтому полностью. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
diwil 0 10 декабря, 2008 Опубликовано 10 декабря, 2008 · Жалоба вот в исходниках. Правда это комбинированое 4-2, т.е. если степень двойки кратна 2м, то чистое ффт по основанию 4, а если не кратно, то сначала 4ка, а последний шаг по основанию 2. Это обратное фурье - экспонента положительна. Чтобы сделать отрицательную экспоненту, надо поменять знаки в операциях умножения на j и знак мнимой части коэффициентов. Масштабирования нет. Чтобы ввести масштабирование на 1/N надо поменять значение SHIFT_AMOUNT_UNSCALED Код писался мной для себя, поэтому комментов не много, уж не обессудте. static void data_swap(fft_inter data[], int np) { fft_inter *pxk, *pxi; fft_inter *px = data; int tmp; int n2, n21; int i, j, k; n2 = np >> 1; /// n2: 512 n21 = np + 1; /// n21: 1025 j = 0; for (i = 0; i < np; i += 4) { /// i,j: 0,0; 4,1024; 8,512; 12,1536; ... 2044,??? [255x] if (i < j) { // swap 4 pairs of values (2 complex pairs with 2 others) // px[i] <-> px[j]; px[i+1] <-> px[j+1] // px[i+1+n21] <-> px[j+1+n21]; px[i+1+n21+1] <-> px[j+1+n21+1] pxi = &px[i]; pxk = &px[j]; tmp = *pxi; *pxi++ = *pxk; *pxk++ = tmp; tmp = *pxi; *pxi = *pxk; *pxk = tmp; pxi += n21; pxk += n21; tmp = *pxi; *pxi++ = *pxk; *pxk++ = tmp; tmp = *pxi; *pxi = *pxk; *pxk = tmp; } // swap 2 pairs of values (1 complex pair with another) // px[i+2] <-> px[j+np]; px[i+3] <-> px[j+np+1] pxi = &px[i + 2]; pxk = &px[j + np]; tmp = *pxi; *pxi++ = *pxk; *pxk++ = tmp; tmp = *pxi; *pxi = *pxk; *pxk = tmp; k = n2; while (k <= j) { /// k: {1024} {1024,512} {1024} {1024,512,256} ... j -= k; k = k >> 1; } j += k; /// j: {1024} {512} {1536} {256} ... } } #define COMPLEX_MUL(a,b,c,d)\ do{ fft_extend x,y; x = MULT16(a,c)-MULT16(b,d); y = MULT16(c,b)+MULT16(a,d); a = x; b = y; }while(0) #define SHIFT_AMOUNT_UNSCALED 0 void radix4_unscaled(struct complex_s *data, int size, fft_inter *tw) { struct complex_s *x = data; int i, l; int N; struct complex_s32 x0,x1,x2,x3, t1,t2,t3,t4,t; fft_inter wre,wim; N = size; while(N > 4) { for (l = 0; l< size; l+= N) { for (i = l; i < l + N/4; i++) { int idx = i; x0.re = x[i].re; x0.im = x[i].im; x2.re = x[N/2+i].re; x2.im = x[N/2+i].im; t1.re = x0.re + x2.re; t1.im = x0.im + x2.im; t2.re = x0.re - x2.re; t2.im = x0.im - x2.im; x0.re = x[i+N/4].re; x0.im = x[i+N/4].im; x2.re = x[N/2+i+N/4].re; x2.im = x[N/2+i+N/4].im; t3.re = x0.re + x2.re; t3.im = x0.im + x2.im; t4.re = x0.re - x2.re; t4.im = x0.im - x2.im; // update x_{0+i} t.re = t1.re + t3.re; t.im = t1.im + t3.im; x[i].re = t.re>>SHIFT_AMOUNT_UNSCALED; x[i].im = t.im>>SHIFT_AMOUNT_UNSCALED; // x_{0+i} updated // update x_{N/4+i} // wre = FR2(cos(2.* M_PI / N * 2. * idx/1)); // wim = FR2(sin(2.* M_PI / N * 2. * idx/1)); //printf("/*N:%2d,i:%2d */ FR2(%+.16f), FR2(%+.16f),\n", N,idx, // cos(2.* M_PI / N * 2. * idx/1), sin(2.* M_PI / N * 2. * idx/1)); wre = *tw++; wim = *tw++; t.re = t1.re - t3.re; t.im = t1.im - t3.im; COMPLEX_MUL(t.re,t.im, wre, wim); x[i+N/4].re = t.re>>SHIFT_AMOUNT_UNSCALED; x[i+N/4].im = t.im>>SHIFT_AMOUNT_UNSCALED; // x_{N/4+i} updated // update x_{N/2} // wre = FR2(cos(2.* M_PI / N * 1. * idx/1)); // wim = FR2(sin(2.* M_PI / N * 1. * idx/1)); //printf("/*N:%2d,i:%2d */ FR2(%+.16f), FR2(%+.16f),\n", N,idx, // cos(2.* M_PI / N * 1. * idx/1), sin(2.* M_PI / N * 1. * idx/1)); wre = *tw++; wim = *tw++; t.re = t2.re - t4.im; /// ??? t.im = t2.im + t4.re; /// ??? COMPLEX_MUL(t.re,t.im, wre, wim); x[N/2+i].re = t.re>>SHIFT_AMOUNT_UNSCALED; x[N/2+i].im = t.im>>SHIFT_AMOUNT_UNSCALED; // x_{N/2} updated // update x_{3N/4+i} // wre = FR2(cos(2.* M_PI / N * 3. * idx/1)); // wim = FR2(sin(2.* M_PI / N * 3. * idx/1)); //printf("/*N:%2d,i:%2d */ FR2(%+.16f), FR2(%+.16f),\n", N,idx, // cos(2.* M_PI / N * 3. * idx/1), sin(2.* M_PI / N * 3. * idx/1)); wre = *tw++; wim = *tw++; t.re = t2.re + t4.im; /// ??? t.im = t2.im - t4.re; /// ??? COMPLEX_MUL(t.re,t.im, wre, wim); x[N/2+i+N/4].re = t.re>>SHIFT_AMOUNT_UNSCALED; x[N/2+i+N/4].im = t.im>>SHIFT_AMOUNT_UNSCALED; // x_{3N/4+i} updated } } N >>= 2; } if ( N != 2 ) // FIXME: N == 4 { for (i = 0; i < size; i+= 4) { x0.re = x[i].re; x0.im = x[i].im; x2.re = x[2+i].re; x2.im = x[2+i].im; t1.re = x0.re + x2.re; t1.im = x0.im + x2.im; t2.re = x0.re - x2.re; t2.im = x0.im - x2.im; x0.re = x[i+1].re; x0.im = x[i+1].im; x2.re = x[3+i].re; x2.im = x[3+i].im; t3.re = x0.re + x2.re; t3.im = x0.im + x2.im; t4.re = x0.re - x2.re; t4.im = x0.im - x2.im; // update x_{0+i} t.re = t1.re + t3.re; t.im = t1.im + t3.im; x[i].re = t.re>>(SHIFT_AMOUNT_UNSCALED/1); x[i].im = t.im>>(SHIFT_AMOUNT_UNSCALED/1); // x_{0+i} updated // update x_{N/4+i} t.re = t1.re - t3.re; t.im = t1.im - t3.im; x[i+1].re = t.re>>(SHIFT_AMOUNT_UNSCALED/1); x[i+1].im = t.im>>(SHIFT_AMOUNT_UNSCALED/1); // x_{N/4+i} updated // update x_{N/2} t.re = t2.re - t4.im; /// ??? t.im = t2.im + t4.re; /// ??? x[2+i].re = t.re>>(SHIFT_AMOUNT_UNSCALED/1); x[2+i].im = t.im>>(SHIFT_AMOUNT_UNSCALED/1); // x_{N/2} updated // update x_{3N/4+i} t.re = t2.re + t4.im; /// ??? t.im = t2.im - t4.re; /// ??? x[3+i].re = t.re>>(SHIFT_AMOUNT_UNSCALED/1); x[3+i].im = t.im>>(SHIFT_AMOUNT_UNSCALED/1); // x_{3N/4+i} updated } } else { // trivial butts at the end for (i = 0; i<size; i+= 2) { x0.re = x[i].re; x0.im = x[i].im; x1.re = x[i+1].re; x1.im = x[i+1].im; t.re = x0.re + x1.re; t.im = x0.im + x1.im; x[i].re = t.re>>1; x[i].im = t.im>>1; t.re = x0.re - x1.re; t.im = x0.im - x1.im; x[i+1].re = t.re>>1; x[i+1].im = t.im>>1; } } data_swap((fft_inter *)data,size); } Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
ivan219 0 12 декабря, 2008 Опубликовано 12 декабря, 2008 (изменено) · Жалоба Спасибо конечно но я в С++ не силён :( может на Пскале есть??? Изменено 12 декабря, 2008 пользователем ivan219 Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
vadkudr 0 16 декабря, 2008 Опубликовано 16 декабря, 2008 · Жалоба http://www.fftw.org/newsplit.pdf Вот здесь есть алгоритм еще быстрее, чем сплит-радикс 4. Проверял на матлабе - работает. И на сам сайт интересно посмотреть. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
diwil 0 17 декабря, 2008 Опубликовано 17 декабря, 2008 · Жалоба http://www.fftw.org/newsplit.pdf Вот здесь есть алгоритм еще быстрее, чем сплит-радикс 4. Проверял на матлабе - работает. И на сам сайт интересно посмотреть. это сплит.... Теоретически он работает несколько быстреев смысле требуемого числа операций. Однако при его реализации оказывается, что помимо выполнения основных операций сложения-умножения, нужно еще крутить большое количество счетчиков. Это можно хорошо делать на жесткой логике (типа ПЛМ). Однако при реализации на процессоре (-рах) большое количество счетчиков не позволяет разместить все промежуточные данные в регистрах и приходится задействовать или стэк или память, что может привести к заметному замедлению вычисления. Не могу говорить про все ДСП... да и вообще про ДСП :) но, скажем, на АРМе радикс4 будет самым быстрым. Хотя, на код матлаба из статьи было бы посмотреть оч интересно. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
shasik 0 17 декабря, 2008 Опубликовано 17 декабря, 2008 · Жалоба Теоретически он работает несколько быстреев... А вдруг у автора топика сигнал полностью реальный (нет комплексной части), тогда можно предложить пару-тройку десятков еще более быстрых алгоритмов. И т.д., и т.п. Все это пустое. Убежден, что автор топика - студент, который получил задание и решил, что он самый хитрый. Вас не настораживают его посты: Уковонибудь есть алгоритм вычесления БПФ по основанию 4 А в исходниках есть??? может на Пскале есть???На мой взгляд все ясно. Пока мы здесь напрягаемся, студент уже давно где-нибудь на форуме программистов: ищет исходники для Delphi. Здесь же ему не помогли. ЗЫЖ Извини, студент, если обидел. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
ivan219 0 28 декабря, 2008 Опубликовано 28 декабря, 2008 · Жалоба ЗЫЖ Извини, студент, если обидел. Не необиделся. А исходники на Delphi не нашол :( Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
DRUID3 0 28 декабря, 2008 Опубликовано 28 декабря, 2008 · Жалоба Не необиделся. А исходники на Delphi не нашол :( как насчет варианта переделать руками из-под C-шного? И еще вопрос - а зачем это Вам? Не успевает отработать между прерываниями? по основанию 4 это в лучшем случае 25% выигрыша по сравнению чем по основанию 2... За что это уже "делфисты" беруЦЦо такое... Если взять RADIX2, написать на asm собрать в библиотеку и вызывать из нее выигрыш будет возможно даже выше... А теперь вопрос - а нужно ли это делать если уже есть готовые оптимизированные и протестированные под x86 аппаратную платформу... Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
vadkudr 0 29 декабря, 2008 Опубликовано 29 декабря, 2008 · Жалоба Хотя, на код матлаба из статьи было бы посмотреть оч интересно. Посмотрите - писался просто для проверки теории в статье. Поэтому заведомо работает медленно. Запуск: >> a=randn(1,32) a = Columns 1 through 8 -0.4326 -1.6656 0.1253 0.2877 -1.1465 1.1909 1.1892 -0.0376 Columns 9 through 16 0.3273 0.1746 -0.1867 0.7258 -0.5883 2.1832 -0.1364 0.1139 Columns 17 through 24 1.0668 0.0593 -0.0956 -0.8323 0.2944 -1.3362 0.7143 1.6236 Columns 25 through 32 -0.6918 0.8580 1.2540 -1.5937 -1.4410 0.5711 -0.3999 0.6900 >> b=Alg_02_NewSplitRadixFFT(a) b = Columns 1 through 4 2.8652 -4.0386 - 2.6955i -3.3874 + 0.9221i -2.1831 - 3.6784i Columns 5 through 8 3.5893 + 5.2095i 0.9802 + 5.1840i -0.8291 + 0.9154i -2.6340 + 2.6884i Columns 9 through 12 -5.0758 - 1.0582i -4.9958 + 0.0682i 7.7442 + 0.5439i 2.8367 + 9.6884i Columns 13 through 16 2.7128 + 4.6691i -1.1440 - 4.9140i 0.4670 + 5.2595i -0.8161 - 2.9031i Columns 17 through 20 -3.1601 -0.8161 + 2.9031i 0.4670 - 5.2595i -1.1440 + 4.9140i Columns 21 through 24 2.7128 - 4.6691i 2.8367 - 9.6884i 7.7442 - 0.5439i -4.9958 - 0.0682i Columns 25 through 28 -5.0758 + 1.0582i -2.6340 - 2.6884i -0.8291 - 0.9154i 0.9802 - 5.1840i Columns 29 through 32 3.5893 - 5.2095i -2.1831 + 3.6784i -3.3874 - 0.9221i -4.0386 + 2.6955i >> c=fft(a) c = Columns 1 through 4 2.8652 -4.0386 - 2.6955i -3.3874 + 0.9221i -2.1831 - 3.6784i Columns 5 through 8 3.5893 + 5.2095i 0.9802 + 5.1840i -0.8291 + 0.9154i -2.6340 + 2.6884i Columns 9 through 12 -5.0758 - 1.0582i -4.9958 + 0.0682i 7.7442 + 0.5439i 2.8367 + 9.6884i Columns 13 through 16 2.7128 + 4.6691i -1.1440 - 4.9140i 0.4670 + 5.2595i -0.8161 - 2.9031i Columns 17 through 20 -3.1601 -0.8161 + 2.9031i 0.4670 - 5.2595i -1.1440 + 4.9140i Columns 21 through 24 2.7128 - 4.6691i 2.8367 - 9.6884i 7.7442 - 0.5439i -4.9958 - 0.0682i Columns 25 through 28 -5.0758 + 1.0582i -2.6340 - 2.6884i -0.8291 - 0.9154i 0.9802 - 5.1840i Columns 29 through 32 3.5893 - 5.2095i -2.1831 + 3.6784i -3.3874 - 0.9221i -4.0386 + 2.6955i >> NewSplitRadix.zip Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
diwil 0 31 декабря, 2008 Опубликовано 31 декабря, 2008 · Жалоба о! пасиб! на досуге попробую посмотреть, что можно сделать. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
