Преобразование фурье

Еще раз о FFT

Может ли спектр быть больше самого сигнала? т.е. есть функция fft(double *In, int Size);

она делает преобразование фурье от буфера In и выдает его в тот же буфер.

Например
входной сигнал имеет амплитуду 32767. Может ли после преобразования фурье буфер In содержать значения больше 32767 (или меньше -32678)

Знаешь, сам когда-то возился с этим. Вообще не может, только в редких случаях, когда колебания накладываются и получается резонанс, но это редкий случай. Это теория, но на практике у меня выходило наоборот. Вот сорец, он  100\% работоспособный, переработанный алгоритм из книги Рабинера и Гоулда.
int math_inverse_bits( int value, int bits )
// инвертируем биты для преобразования Фурье
{

   int result = 0 ;
   int mask = 1 << (bits-1) ;

   for ( int i=0; i<bits; i++ )
   {
       
       if ( value & mask )
           result |= 1 << i ;

       mask = mask >> 1 ;

   }

   return ( result ) ;

}

BOOL math_fft( double* dbl_array, int* nSize )
{

   // определяем длину для преобразования фурье
   int tmp_size = *nSize ;
   for( int M=0; tmp_size>1; tmp_size/=2,M++ ) ;

   int fft_size = 1 << M ; // 1<<M == 2^M

   // подготавливаем массив
   std::complex<double>* fft_array =
                           new std::complex<double>[ fft_size ] ;
   ASSERT( fft_array ) ;

   // устанавливаем порядок для fft
   for ( int i=0; i<fft_size; i++ )
   {
       fft_array[ math_inverse_bits(i,M) ] =
           std::complex<double>(dbl_array[ i ],0.0) ;
   }

   double pi = 3.141592653589793 ;

   // M этапов
   for ( int l=0; l<M; l++ )
   {

       int le = 1 << (l+1) ;  // le - смещение между бабочками
       int le1 = le >> 1 ;    // le1 - размер бабочки
       std::complex<double> U ( 1.0, 0.0 ) ;
       std::complex<double> W ( cos( pi / le1  ), sin( pi / le1 ) ) ;

       for ( int j=0; j<le1; j++ )
       {
           for ( int i=j; i<fft_size; i+=le )
           {
               int ip = i + le1 ;
               std::complex<double> T = fft_array[ ip ] * U ;
               fft_array [ ip ] = fft_array[ i ] - T ;
               fft_array [ i ] = fft_array[ i ] + T ;
           }

           U *= W ;

       }

   }

   for ( i=0; i<fft_size / 2 ; i++ )
   {
       dbl_array[ i ] = std::abs( fft_array[ i ] ) ;
   }

   *nSize = fft_size / 2 ;

   delete[] fft_array ;   fft_array = NULL ;

   return ( TRUE ) ;

}

Исходник мне не нужен, мне нужно знать может ли фурье образ быть больше самого сигнала и в каких случаях это может быть?

и на сколько (или во сколько) фурье образ может быть больше сигнала

безусловно может... математическое доказательство - это например дельта-функция... что, в свою очередь, дает ответ на второй вопрос - не просто может а может быть больше в любое число раз (и в бесконечное тоже)...

А разве фурье от дельта функции не 1?
если берем дискретную послелдовательность то дельта функция будет выглядеть так.
0, 0, 0, 0, 0, .... 1, 0, 0, 0, 0,

дискретное фурье от нее разве будет больше 1?

Мне о реальных сигналах? Может ли так получится что из 16 битного сигнала получится после фурье 32 битный спектр ( к примеру)?

я наверное не совсем ясно выразился (продолжая тему матемитики). фурье от дельты - единица. но есть преобразование лапласа (взаимнооднозначное), которое в частности, используется в задачах на исследование поведения распределеия Гиббса в предельных случаях.
с помошью него показывается, это самое распределение и дельта являются образами/прообразами друг друга. Но распределение Гиббса конечно. Это и есть ответ на твой вопрос. :-[ пример реального сигнала - цуг синуса и его спектр.

Вот это мне понравилось:

Цитата

Вообще не может, только в редких случаях

Напоминает "В военное время значение синуса может достигать двух!" ))

Цитата vot, 04.04.02, 16:30:07

Вот это мне понравилось:
Напоминает "В военное время значение синуса может достигать двух!" ))

