|
c28х FPU library beta1, помогите разобраться с real FFT |
|
|
|
Sep 16 2009, 20:35
|
Местный
  
Группа: Участник
Сообщений: 336
Регистрация: 7-03-07
Из: Петербург
Пользователь №: 25 961

|
Цитата(sigmaN @ Sep 16 2009, 22:33)  a. Почему на выходе _REAL_ FFT видим imag, т.е. мнимую часть? b. Каков смысл FFT magnitude и phase? c. Как используя эти функции реализовать обратное преобразование? а. мнимые части появляются из-за комплексных матричных факторов в. mag = sqrt (re**2 + im**2), phase = atan(im/re) с. только специальной отдельной подпрограммой, у которой входные данные расположены в порядке выходных данных прямого преобразования.
|
|
|
|
|
Sep 16 2009, 20:53
|

I WANT TO BELIEVE
     
Группа: Свой
Сообщений: 2 617
Регистрация: 9-03-08
Пользователь №: 35 751

|
Ответ по пункту c. подразумевает редактирование исходников либы или же имеется ввиду подпрограмма-обёртка, формирующая вход? Мне кажется realFFT на то и real, что на входе у него только real, без мнимых частей. Тогда становится невозможным подать выход на вход и получить обратное преобразование. Следовательно, необходимо редактирование исходников либы. Что не очень-то и радует. Или я чего-то не понимаю. P.S. да нет, тут как ни крути, а придётся редактировать....  Но не верится, что люди выпустили оптимизированную либу без такой, казалось-бы очевидной и нужной вещи как обратное преобразование.
--------------------
The truth is out there...
|
|
|
|
|
Sep 17 2009, 03:16
|
Вечный ламер
     
Группа: Модераторы
Сообщений: 7 248
Регистрация: 18-03-05
Из: Томск
Пользователь №: 3 453

|
Цитата(sigmaN @ Sep 16 2009, 14:53)  Ответ по пункту c. подразумевает редактирование исходников либы или же имеется ввиду подпрограмма-обёртка, формирующая вход? Мне кажется realFFT на то и real, что на входе у него только real, без мнимых частей. Тогда становится невозможным подать выход на вход и получить обратное преобразование. Следовательно, необходимо редактирование исходников либы. Что не очень-то и радует. Или я чего-то не понимаю. P.S. да нет, тут как ни крути, а придётся редактировать....  Но не верится, что люди выпустили оптимизированную либу без такой, казалось-бы очевидной и нужной вещи как обратное преобразование. мой английский слаб, но мне кажеться что речь идет о реализации вычислений в формате с плавающей точкой при длине слова 32 бита. Насчет как посчитать IFFT через FFT смотрите http://www.dspsystem.narod.ru/content/fft/fft/fft.html
--------------------
|
|
|
|
|
Sep 17 2009, 05:34
|

Местный
  
Группа: Участник
Сообщений: 240
Регистрация: 20-09-08
Пользователь №: 40 347

|
Цитата(sigmaN @ Sep 16 2009, 23:33)  Почему на выходе _REAL_ FFT видим imag, т.е. мнимую часть? Каков смысл FFT magnitude и phase? Как используя эти функции реализовать обратное преобразование? а есть какое нить описание остальных буферов: float32 *InBuf;//Input buffer float32 *CosSinBuf; float32 *MagBuf; float32 *PhaseBuf; Для обратного преобразования надо чтобы в InBuf можно было положить как реальную так и мнимую части спектра.
|
|
|
|
|
Sep 17 2009, 05:41
|

I WANT TO BELIEVE
     
Группа: Свой
Сообщений: 2 617
Регистрация: 9-03-08
Пользователь №: 35 751

|
Собственно говоря у меня тот-же вопрос и в сопроводительной доке ответа на него нет. http://www.ti.com/litv/zip/sprc624 - либа вместе с докой Вечерком ковырну исходники и посмотрю. Ну или эксперимент организую..... Ну или накрайняк сапорт озадачить попробую...... дока отдельно:
--------------------
The truth is out there...
|
|
|
|
|
Sep 17 2009, 10:46
|

