c28х FPU library beta1, помогите разобраться с real FFT Опции

[sigmaN]

Имеется оптимизированная либа C28x Floating Point Unit Library, содержащая следующие функции

RFFT_f32 void RFFT_f32(RFFT_F32_STRUCT *);
This module computes a 32-bit floating-point real FFT including input bit
reversing. This version of the function requires input buffer memory alignment. If
you do not wish to align the input buffer, then use the RFFT_f32u function.

RFFT_f32u void RFFT_f32u(RFFT_F32_STRUCT *);
This module computes a 32-bit floating-point real FFT including input bit
reversing. This version of the function does not have any buffer alignment
requirements. If you can align the input buffer, then use the RFFT_f32 function
for improved performance.

RFFT_f32_mag void RFFT_f32_mag(RFFT_F32_STRUCT *);
This module computes the real FFT magnitude.

RFFT_f32s_mag void RFFT_f32s_mag(RFFT_F32_STRUCT *);
This module computes the scaled real FFT magnitude.

RFFT_f32_phase void RFFT_f32_phase(RFFT_F32_STRUCT *);
This module computes FFT Phase

RFFT_f32_sincostable void RFFT_f32_sincostable(RFFT_F32_STRUCT *);
This module generates the twiddle factors used by the real FFT.


typedef struct {
float32 *InBuf;//Input buffer
float32 *OutBuf; //Output buffer
float32 *CosSinBuf;
float32 *MagBuf;
float32 *PhaseBuf;
Uint16 FFTSize;
Uint16 FFTStages;
} RFFT_F32_STRUCT;
По поводу output buffer сказано:
Result of RFFT_f32. This buffer is used as
the input to the magnitude and phase
calculations. The output order for FFTSize
= N is:
OutBuf[0] = real[0]
OutBuf[1] = real[1]
OutBuf[2] = real[2]
………
OutBuf[N/2] = real[N/2]
OutBuf[N/2+1] = imag[N/2-1]
………
OutBuf[N-3] = imag[3]
OutBuf[N-2] = imag[2]
OutBuf[N-1] = imag[1]

Почему на выходе _REAL_ FFT видим imag, т.е. мнимую часть?
Каков смысл FFT magnitude и phase?
Как используя эти функции реализовать обратное преобразование?

[AndrewN]

a. Почему на выходе _REAL_ FFT видим imag, т.е. мнимую часть?
b. Каков смысл FFT magnitude и phase?
c. Как используя эти функции реализовать обратное преобразование?

а. мнимые части появляются из-за комплексных матричных факторов
в. mag = sqrt (re**2 + im**2), phase = atan(im/re)
с. только специальной отдельной подпрограммой, у которой входные данные расположены в порядке выходных данных прямого преобразования.

[sigmaN]

Ответ по пункту c. подразумевает редактирование исходников либы или же имеется ввиду подпрограмма-обёртка, формирующая вход?

Мне кажется realFFT на то и real, что на входе у него только real, без мнимых частей. Тогда становится невозможным подать выход на вход и получить обратное преобразование. Следовательно, необходимо редактирование исходников либы. Что не очень-то и радует.

Или я чего-то не понимаю.

P.S. да нет, тут как ни крути, а придётся редактировать....
Но не верится, что люди выпустили оптимизированную либу без такой, казалось-бы очевидной и нужной вещи как обратное преобразование.

[des00]

Ответ по пункту c. подразумевает редактирование исходников либы или же имеется ввиду подпрограмма-обёртка, формирующая вход?

Мне кажется realFFT на то и real, что на входе у него только real, без мнимых частей. Тогда становится невозможным подать выход на вход и получить обратное преобразование. Следовательно, необходимо редактирование исходников либы. Что не очень-то и радует.

Или я чего-то не понимаю.

P.S. да нет, тут как ни крути, а придётся редактировать....
Но не верится, что люди выпустили оптимизированную либу без такой, казалось-бы очевидной и нужной вещи как обратное преобразование.