Ой-ой, больно ты умен :-) , ну подумаешь, потрещал я , вот.
На самом деле, ФФТ - загадка природы и доконца еще не изучена человечеством,
так что вы, ребята, зря паритесь, я сорец залил, вот и юзайте. А вообще я знаю человека, СУПЕРПРОФЕССИОНАЛА по преобразованию Фурье и вообще по обработке сигналов в цифре, могу порекомендовать в индивидуальном порядке, правда...он не с этой планеты, но думаю это вас не остановит.
поки

2Melan
Порекомендуй человека плиз!

Цитата Melan, 08.04.02, 19:41:28

могу порекомендовать в индивидуальном порядке, правда...он не с этой планеты, но думаю это вас не остановит.
поки

И мне плиз.

Если у кого есть вопросы то могу немного помочь про фурье, я вроде как достаточно уже с ним разобрался

Цитата Z@, 03.05.02, 23:48:21

Если у кого есть вопросы то могу немного помочь про фурье, я вроде как достаточно уже с ним разобрался

Здравствуй!
Мне нужны исходники преобразования Фурье на jave, где я мог бы их достать?

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

Цитата

Z@
13.03.02 19:07:22 Еще раз о FFT  
Может ли спектр быть больше самого сигнала? т.е. есть функция fft(double *In, int Size);
она делает преобразование фурье от буфера In и выдает его в тот же буфер.  
Например
входной сигнал имеет амплитуду 32767. Может ли после преобразования фурье буфер In содержать значения больше 32767 (или меньше -3267

Ответ может на пример чистая синусоида, при обратном приобразованиии делят на  Size

Блин, хочу  сделать преобразование Фурье.... Нифига не получается, вернее получается, но какая-то фигня. Числа охрененно большие, рисует полную парашу. Конечно спектр это напоминает, но совсем отдаленно.  Я вот вообще не  рублю: это преобразование  ничего "не знает ни о частоте дискретизиции сигнала, ни о ширине спектра, так что же означают выходные данные, амплитуду какой частоты они  показывают, и вообще, какая это нахрен амплитуда 10^10.  
Я вот думаю.... преобразование Фурье - это ведь по сути скалярное произведение векторов (т.е. сигнала на комплексные экспоненты) которое позволяет выявить меру присутствия комплексной экспоненты той или иной частоты в сигнале. Если я не прав, то поправьте меня скорее, пока не наделал глупостей и не погубил весь мир. Так вот, модет, чтобы не париться с выходными величинами нужно просто нужно привести весь сиггнал к амплитуде
от -1 до 1, потом сделать ффт и получить приоичный  результат.
 Кстати, кто-нибудь в кусре - неужели программы типа Winamp делают "честное" FFT для визуализации спектра или же там делаются какие-то упрощения и допущения, чтобы ускорить процесс и сделать картинку более красивой?

Цитата

Melan
07.10.02 20:43:44 Блин, хочу  сделать преобразование Фурье.... Нифига не получается, вернее получается, но какая-то фигня. Числа охрененно большие, рисует полную парашу. Конечно спектр это напоминает, но совсем отдаленно.  Я вот вообще не  рублю: это преобразование  ничего "не знает ни о частоте дискретизиции сигнала, ни о ширине спектра, так что же означают выходные данные, амплитуду какой частоты они  показывают, и вообще, какая это нахрен амплитуда 10^10.  
 

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

Однозначный ответ - может. И если в твоем сигнале какя-то частота будет сильно выражена(переносить большую мощность), то это будет выражено большим пиком в спектре. Связано это с тем что спектр характеризует мощности представленые различными частотами сигнала. И если сигнал достаточно гармоничный, то в его спектре обязательно будет пик, и очень большая вероятность, что амплитуда данного пика больше значения сигнала.
Ничего не мешает уменьшить на определенный коэфициэнт фурье образ, но при этом можно потерять мелкие детали спектра.

Цитата rodion, 08.10.02, 16:01:19

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

Ой, умник, только посмотрите на него, будешь ты мне еще впаривать что такое частота дискретизации, я вообще-то и  не спрашивал такую глупость. Я говорю, что FFT принимает в качекстве входного параметра  только массив отсчетов и  ничего не знает  о характеристиках сигнала.

Цитата Cubloid, 08.10.02, 16:32:15

И если сигнал достаточно гармоничный...

Как это... что сигнал может быть немного гармоничным, достаточно гармоничным и очень гармоничным? ;D

Сообщение отредактировано: .alex - 09.10.02, 08:04

Синус и косинус достаточно гармоничный сигнал. В нем вся мощность представлена в одной частоте.
Возьми разночастотный шум или сложи кучу сигналов синуса одинаковой мощности с рандомной частотой и получишь фурье образ в котором нельзя выделить основную гармонику.

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

2Cubloid: Это всё понятно, но выражение типа "достаточно гармоничный", ИМХО не корректное, можно сказать, что гармоническое колебание с большей амплитудой... А вообще для радиосигналов, видео и д.р. в спектре обычно присутствует основная гармоника (несущая частота), а весь полезный сигнал находится симметрично относительно неё.

В некоректности данного высказывания не могу не согласиться. Но вы же меня поняли, я думаю. Мне все таки позволительно делать ошибки, я всего лишь студент недоучка и двоешник к тому же.

2Cubloid: я тоже студент, но не двоешник.

Уважаемые, а как опредилить основную гармонику?

??? ??? ??? ??? ??? ???
Уважаемый школьники, студенты и просто умные люди, вы гуру в  DSP, так помогите же мне разобраться как это работает исходниками, инфой и прочей фигней. Хватит попусту базарить всякую ерунду, давайте займемся делом, раз уж тут всех нас интересует одна тема

function [Y,Flopy] = mybpf(X)
N=length(X);
k=0:N-1;
n=0:N-1;
n'*k;
EXPON=exp(-i*2*pi*n'*k/N);
Y=X*EXPON;
Программа для matlab реализующая просто дискретное преобразование Фурье.
А что еще вам интресно?

Цитата Melan, 14.10.02, 18:05:01

??? ??? ??? ??? ??? ???
Уважаемый школьники, студенты и просто умные люди, вы гуру в  DSP, так помогите же мне разобраться как это работает исходниками, инфой и прочей фигней. Хватит попусту базарить всякую ерунду, давайте займемся делом, раз уж тут всех нас интересует одна тема

dsp что это такое никогда не слышал

Цитата esperanto, 15.10.02, 13:52:00

dsp что это такое никогда не слышал

digital sound proccessing.....вроде

Digital Signal Processing вернее...

Цитата Alexei, 11.10.02, 18:31:51

Уважаемые, а как опредилить основную гармонику?

Не знаю что такое "основная гармоника", но если это то же самое что доминантная частота, то есть способ,который считает ее как положение "центра тяжести спектра".
Я когда-то считал.Получалось очень неплохо.Довольно точно.
Если нужно,пришлю исходник с-шный по почте.

Я несколько дней назад выводил дискретное ПФ из обычного.
И могу помочь в вашем вопросе.

Дискретное преобразование Фурье - это
    k-1    
g(u)=SUM (f(n)*cos(2*pi*n*u/k)-j*f(n)*sin(2*pi*n*u/k));
    n=0

k - число квантов и в сигнале и в спектре, оно одинаковое.

И максимум она достигает, когда подынтегральная функция имеет максимальный модуль, а её знак совпадает со знаком косинуса или синуса(смотря под что максимизируем, тут же и определим максимизованную частоту).
Косинус на его знак на максимальный модуль равно максимум умножить на косинус квадрат.
А интеграл по периоду от косинуса квадрат(с любой частотой) равно 4
В результате сумма получается равна 4*max*k

Этот результат сразу делится на k, но это уже на практике, а не по формуле

получается, что спектр превышает сигнал не более чем в четыре раза, да и этого значения она достигает в редких случаях.

Спасибо, наконец-то нашелся умный человек. А то тут  не форум, а базар  кокай-то был.

Извините, но в прошлом сообщении я допустил ошибку.
Дискретное преобразование, как я говорил, - это
      k-1
g(u)=SUM cos(2*pi*u*n/k)*f(n)-j*sin(2*pi*u*n/k)*f(n)
      n=k

сразу хочу привести формулу обратного преобразования фурье

              k-1
f(n)=1/k * SUM cos(2*pi*u*n/k)*g(u)+j*sin(2*pi*u*n/k)*g(u)
              u=k
(разница в замене n<->u g<->f, в изменении знака перед синусом и в ведении коэффициента)


Первая ошибка - в том, что чатоты равнозначны. Если частота равна половине от числа квантов, то значение в косинуса будут принимать не плавные значения, а лишь +- 1. Теперь, если сопоставить функцию принимающую значение +- u (где u-максимальное значение(примем за еденицу)) так, что в сумме умножая на синус всегда получается u. Тогда вся сумма будет равна k*u
Если затем спектр поделить на k, то получится u. Т.е. тогда размерность спектра будет такая же РАДУЙТЕСЬ!

Только тогда не сработает обратное ДПФ, приведённое выше, так как деление уже произведено.

Вторая ошибка была - в цифре 4. Она возникла из-за того, что я неправильно перешёл от интеграла к сумме - надо было ещё поделить на 2pi. Получилось бы вместо 4-х 0,63661977. Это - для первой гармоники.

Но так как я исправил и первую ошибку, и считал для k/2-й гармоники, то получилось ровно 1.

Я проверил в МАТЛАБе:
» a=0:999;                      -нумеруем кванты
» s=sign(cos(2*pi*500*a/1000)); -содержит +1 -1 +1 -1 +1 -1  и так 1000 квантов
» plot(a,s);                    -на графике сплошная стена - слишком частый меандр
» f=fft(s,1000);                -делаем ДПФ(в реализация - быстрое ПФ)
» plot(a,real(f));              -на графике все кванты примерно 0, кроме пятисотого, он - 1000
» hold                          -здесь же строим мнимую часть
Current plot held                  
» plot(a,imag(f),'-r');         -Красным цветом, чтоб отличить. Получается примерно ноль везде.
»        

Вся "энергия" меандра( sign(cos(t)) ) уровня U превращается в один импульс размером в число квантов умножить на U. Так что, разделив на число квантов, получим U

И последнее: если делаешь преобразование сам, то сохраняй битов сколько хочешь - всё равно информация теряется, так что можешь сохранять хоть больше, хоть, если припрёт, меньше.

Ещё совет:
Предельный спектр( соответствие частоте значения ПФ при сигнале меандр с этой же частотой
c(i+1)=SUM (cos(2*pi*i*(0:999)/1000)*sign(cos(2*pi*i*(0:999)/1000)))
)
который я считал для кванта 1 (=2/pi)и кванта 500 (=1)
кроме 0-го и k/2-квантов не превышает 0,65
На этом можно или сэкономить место, или сохранить больше битов.

У меня вопрос :
уважаемые кто нить занимался ФП под дельфи ? только не дискретным а полным (N) ?
если да поделитесь пожалуста алгоритмами и насколько медленнее БПФ это работает ?
есть ли смысл вообще делать ПФП или нет ?

Обычно преобразование - N^2 операций
БПФ - N*log2(N)
А какая разница где в Дельфи или где-то еще?
Алгоритм везде одинаков, синтаксис написания разный.

а чо есть бпф? или я что-либо пропустил в вашем обсуждении?

в смысле в чём заключается алгоритм, делающий пф за столь хорошее время?

Простите я наверное немного невнимательно прочитал вопрос от GEO.
Все написаное мною выше про число операций, касается ДПФ(Дискретного Преобразования Фурье), а БПФ есть ничто иное как алгоритм вычисления ДПФ используя меньшее число операций.
Смысл примерно такой БПФ: если полностью расписать 8-ми точечное ДПФ, то среди слагаемых будет ab+ac (3 операции), что можно вычислить a(b+c) (2 операции).
У меня сейчас нет под рукой конспекта лекций, но я думаю что это есть в Инете.
Вот собственно здесь:
http://www.dsp.sut.ru/book/lections/l4/l4_3.htm

2Cubloid

там слишком мало воды, непонятно, что какая переменная означает

Что я вам могу сказать на это:
Там все ясно по теории, но нет конечной реализации.
Только надо взять бумагу и ручку или пойти на Я.ру
Последний пункт я выбрал и на второй странице поиска нашел вот это.
http://home.tula.net/avm/fft/index.html
И там же ссылка на www.fftw.org
)) На последнем сайте представлена библиотека, которую можно использовать и большой набор ссылок на ресурсы по теме БПФ.
Что еще вам надо по этому вопросу? Все есть в Инете.

