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?
Как используя эти функции реализовать обратное преобразование?

[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

[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).

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

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

Спасибо.