мой английский слаб, но мне кажеться что речь идет о реализации вычислений в формате с плавающей точкой при длине слова 32 бита. Насчет как посчитать IFFT через FFT смотрите http://www.dspsystem.narod.ru/content/fft/fft/fft.html

[sigmaN]

Английский у всех у нас слаб

Насчёт комплексного сопряжения я писал в другой ветке и получил вот какой ответ
Идея правильная, но, увы, абсолютно неприменимая в
случае вещественнозначного преобразования. Входные
данные такая процедура интерпретирует как строго
вещественные, следовательно они тождественно равны
своим сопряжённым. Поэтому и конструируют две подпрограммы,
для прямого и обратного преобразований.

[bahurin]

Почему на выходе _REAL_ FFT видим imag, т.е. мнимую часть?
Каков смысл FFT magnitude и phase?
Как используя эти функции реализовать обратное преобразование?

а есть какое нить описание остальных буферов:
float32 *InBuf;//Input buffer
float32 *CosSinBuf;
float32 *MagBuf;
float32 *PhaseBuf;

Для обратного преобразования надо чтобы в InBuf можно было положить как реальную так и мнимую части спектра.

[sigmaN]

Собственно говоря у меня тот-же вопрос и в сопроводительной доке ответа на него нет.

http://www.ti.com/litv/zip/sprc624 - либа вместе с докой

Вечерком ковырну исходники и посмотрю. Ну или эксперимент организую.....


Ну или накрайняк сапорт озадачить попробую......

дока отдельно:

Прикрепленные файлы

 C28x_FPU_Library_Beta1.pdf ( 112.01 килобайт ) Кол-во скачиваний: 733

 

[bahurin]

Собственно говоря у меня тот-же вопрос и в сопроводительной доке ответа на него нет.

http://www.ti.com/litv/zip/sprc624 - либа вместе с докой

Вечерком ковырну исходники и посмотрю. Ну или эксперимент организую.....


Ну или накрайняк сапорт озадачить попробую......

дока отдельно:

В OutBuf лежит N чисел, причем первые N/2+1 чисел а именно OutBuf[0:N/2] - реальные части половины БПФ, а вторые N/2-1 точек OutBuf[N/2+1:N-1] - мнимые части половины БПФ отраженные зеркально. Поскольку сигнал реальный, то вторая половина БПФ симметрична первой и ее можно не учитывать. Теперь если надо делать IFFT, то надо поступить так:
R = OutBuf[0:N/2] N/2+1 реальных частей
I = OutBuf[N/2+1:N-1] N/2-1 МНИМЫХ частей

R1 = [R fliplr(R(1:N/2-1))] формируем массив N точек путем симметричного отображения
I1 = [0 fliplr(I) 0 -I] формируем массив N точек путем симметричного отображения c инверсией
fliplr отражает массив слева направо

таким образом R1+j*S1 - комплексный спектр, т.е. точный результат FFT теперь надо вязть IFFT:
IFFT(R1+j*S1) = IFFT(R1)+j*IFFT(S1) = (RR-II)/N, где RR = real(FFT(R1)), т.е. реальная часть FFT(R1) а II = imag(FFT(I1)) - мнимая часть FFT(I1)

Таким образом надо взять БПФ от реальной и мнимой части спектра

вот пример в матлабе

Код

clear all
j = sqrt(-1);
N = 64;
s  = rand(1,N);
S = fft(s);
r = real(S(1:N/2+1)); % это реальная часть половины спектра в OutBuf
i = fliplr(imag(S(2:N/2))); % это мнимая часть половины спектра в OutBuf

r = [r fliplr(r(2:N/2))];% ОТРАЖАЮ зеркально
i = [0 fliplr(i) 0 -i]; % ОТРАЖАЮ зеркально с минусом


R = fft(r);
RR = real(R(1:N/2+1)); % это реальная часть половины спектра в OutBuf при fft(r)  
RI = fliplr(imag(R(2:N/2))); % это мнимая часть половины спектра в OutBuf при fft(r)