2Cubloid

ура! я понял! там не N log N умножений, а те же самые N^2, просто их в 2 раза меньше. ради этого и не стоило огород городить, в 2 раза я и сам могу уменьшить.....

8) Иди ты на хуй, придурок. Что значит N^2, только в два раза меньше. Ты долбоеб, не видишь разницы между NlogN и N^2 ????? Тогда вали отсюда, сука и больше здесь  не появляйся!!!!!! Нехуй людям попусту  голову морочить, пиздюк!

2melan

я так тоже могу :

2cubloid

я тут на досуге подумал и пришёл к следующим выводам:

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

2. там в натуре n log n умножений, однако там такое же количество сложений, даже больше. а вот у тебя в лекции приведён пример, который со всей присущей ему наглядностью показывает, как замечательно этот алгоритм делает 500 000 умножений вместо миллиона, что, конечно же, неправда.

3. всё это чушь.

как видно из рисунка, если мы умножим am или am+1 на тот же коэффициент k, что и an, то мы получим неверный результат.

Сообщение отредактировано: wormball - 25.03.04, 15:51

0. Наиболее авторитетный, ИМХО, источник - Рабинер и Голд. Выложил его на dsp-book.narod.ru
1. При выполнении Фурье по-разному решают вопрос нормирования. Математики чаще делят на нормирующий множитель и прямое, и обратное преобразования, а при программировании зачастую, экономя операции, нормируют (деля уже не на sqrt(n), а на n) только обратное преобразование. Но даже при математическом варианте - возможно драматическое возрастание элементов Фурье-образа по сравнению с исходными. В крайнем случае, согласно теореме Парсеваля, возможно возрастание в sqrt(n/2) раз, т.е. для того, чтобы пользоваться только 16-битной арифметикой при векторе из 512 элементов, исходный вектор должен быть не более чем 12-разрядный, а если брать "программистскую" нормировку - то 8-разрядный. Поэтому используется 32-битная или же плавающая (кстати, нынче это может быть и дешевле фиксированной, но с дополнительными операциями по нормировке на каждом шаге) арифметика.
2. Вычислять ДПФ "в лоб" стоит только при ОЧЕНЬ малых размерностях исходного вектора. В точной арифметике они дают одинаковый результат, а в реальной - БПФ точнее. Что до числа операций - "в лоб" надо O(n^2), а через БПФ не O(n^2/2) - понимающие поймут шутку - а O(n log n). При длине вектора порядка 100 и более - выигрыш в разы и десятки раз. (Лично я, искусно применив БПФ, смог ускорить вычисление корреляции в 300 раз...) Но, скажем, в JPEG считают прямо - поскольку там матрица 8х8 и вектора, соответственно, по 8 элементов.
3. Во многих случаях можно обойтись без Фурье. Скажем, если нужно установить наличие в сигнале строго определенных частот, поможет алгоритм Герцеля, если мощность в полосе - фильтрация, спектры могут вычисляться через авторегрессию или методом Писаренко и т.п.
4. В качестве доминирующей частоты сигнала можно брать максимальный пик в спектре, среднюю частоту, медианную частоту. Зависит от конкретной задачи.