Местный
  
Группа: Участник
Сообщений: 240
Регистрация: 20-09-08
Пользователь №: 40 347

|
Цитата(sigmaN @ Sep 17 2009, 09:41)  Собственно говоря у меня тот-же вопрос и в сопроводительной доке ответа на него нет. 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
|
|
|
|
|
Sep 17 2009, 13:03
|

I WANT TO BELIEVE
     
Группа: Свой
Сообщений: 2 617
Регистрация: 9-03-08
Пользователь №: 35 751

|
Нашел пост на техасовском форуме Код 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 отсчёта вправо..... Но картина уже выясняется потихоньку...скоро буду перехдить к "полевым испытаниям" "боевых" алгоритмов
--------------------
The truth is out there...
|
|
|
|
|
Sep 17 2009, 23:09
|

I WANT TO BELIEVE
     
Группа: Свой
Сообщений: 2 617
Регистрация: 9-03-08
Пользователь №: 35 751

|
Код этих техасовских парней у меня не заработал. Уж не знаю почему они его завалидировали и поставили там галочку. Результат работы RFFT_f32u(); сверил с выводом KISS FFT - всё Ok. Теперь, что касается обратного преобразования. Как я размышляю(а вы поправьте где я ошибаюсь) На выходе преобразования имеем N комплексных числа. вещественные части представляют собой амплитуды косинусов, мнимые части - синусов соответствующих частот. Хорошо, но нам нужно подать эту мешанину на вход преобразвателя, ожидающего вещественный input, т.е. с этим комплексным output'ом нужно что-то делать. Техасовские парни колбасят эти массивы, отражая, складывая и снова отражая и в конце концов у них что-то получается, но мне не понятен сам смысл. В чём суть. Как из комплексного вывода получают "правильный" вещественный ввод, а потом снова из комплексного вывода получают "правильный" вещественный вывод??? Ведь чую, что не хватает мне нескольких слов, чтобы всё в голове уложилось. Added: последние мысли. 1. БПФ вещественного сигнала сопряженно-симметрично. И моя функция выдаёт мне только половину(и в комплексном виде). Верно ли это? Думаю, что да. 2. Для того, чтобы подготовить сигнал к обратному БПФ - его вначале нужно дополнить этими недостающими комплексными сопряжениями первой половины. Вроде как тоже верно. 3. А вот теперь как из этого сделать ввод для вещественного БПФ? Упаковать нужно как-то эти комплексные числа в вещественный ввод HEEEELP!
--------------------
The truth is out there...
|
|
|
|
|
Sep 18 2009, 04:11
|

Местный
  
Группа: Участник
Сообщений: 240
Регистрация: 20-09-08
Пользователь №: 40 347

|
Цитата(sigmaN @ Sep 18 2009, 03:09)  Код этих техасовских парней у меня не заработал. Уж не знаю почему они его завалидировали и поставили там галочку. Не знаю как код а идея у них правильная. Давайте разбираться. Есть сигнал, и мы получили его комплексный спектр 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
|
|
|
|
|
Sep 18 2009, 05:48
|
Местный
  
Группа: Участник
Сообщений: 336
Регистрация: 7-03-07
Из: Петербург
Пользователь №: 25 961

|
Цитата(bahurin @ Sep 18 2009, 08:11)  Есть сигнал, и мы получили его комплексный спектр 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). Бинго! Беру своё предыдущее утвердение обратно. Ну до чего симметрия хороша :) И всё равно, лучше запрограммировать специальную процедуру, что бы не возиться с суммами и перестановками. Спасибо.
|
|
|
|
2 чел. читают эту тему (гостей: 2, скрытых пользователей: 0)
Пользователей: 0
|
|
|