I = fft(-i);
IR = real(I(1:N/2+1)); % это реальная часть половины спектра в OutBuf при fft(i)  
II = fliplr(imag(I(2:N/2))); % это мнимая часть половины спектра в OutBuf при fft(i)

%отражаю вторые половины спектра:
A = [RR fliplr(RR(2:N/2))];% ОТРАЖАЮ зеркально
B = [0 fliplr(II) 0 -II]; % ОТРАЖАЮ зеркально с минусом

s0 = (A-B)/N;
plot(1:N,s,'-o',1:N,s0)

Сообщение отредактировал bahurin - Sep 17 2009, 10:58

[sigmaN]

Нашел пост на техасовском форуме

Код

There are two things why this won't work

-You will need a complex FFT to do the IFFT. Because your trying to place imaginaire data in the real FFT.
- You need mirror the imaginaire part.

(1)

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]

А вот отзыв человека это испробовавшего:

Код

Thanks, Leo. It's realy helpful.

I optimize you code for pre-IFFT computation and I need to change last "add" section for post-IFFT computation. And IFFT result is shifted for 3 samples rigth from source signal.

This is my code:

  RFFT_f32u(&fft1);
  
    for (int16 i=0; i<fftSize/2; i++){
        ifft.InBuf[i]=-2*fft1.OutBuf[fftSize-i];
        ifft.InBuf[fftSize-1-i]=2*fft1.OutBuf[i+1];
    }
    ifft.InBuf[0]=0;
    ifft.InBuf[fftSize/2]=fft1.OutBuf[fftSize/2];
      
    RFFT_f32u(&ifft);
      
// Split  
    Re[0] = ifft.OutBuf[0];
    Im[0] = 0;
    for (int16 i=1; i<(fftSize/2+1); i++) {
        Re[i] = ifft.OutBuf[i];            // split real
        Im[i] = ifft.OutBuf[fftSize-i]; // split imag
    }
    Im[fftSize/2] = 0;