Еще про БПФ.
Число умножений равно числу сложений и в БПФ и в "тупом" алгоритме. Так что экономия и на умножениях, и на сложениях.
А основная идея, скажем, алгоритма Кули-Тьюки проста.
ДПФ напрямую считается за квадратичное время. Но, имея два отрезка и их ДПФ, можно получить ДПФ объединения этих двух отрезков (в данном случае точки одного ставятся на четные, второго на нечетные позиции) за линейное время, поскольку попросту коэффициенты умножаются на множители поворота и складываются. Поэтому можно сделать вместо одной ДПФ за полное время два за четверть - и еще немного на линейноую по времени операцию слияния. Деля далее - обнаружим, что ДПФ делается от очень маленьких отрезков, время пренебрежимо, и основное время уходит на слияние. Но слияний нам нужно логарифм двоичный от Эн. Вот и получаем сложность.

2СанитарЖеня

непонял???????

Вот пример реализации на Паскале. Здесь все прозрачно.

procedure fft(var Ar:array of double; var Ai:array of double);
// На входе - действительная и мнимая части, на выходе они же, но для образа
var
Le,Le1,jt,it,ip,lt,kt:Integer;
Tr,Ti,Ur,Ui,Wr,Wi:Double;
// Предполагается, что точек 256, отсюда константы 255, 128 и 8.
// Изменения и обобщения очевидны
begin
    for it:=0 to 255 do
    begin
    kt:=128;
    jt:=0;
    Le:=1;
    Le1:=it;
    for lt:=1 to 8 do
    begin
     if kt<=Le1 then begin
     jt:=jt+Le;
     Le1:=Le1-kt;
     end;
     Le:=Le*2;
     kt:=kt div 2;
    end;
