Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > C/C++: Для новичков > Разложение в ряд Фурье

Автор: unitto 25.4.2010, 21:23

Доброго времени суток. Полтора часа насилования гугла не помогло найти теоретические сведенья для решения задачи:

Виконати чисельний спектральний аналіз і синтез періодичної функції f(t) (за n першими гармоніками розкладання в ряд Фур’є). Результат спектрального синтезу функції та вихідну функцію зобразити графічно.

В частности меня интересует как в С++ можно заданую функцию разложить на гармоники ряда Фурье, если можно с сылками на источники.
Заранее благодарен.

Автор: GoldFinch 25.4.2010, 21:25

гугли FFT

Автор: ToTToRo 26.4.2010, 05:40

google в помощь.

Автор: xvr 26.4.2010, 12:36

Плохо искал - http://ru.wikipedia.org/wiki/%D0%94%D0%B8%D1%81%D0%BA%D1%80%D0%B5%D1%82%D0%BD%D0%BE%D0%B5_%D0%BF%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A4%D1%83%D1%80%D1%8C%D0%B5

Автор: unitto 26.4.2010, 22:05

Цитата(xvr @ 26.4.2010, 12:36)

есть предложения по работе с формулами на этой странице в с++ ?

Автор: bilbobagginz 26.4.2010, 22:18

unitto,

если ты не понимаешь математику этого действия, поищи готовые реализации DFT (и в частности FFT) в сети с открытым кодом, и попробуй понять как пример.

надо понять общую идею перехода из измерения времени в измерение частоты, http://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A4%D1%83%D1%80%D1%8C%D0%B5
надо понять бесконечность ряда произвольной и непрерывной функции.

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

Автор: xvr 26.4.2010, 22:30

Цитата(unitto @ 26.4.2010, 22:05

)

есть предложения по работе с формулами на этой странице в с++ ?

Есть.
#include <complex> далее std::complex<double> для всех величин с 'этой страницы'
Или не понятно как перемножать и складывать ряды чисел?

Автор: unitto 27.4.2010, 00:00

Цитата(bilbobagginz @ 26.4.2010, 22:18)

ок, пасибо

Добавлено через 1 минуту и 42 секунды

Цитата(xvr @ 26.4.2010, 22:30)

Цитата(unitto @ 26.4.2010, 22:05

)

есть предложения по работе с формулами на этой странице в с++ ?

не работал с комплексными числами в с++
и уверен что в этой работе можно обойтись без них...

Автор: xvr 27.4.2010, 09:17

Цитата(unitto @ 27.4.2010, 00:00

)

и уверен что в этой работе можно обойтись без них...

Преобразование Фурье делается в комплексных числах. Но если хочется самому написать всю арифметику комплексных чисел - флаг в руки smile

Код


#include <complex>
#include <vector>

#define _USE_MATH_DEFINES
#include <cmath>

using namespace std;

vector<complex> DFT(vector<complex>& x)
{
 const size_t N=inp.size();
 vector<complex> rv;
 for(size_t k=0;k<N;++k)
  {
   complex acc(0,0);
   for(size_t n=0;n<N;++n)
    acc+=x[n]*exp(complex(0,-2*M_PI*k*n/N));
   rv.push_back(acc);
  }
 return rv;
}

Автор: Alexeis 27.4.2010, 09:42

Чисто случайно, мне удалось добыть фрагмент ультра секретного кода целочисленного FFT под кодовым названием "бабочка" smile

К сожалению, в коде отсутствуют некоторые фундаментальные вещи, которые делают его невозможным для использования. Секретная константа L скорее всего означает двоичный логарифим по размеру массива пример log2(256)=8 . CosQuot = size >> 2; (size - число точек фурье) . SIN_COS_AMPLITUDE - амплитуда синуса в таблице синусов. Таблица синусов tbl->arraySinCos[] представляет собой целочисленный массив размера size, содержащая полный период синуса значения которого умножены на SIN_COS_AMPLITUDE.

Код


int CFourier::GetCos(int index)
{
    index += CosQuot;
    if(index > size)
        index -= size;
    return tbl->arraySinCos[index];
}