// Combine and split imag/real
    for(int16 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
    for (int16 i=0; i<(fftSize/2); i++){
        Re[fftSize-1-i] = Re[i+1];
        Im[fftSize-1-i] = -Im[i+1];
    }
    // Modifed
    for (int16 i=0; i<fftSize/2+1; i++) fft2.InBuf[i] = Re[i]/fftSize;
    for (int16 i=fftSize/2+1; i<fftSize; i++) fft2.InBuf[i] = -Im[i]/fftSize;

Т.е. обходятся они одним FFT для реализации обратного преобразования....А ещё не очень я пока понял почему у него результат сдвигается на 3 отсчёта вправо.....

Но картина уже выясняется потихоньку...скоро буду перехдить к "полевым испытаниям" "боевых" алгоритмов

[sigmaN]

Код этих техасовских парней у меня не заработал.
Уж не знаю почему они его завалидировали и поставили там галочку.

Результат работы RFFT_f32u(); сверил с выводом KISS FFT - всё Ok.

Теперь, что касается обратного преобразования.
Как я размышляю(а вы поправьте где я ошибаюсь)
На выходе преобразования имеем N комплексных числа. вещественные части представляют собой амплитуды косинусов, мнимые части - синусов соответствующих частот. Хорошо, но нам нужно подать эту мешанину на вход преобразвателя, ожидающего вещественный input, т.е. с этим комплексным output'ом нужно что-то делать. Техасовские парни колбасят эти массивы, отражая, складывая и снова отражая и в конце концов у них что-то получается, но мне не понятен сам смысл.

В чём суть. Как из комплексного вывода получают "правильный" вещественный ввод, а потом снова из комплексного вывода получают "правильный" вещественный вывод??? Ведь чую, что не хватает мне нескольких слов, чтобы всё в голове уложилось.

Added: последние мысли.
1. БПФ вещественного сигнала сопряженно-симметрично. И моя функция выдаёт мне только половину(и в комплексном виде). Верно ли это? Думаю, что да.
2. Для того, чтобы подготовить сигнал к обратному БПФ - его вначале нужно дополнить этими недостающими комплексными сопряжениями первой половины. Вроде как тоже верно.
3. А вот теперь как из этого сделать ввод для вещественного БПФ? Упаковать нужно как-то эти комплексные числа в вещественный ввод
HEEEELP!

[bahurin]

Код этих техасовских парней у меня не заработал.
Уж не знаю почему они его завалидировали и поставили там галочку.

Не знаю как код а идея у них правильная. Давайте разбираться. Есть сигнал, и мы получили его комплексный спектр A+j*B. Надо взять IFFT и получим снова вещественный сигнал s = IFFT(A+j*B ) = IFFT(A)+j*IFFT(B ). Поскольку А и B вещественные, то IFFT можно заменить на FFT, только на выходе необходимо поделить результат на N (нормировка). Тогда получили s=FFT(A)+j*FFT(B ). А теперь самое главное. Поскольку А - разложение по косинусам а B - разложение по синусам, то FFT(A) = SA - чисто вещественно, а FFT(B ) = j*SB - чисто мнимое. Это значит, что можно заменить два FFT одним следующим образом: S = FFT(A+B ). Тогда S будет содержать комплексные значения и для расчета вещественного сигнала надо s = real(S)+imag(S). Вот матлабовский пример

Код

N = 16;
s  = rand(1,N); %случайный входной сигнал
S0=fft(s);      %FFT вещественного входного сигнала
A = real(S0);   %реалная часть
B = imag(S0);   %мнимая часть
%расчет ifft на осове вещесвенного fft
S=fft(A+B); % A+B - чисто вещественно, S - комплексное    
s0 = (real(S)+imag(S))/N; % результат IFFT
plot(1:N,s,'-o',1:N,s0)

в вашем случае необходимо учесть, что результат OutBuf содержит только половины спектра, которые надо шаманить зеркально отражая. Но если это все аккуратно сделать то должно работать. Вот матлабовский пример как надо вертеть массивы в вашем случае.

Код

clear all
j = sqrt(-1);
N = 16;
s  = rand(1,N); %случайный входной сигнал
S = fft(s);
r = real(S(1:N/2+1)); % это реальная часть половины спектра в OutBuf
i = fliplr(imag(S(2:N/2))); % это мнимая часть половины спектра в OutBuf

r = [r fliplr(r(2:N/2))];% ОТРАЖАЮ зеркально
i = [0 fliplr(i) 0 -i]; % ОТРАЖАЮ зеркально с минусом

s1 = fft(r+i)/N;

s1R = real(s1(1:N/2+1)); % это реальная часть половины спектра в OutBuf при fft(r+i)  
s1I = fliplr(imag(s1(2:N/2))); % это мнимая часть половины спектра в OutBuf при fft(r+i)

SR = [s1R fliplr(s1R(2:N/2))];% ОТРАЖАЮ зеркально
SI = [0 fliplr(s1I) 0 -s1I]; % ОТРАЖАЮ зеркально с минусом
s1 = SR+SI;

plot(1:N,s,'-o',1:N,s1)

Ищите ошибку в своем коде должно работать.

Сообщение отредактировал bahurin - Sep 18 2009, 04:12

[sigmaN]

Ну я был где-то рядом, но толком всё это разложить по полочкам не хватало ещё CPU power

Щас я уже с пониманием к вопросу подойду, думаю что всё заработает.


Огромное спасибо за консультэйшн


Кстати,а что у него там за проблема со сдвигом выходного сигнала на три отсчёта вправо?
I optimize you code for pre-IFFT computation and I need to change last "add" section for post-IFFT computation. And IFFT result is shifted for 3 samples rigth from source signal.
Хотя это уже детали, наверное сам разберусь когда толком всё отлажу и проверю.

[AndrewN]

Есть сигнал, и мы получили его комплексный спектр A+j*B. Надо взять IFFT и получим снова вещественный сигнал s = IFFT(A+j*B ) = IFFT(A)+j*IFFT(B ). Поскольку А и B вещественные, то IFFT можно заменить на FFT, только на выходе необходимо поделить результат на N (нормировка). Тогда получили s=FFT(A)+j*FFT(B ). А теперь самое главное. Поскольку А - разложение по косинусам а B - разложение по синусам, то FFT(A) = SA - чисто вещественно, а FFT(B ) = j*SB - чисто мнимое. Это значит, что можно заменить два FFT одним следующим образом: S = FFT(A+B ). Тогда S будет содержать комплексные значения и для расчета вещественного сигнала надо s = real(S)+imag(S).

Бинго! Беру своё предыдущее утвердение обратно. Ну до чего симметрия хороша :)

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