// Это предобработка - перестановка элементов в порядке инвертированных битов №
    if it<jt then begin
       Tr:=Ar[jt];Ti:=Ai[jt];
       Ar[jt]:=Ar[it];Ai[jt]:=Ai[it];
       Ar[it]:=Tr;Ai[it]:=Ti;
       end;
    end;

// Основной цикл
    for lt:=1 to 8 do
    begin
     Le:=1;
     for jt:=1 to lt do Le:=Le*2;
     Le1:=Le div 2;
     Ur:=1.0;
     Ui:=0.0;
// Значение синуса и косинуса вычисляется по формуле С. и К. суммы углов
     Wr:=cos(Pi/Le1);
     Wi:=sin(Pi/Le1);
// Проход по элементам, с постоянно уменьшающимся шагом
     for jt:=0 to Le1-1 do
     begin
      it:=jt;
      while it<256-Le1 do
      begin
       ip:=it+Le1;
// Бабочка. Оптимизирована по скорости
       Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
       Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
       Ar[ip]:=Ar[it]-Tr;
       Ai[ip]:=Ai[it]-Ti;
       Ar[it]:=Ar[it]+Tr;
       Ai[it]:=Ai[it]+Ti;
       it:=it+Le;
      end;
// Вычисление новых значений синуса и косинуса
     Tr:=Ur*Wr-Ui*Wi;
     Ti:=Ur*Wi+Ui*Wr;
     Ur:=Tr;
     Ui:=Ti;
     end;
    end;
