sigmaN 0 21 сентября, 2009 Опубликовано 21 сентября, 2009 · Жалоба bahurin, БОЛЬШОЕ Вам человеческое спасибо за вразумление :) Вопрос решен в таком виде: //это прямое преобразование RFFT_F32_STRUCT *t=(RFFT_F32_STRUCT*)table; scale = 1./t->FFTSize; t->InBuf=in; t->OutBuf=out; RFFT_f32u(t); //делаем БПФ for(i=0;i<t->FFTSize;i++) t->OutBuf[i] *= scale; //ОБРАТНОЕ //готовимся к обратному преобразованию //Re[i] содержит вещественную часть Im[i] - мнимую //отразим мнимую часть и дополним где надо нулями Re[0] = in[0]; Im[0] = 0; for (i=1; i<(t->FFTSize/2+1);i++){ Re[i] = in[i]; Im[i] = in[t->FFTSize-i]; } Im[t->FFTSize/2] = 0; //теперь восстанавливаем вторую половину спектра // ************* mirror *************** for (i=0; i<(t->FFTSize/2);i++) { Re[t->FFTSize-1-i] = Re[i+1]; //Re просто зеркально Im[t->FFTSize-1-i] = - Im[i+1]; //Im зеркально с минусом } //Сложить for (i=0; i<(t->FFTSize);i++){ t->InBuf[i]=Re[i]+Im[i]; } RFFT_f32u(t); //делаем БПФ, теперь в OutBuf имеем time domain Re[0] = t->OutBuf[0]; Im[0] = 0; for (i=1; i<(t->FFTSize/2+1);i++){ Re[i] = t->OutBuf[i]; Im[i] = t->OutBuf[t->FFTSize-i]; } Im[t->FFTSize/2] = 0; // ************* mirror *************** for (i=0; i<(t->FFTSize/2);i++) { Re[t->FFTSize-1-i] = Re[i+1]; //Re просто зеркально Im[t->FFTSize-1-i] = - Im[i+1]; //Im зеркально с минусом } //Сложить for (i=0; i<(t->FFTSize);i++){ out[i]=Re[i]+Im[i]; } //на FFTSize делить не нужно, т.к. это делает функция прямого преобразования Единственное, что не даёт покоя - это кусок кода парней с Техасовского форума ---combine and split imag/real--- for(i=0; i< FFTsize/2;i++) { tempR = Re[i] - Im[i]; tempI = Re[i]] + Im[i]; Re[i] = tempR; //can also be done in one step Im[i] = tempI; } Зачем, если работает и без этого!?!?! Вот полный пример You can use the Real FFT to do IFFT, but you have to the following steps --Split-- Re[0] = fft.OutBuf[0]; Im[0] = 0; for (i=1; i<(FFTsize/2+1);i++) { Re[i] = fft.OutBuf[i]; // split real Im[i] = fft.OutBuf[FFTsize-i]; // split imag } Im[FFTsize/2] = 0; ---combine and split imag/real--- //ЧТО ЭТО И ЗАЧЕМ for(i=0; i< FFTsize/2;i++) { tempR = Re[i] - Im[i]; tempI = Re[i]] + Im[i]; Re[i] = tempR; //can also be done in one step Im[i] = tempI; } --mirror--- // ************* mirror *************** for (i=0; i<(FFTsize/2);i++) { Re[FFTsize-1-i] = Re[i+1]; Im[FFTsize-1-i] = - Im[i+1]; } // ************ add ******************* //А ТУТ "-" ВМЕСТО "+" for (i=0; i<FFTsize;i++) { fft.InBuf[i] = Re[i] - Im[i]; } And now you can do the FFT again to get the IFFT After this, you will need to the Split, Mirror routine again The last step is to divide the all the result by [fftSize] А вообще всё хорошо. Завтра на асме все эти отражения реализую и будет супер :) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
andrewn 0 21 сентября, 2009 Опубликовано 21 сентября, 2009 (изменено) · Жалоба Вопрос решен в таком виде: Сложно. Можно сделать проще. Я не стал связываться с кодом для С283хх, а использовал свою подпрограмму для С67хх, но всё так просто, что копируется в структуру 1 в 1. Этот способ - чистой воды overkill и упражнение в алгебре. Для проверки правильности я использовал rffti(), обратное преобразование, которое работает напрямую с выходом прямого. // --------------------------------------------------------------------------- // // InvRFFT() - simulates inverse real FFT // // all input/output arrays are of length 2**log2n real entries // // on input // // log2n - base 2 log of Real FFT length // w - twiddle factors table // x - input array for Inverse Real FFT (stored in the order // of Forward Real FFT output array) // f - scratch array // // on input // // s - output array of inverse coeffs stored in natural order // // notes // // this is a math exercise, not suitable for real work // // --------------------------------------------------------------------------- #define real float void InvRFFT (int log2n, const real *w, real *x, real *f, real *s) { int i, n, n2; real invn; n = 1 << log2n; // RFFT length n2 = n >> 1; // half length invn = 1.0E0F / (real)(n); // 1/N scale factor // add real and imag parts f[0] = x[0]; for (i = 1; i < n2; i++) { f[i] = x[i] + x[n-i]; // add real and mirror imag f[n-i] = x[i] - x[n-i]; // add mirror real and conj imag } f[n2] = x[n2]; rfftf (log2n, (real *)w, (real *)f); // FORWARD RFFT // add real and imag parts and scale down to form inverse real fft // get rid of scaling if done inside forward real fft s[0] = f[0] * invn; // save a few clocks on division for (i = 1; i < n2; i++) { s[i] = (f[i] + f[n-i]) * invn; s[n-i] = (f[i] - f[n-i]) * invn; } s[n2] = f[n2] * invn; return; // all done, exit } Изменено 21 сентября, 2009 пользователем AndrewN Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
sigmaN 0 22 сентября, 2009 Опубликовано 22 сентября, 2009 · Жалоба Даа, я естественно тоже пришел к выводу, что там возня лишняя присутствует. И ресурсы как памяти так и проца расходуются непонятно куда. t->InBuf - входной буфер RealFFT in - комплексный вывод того самого RealFFT in в формате принятом в Speex. А именно R[0],R[1],I[1],R[2],I[2]....R[N/2] t->InBuf[0]=in[0]; for (i=1, j=1; i<(t->FFTSize/2); i++){ t->InBuf[i]=in[j]+in[j+1]; j+=2; } t->InBuf[t->FFTSize/2]=in[t->FFTSize-1]; for (i=t->FFTSize/2+1, j=(t->FFTSize-3); i<(t->FFTSize); i++){ t->InBuf[i]=in[j]-in[j+1]; j-=2; } А это из техасовского комплексного вывода делает вещественный. Как последний шаг моей IFFT. t->OutBuf - комплексный вывод техасовского RFFT в формате R[0],R[1],R[2]....R[N/2],I[N/2-1],I[N/2-2]....I[1] out - вывод функции, где в итоге получается восстановленый сигнал(in time domain) out[0]=t->OutBuf[0]; for (i=1; i<(t->FFTSize/2);i++) out[i]=t->OutBuf[i] + t->OutBuf[t->FFTSize-i]; out[i++]=t->OutBuf[i]; for (; i<(t->FFTSize);i++) out[i]=t->OutBuf[t->FFTSize-i] - t->OutBuf[i]; Функция "перепаковывающая" техасовский вывод в формат, принятый в Speex ;static void _TMStoSPEEXrepack(spx_word16_t *SPXpacked, spx_word16_t *TMSpacked, int N){ ;//перепаковывает рузультаты FFT в формат принятый в speex ;//TMSpacked[0] = real[0] TMSpacked[1] = real[1].... TMSpacked[N/2] = real[N/2] ;//TMSpacked[N-1] = imag[1] TMSpacked[N-2] = imag[2].... TMSpacked[N/2+1] = imag[N/2-1] ;//SPXpacked = R(0),R(1),I(1),R(2),I(2)....R(N/2) ; ;SPXpacked XAR4 ;TMSpacked XAR5 ;N AL .global __TMStoSPEEXrepack .sect "ramfuncs" __TMStoSPEEXrepack: PUSH AL ;расчитаем указатель на imag[N/2-1] MOVL ACC,XAR5 ADD ACC,*-SP[1] ADD ACC,*-SP[1];второй раз потому что float занимает 2слова MOVL XAR6,ACC;XAR6 указывает на TMSpacked[N-1] т.е. на imag[1] ;подготовить счётчик цикла POP AL LSR AL#1 ADD AL,#-2;цикл должен пройти N/2-1 раз ;до начала цикла делаем ;SPXpacked[0]=TMSpacked[0]; MOV32 R0H,*XAR5++ MOV32 *XAR4++,R0H .align 2 NOP RPTB fftrepack_loop,AL MOV32 R0H,*XAR5++;real[i] MOV32 *XAR4++,R0H MOV32 R0H,*--XAR6;imag[i] MOV32 *XAR4++,R0H fftrepack_loop: MOV32 R0H,*XAR5++;real[i] MOV32 *XAR4++,R0H LRETR А это чтбы быстро перемножить буфер на множитель, попутно скопировав результат в другой буфер. Также работает "на месте" - т.е. с одним буфером. ;void asm_mul_by_const(float* dst, float* src, float mul, int len); ;dst XAR4 ;src XAR5 ;mul R0H ;len AL ;умножает каждый элемент src на mul результат ложится в dst ;dst[i]=src[i]*mul .global _asm_mul_by_const .sect "ramfuncs" _asm_mul_by_const: LSR AL,#1;цикл unrolled, так что делим счётчик на 2 ADDB AL,#-1; чтобы RPTB не перестарался RPTB asm_mul_by_const_loop,AL MOV32 R1H,*XAR5++;R1H=src[i] MPYF32 R1H,R1H,R0H;R1H=src[i]*mul MOV32 R2H,*XAR5++;R2H=src[i+1] MPYF32 R2H,R2H,R0H;R2H=src[i+1]*mul MOV32 *XAR4++,R1H;dst[i]=R1H MOV32 *XAR4++,R2H;dst[i+1]=R1H asm_mul_by_const_loop: LRETR Остались эти speexComplex_to_real() и TMScomplex_to_real() оптимизировать. На сях они кушают непомерно много даже со всеми оптимизациями компилятора и прочими премудростями. Обидно, когда тупое копирование буферов занимает по тактам половину времени расчёта FFT :) Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
sigmaN 0 23 сентября, 2009 Опубликовано 23 сентября, 2009 · Жалоба Воистину ассемблер и трпение - великая сила! В итоге Speex echo canceller при праметрах ECHO_FRAME_LEN = 32 filter tail size = 128 удалось вписать в 21MIPS! Это ещё осталось 2 кандидата на оптимизацию. Думаю после их оптимизации можно будет вписаться в 17-18, но это уже больше ради спортивного интереса.... Даже не верится, что до начала оптимизации с этими параметрами получалось около 70MIPS Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