Спасибо.

[sigmaN]

Чё-то тут не сходится по-моему
Согласно доке - формат выходного буфера:
OutBuf[0] = real[0]
OutBuf[1] = real[1]
OutBuf[2] = real[2]
………
OutBuf[N/2] = real[N/2]
OutBuf[N/2+1] = imag[N/2-1]
………
OutBuf[N-3] = imag[3]
OutBuf[N-2] = imag[2]
OutBuf[N-1] = imag[1]
Всё это хорошо.
imag[0] всегда == 0 его выводить в буфер не имеет смысла
а вот как быть с real[N/2].
real[0]!=0 потому что cos(0)==1 и тут получается DC offset
A real[N/2] по логике должен быть == 0, потому что cos(pi/2)==0 но у меня он !=0 и мне не совсем понятно why?
Хотя скорее всего я что-то не догоняю, потому что вывод всё-таки соответствует KISS FFT и тут не может быть ошибки.

[bahurin]

Чё-то тут не сходится по-моему
Согласно доке - формат выходного буфера:
OutBuf[0] = real[0]
OutBuf[1] = real[1]
OutBuf[2] = real[2]
………
OutBuf[N/2] = real[N/2]
OutBuf[N/2+1] = imag[N/2-1]
………
OutBuf[N-3] = imag[3]
OutBuf[N-2] = imag[2]
OutBuf[N-1] = imag[1]
Всё это хорошо.
imag[0] всегда == 0 его выводить в буфер не имеет смысла
а вот как быть с real[N/2].
real[0]!=0 потому что cos(0)==1 и тут получается DC offset
A real[N/2] по логике должен быть == 0, потому что cos(pi/2)==0 но у меня он !=0 и мне не совсем понятно why?
Хотя скорее всего я что-то не догоняю, потому что вывод всё-таки соответствует KISS FFT и тут не может быть ошибки.

real[0] - постоянная составляющая может быть ллюбой, OutBuf[N/2] = real[N/2] также может быть любой. imag[0] =imag[N/2] = 0. В буфере real занимает адреса с 0 до N/2 всего N/2+1 значение, а imag c OutBuf[N/2+1] до OutBuf[N-1] всего N/2-1 значение при этом imag[0] =imag[N/2] = 0 не присутсвуют в буфере.

[sigmaN]

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]

Вопрос решен в таком виде:

Сложно. Можно сделать проще. Я не стал связываться с
кодом для С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
}

Сообщение отредактировал AndrewN - Sep 21 2009, 15:47

[sigmaN]

Даа, я естественно тоже пришел к выводу, что там возня лишняя присутствует.
И ресурсы как памяти так и проца расходуются непонятно куда.
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]

Воистину ассемблер и трпение - великая сила!

В итоге Speex echo canceller при праметрах
ECHO_FRAME_LEN = 32
filter tail size = 128
удалось вписать в 21MIPS! Это ещё осталось 2 кандидата на оптимизацию. Думаю после их оптимизации можно будет вписаться в 17-18, но это уже больше ради спортивного интереса....

Даже не верится, что до начала оптимизации с этими параметрами получалось около 70MIPS