end;

2СанитарЖеня

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

Эээ... А понять ответ лень не помешает?
Поворачивающие множители - связаны с известным свойством Фурье, а именно, что при сдвиге (а объединяемые отрезки сдвинуты относительно друг друга) Фурье-образ претерпевает поворот (по фазе сдвигается, попросту).
А откуда новые частоты - оттуда, где и были. Но, для короткого отрезка, они не могли быть выделены явно, а искажали образы отрезков. По разному искажали. Сумма образов отрезков (поправленных корректирующими множителями) дает уточненную величину тех частот, которые были явно видимы в образах, а разность - дает это искажение, по которому и восстанавливается образ в промежуточных точках.

2СанитарЖеня

сколько ни силился, всё равно не понял. твоё объяснение похоже на древний алхимический текст. ты бы чтоль формулу написал, как выходят коэффициенты при нечётных гармониках?

Все к Рабинеру и Голду (или, кто не желает - к Рабинеру и Гоулду). dsp-book.narod.ru
Там подробно изложено. Можно и Гольденберга, Матюшкина и Поляка.
Вкратце - вычитаем две оценки четных компонент, и их разность дает нечетную компоненту. Ну и те же поворачивающие множители.

2СанитарЖеня

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

Ну вот же они, формулы. Вот она, бабочка. Что именно неясно?        
       Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
       Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
       Ar[ip]:=Ar[it]-Tr;
       Ai[ip]:=Ai[it]-Ti;
       Ar[it]:=Ar[it]+Tr;
       Ai[it]:=Ai[it]+Ti;

Цитата СанитарЖеня, 31.03.03, 11:20:58

Ну вот же они, формулы. Вот она, бабочка. Что именно неясно?        
       Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
       Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
       Ar[ip]:=Ar[it]-Tr;
       Ai[ip]:=Ai[it]-Ti;
       Ar[it]:=Ar[it]+Tr;
       Ai[it]:=Ai[it]+Ti;