void CFourier::UltraNFFT(int *re, int *im)
{
  int i, j, k;
  int ind, md2, mdd, mdd1, i1, i2;
  int a, b, c, d, batr, bati;

  ind = 0;        //перестановка отсчетов в битреверсном порядке как для re[] так и для im[]
  nm[0] = 0;

  for (i = 0; i < L; i++) 
  {
    for (j = 0; j < (1 << i); j++)
    {
      if (ind > nm[ind])
      {
        a           = re[ind];
        re[ind]     = re[nm[ind]]; 
        re[nm[ind]] = a;
        a           = im[ind];
        im[ind]     = im[nm[ind]];
        im[nm[ind]] = a;
      }
      ind++;
      nm[ind] = nm[j] + (1 << (L-1-i)); 
    }
  }

  md2 = size >> 1;
  for (k = 0; k < L; k++ )  //цикл по стадиям  (10 стадий)
  {                   
    mdd  = md2 >> k;
    mdd1 = md2 >>(L - k - 1);
    for (i = 0; i < mdd; i++ )  
    {                        //цикл по группам в каждой стадии  (512x1, 256x2, 128x4, 64x8, 16x32,…, 1x512)
      i1 = 2 * mdd1 * i;
      i2 = 2 * mdd1 * i + mdd1;
      for (j = 0; j < mdd1; j++)   
      {                      //цикл по парам внутри группы (выполнение бабочки для каждой пары в группе)
        a = re[j + i1];
        b = im[j + i1];
        c = re[j + i2];
        d = im[j + i2];
        ind = j * mdd;
        batr = (c * GetCos(ind) + 
                d * tbl->arraySinCos[ind]) / SIN_COS_AMPLITUDE;
        bati = (d * GetCos(ind) - 
                c * tbl->arraySinCos[ind]) / SIN_COS_AMPLITUDE;

        re[j + i1] = a + batr;
        im[j + i1] = b + bati;
        re[j + i2] = a - batr;
        im[j + i2] = b - bati;
      }
    }
  }
}

Добавлено через 2 минуты и 53 секунды
Забыл упомянуть что массив im должен быть заполнен нулями перед вызовом функции. В нем будет возвращена комплексная часть спектра.

Автор: unitto 27.4.2010, 09:46

правильно ли я понял:
$user posted image$
наша заданая функция f(x) есть сума гармоник fn(x)
тогда гармонику можно будет "изьять" "лобовым" делением на e^inx ?

Автор: Alexeis 27.4.2010, 10:03

Цитата(unitto @ 27.4.2010, 08:46

)

правильно ли я понял:

Это аналитическая формула. Формула дискретного фурье выглядит так.

$user posted image$

Автор: unitto 30.4.2010, 08:33

я так понимаю делать надо с помощью простого фурье с плавающей точкой?

Автор: Alexeis 30.4.2010, 09:15

Цитата(unitto @ 30.4.2010, 07:33

)

я так понимаю делать надо с помощью простого фурье с плавающей точкой?

Это проще. Не нужно следить за выходом из диапазона. Собственно xvr уже представил алгоритм в наиболее простой и наглядной его форме. Чтобы получить амплитуды гармоник нужно взять модуль от каждого полученного комплексного числа. Тут 2 варианта. Либо использовать 2 массива чисел с плавающей точкой, либо 1 массив комплексных. С точки зрения математики нагляднее всего работать именно с комплексными числами, тем более что для вас производительность не играет совершенно ни какой роли.

Автор: xvr 30.4.2010, 10:36

Цитата(unitto @ 30.4.2010, 08:33)

я так понимаю делать надо с помощью простого фурье с плавающей точкой?

Видите в формуле в показателе степени у 'e' стоит такая буква 'i'? Это мнимая единица. Так что с 'простой плавающей точкой' - облом. Не простая она, а комплексная smile

Если хотите простоты - берите мой код, если скорости - код от Alexis, в нем комплексная арифметика расписанна явно (в виде операций отдельно над действительными и мнимыми частями) и измененн алгоритм (вычисляется не Дискретное Преобразование Фурье а Быстрое Преобразование Фурье, результат тот же, но быстрее и запутаннее)

Добавлено через 8 минут и 47 секунд
Да, если вам просто нужен спектр и не обязательно Фурье - делайте преобразование Хартли, оно без комплексных чисел