я ща заплачу (заплАчу :'(). чтотакое tr, ar, ur, далее везде???????????????????

Господи... Ну как можно объяснить проще? Ну нет в ДСП объяснений на уровне "Ты мне пальцем ткни!" - считать надо, и думать...
Кристально ясный же алгоритм...
Ну почитайте же Рабинера и Голда (АКА Гоулд).
dsp-book.narod.ru
Или Гольденберга, Матюшкина и Поляка
dspbook.km.ru.

Что есть Ar и Ai - в комментарии вначале программы написано, Ur и Ui - действительная и мнимая часть поворачивающего множителя, Tr и Ti - промежуточные переменные, необходимость их вызвана тем, что операции выполняются последовательно.

Очень нужна помощь в организации Фурье-фильтра в Matlab.
Смысл задачи:имеется 2 массива данных, один - частоты, второй - поглощение по ним (фурье-спектроскопия веществ). Необходимо, выделив одну частоту получить, на выходе поглощения по этой частоте.
Заранее благодарен, это очень сторочно!
P.S. Сам я в Матлабе ламер, если есть возможность, то опишите подробно.

to Ciceron
Нужно делать цифровой фильтр (полосовой), по-моему в MatLab есть примеры.

to All
Есть несколько приемов получения FFT (sic!) для призвольного числа точек (<>2^N).
Самый простый -- дополнение 0-ми до 2^N (перерасход памяти и снижение скорости). Посложнее -- разложение на простые множители для, которые и вычисляются аналитически, в случае =2^N этим простым множителем является 2 -->"бабочка", аналогично можно расписать преобразование для 3, 5, 7 и  т.д. Самый   худший случай когда число точек простое и необходимо выполнять DFT.

Чуваки, дайте лучше сцылку на оптимизированное 16бит в 16ит преобразование без использования FPU

Разобрался

Сообщение отредактировано: greek - 03.01.06, 19:00

У меня возник вопрос по-поводу функции окна. при расчете ее значений используется ограничение |nTd|<NTd/2, как и при расчете ДПФ, для неограниченого во времени сигнала, но там оно не играет роли поскольку не влияет на значение функции f(nTd), тут же оно оспользуеться, насколько я понимаю для того, что бы придать функции вид треугольника симметричного относительно 0. то есть для рачета нужно пробегать значения в диапазоне [-(N-1/2)..(N-1)/2], но в этом случае количество точек на 1 меньше реального, и что делать, когда N-нечетно.Насколько я себе представляю, рассчитывать функцию окна проще всего для первых -(N-1)/2 а потом записывать его в формулу симметрически с конца. И еще вопрос, в ДПФ нужно пермножать значения функции окна и значения функции на отсчетах, а в ОДПФ-их спектры. Если я не прав, то обьясните, как ей грамотно пользоваться

Здравствуйте уважаемые форумчане.

У меня возникла маленькая проблемка. Необходимо написать программу по расчету БПФ. Брать уже имеющиеся исходники я не хочу, поскольку хочу разобраться во всем сам.

Если я все правильно понимаю ДПФ выражается через формулу:

Xj = SUM (Xk*e^(-i*(2Pi/N)*kj)) (пределы суммирования по k от 0 до N-1)

e^(-i*(2Pi/N)*kj) = Cos((2Pi/N)*kj)+iSin((2Pi/N)*kj) (1)
где:
Cos((2Pi/N)*kj) = Re
iSin((2Pi/N)*kj) = Im

Так вот, вычисление действительной (Re) части выражения (1) не составляет труда, а вот с вычислением мнимой части (Im) у меня большие затруднения. А именно, я не могу понять как вычислить iSin((2Pi/N)*kj), как избавиться от этой проклятущей i и необходимо ли вообще вычислять мнимую часть выражения (1).

Люди добрые! Помогите, объясните дауну, что к чему!!!

Сообщение отредактировано: pmax - 12.01.06, 09:39

люди, помогите ламеру! у меня имеется ф-ция от времени x(t)=sin(20*t)+3*cos(30*t)

Необходимо в Матлабе:

1.исользовать бпф и построить частотный спектр
2.применить оконную ф-цию(любую)+бпф и посторить спектр!
3.применить фильтр + бпф и построить спектр
4.достроить исходный сигнал нулями и построить его спектр

Выручайте..!!

Заранее огромное спасибо

4elovek, ты хочешь чтоб тебе готовый код Matlab выложили?

я в принципе со всем разобрался! вопрос только в применении окна! Допустим Yn - значения БПФ! W - оконная ф-ция (пусть это будет прямоугольное окно boxcar)

Yn=fft(xn); считаем бпф для исходного дискретного сигнала xn.
W=boxcar(n); оконная ф-ция в n кол-ве отсчетов.

Вопрос:

окно применяется как G=Yn*W?

Если да, то при умножении получается число, т.к. Yn - вектор строка, а W - столбец! как быть? что делать?

!!!!!!И если кому не трудо приведите плиз элементарный пример применения любого окна в MAtlab.!!!!!

СанитарЖеня, а каким образом в твоем примере реализовать обратное преобразования.. а то чет туплю никак не соображу

У него нет обратного преобразования.

Вообще, обратное знаком перед мнимой составляющей (синусной) поворачивающего множителя.
Плюс результат умножается на другой множитель