Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум программистов > Алгоритмы > Обратное быстрое преобразование Фурье...


Автор: Proger10 17.7.2010, 18:13
В википедии есть код http://ru.wikipedia.org/wiki/%D0%91%D1%8B%D1%81%D1%82%D1%80%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#.D0.9E.D0.B1.D1.80.D0.B0.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 на С++.

Насколько понимаю, они программировали эту формулу: http://upload.wikimedia.org/math/f/6/9/f69747fba84996b2d75638608c59a2bb.png

Постоянно запутываюсь в этом коде) может кто-то сможет отследить как из этого прямого сделать обратное преобразование Фурье. По теории, чтобы это сделать нужно экспоненту возвести в степень -1 (не ошибаюсь я?). Т.е. остаётся найти где они считают "exp", ну и сделать: "1/exp(...)". Только вот в этом поиске и заключается проблема.. smile) Может кто-нибудь поможет?

Вот рабочий код прямого преобразования Фурье:
Код

#define _USE_MATH_DEFINES
#include <math.h>
#include <stdio.h>
 
double *FFT( int *dIn, int nn )
{
    int i, j, n, m, mmax, istep;
    double tempr, tempi, wtemp, theta, wpr, wpi, wr, wi;
 
    int isign = -1;
    double *data = new double [ nn * 2 + 1 ];
 
    for( i = 0; i < nn; i++ )
    {
        data[ i * 2 ] = 0;
        data[ i * 2 + 1 ] = dIn[ i ];
    }
 
    n = nn << 1;
    j = 1;
    i = 1;
    while( i < n )
    {
        if( j > i )
        {
            tempr = data[ i ]; data[ i ] = data[ j ]; data[ j ] = tempr;
            tempr = data[ i + 1 ]; data[ i + 1 ] = data[ j + 1 ]; data[ j + 1 ] = tempr;
        }
        m = n >> 1;
        while( ( m >= 2 ) && ( j > m ) )
        {
            j = j - m;
            m = m >> 1;
        }
        j = j + m;
        i = i + 2;
    }
    mmax = 2;
    while( n > mmax )
    {
        istep = 2 * mmax;
        theta = 2.0 * M_PI / ( isign * mmax );
        wtemp = sin( 0.5 * theta );
        wpr = -2.0 * wtemp * wtemp;
        wpi = sin( theta );
        wr = 1.0;
        wi = 0.0;
        m = 1;
        while( m < mmax )
        {
            i = m;
            while( i < n )
            {
                j = i + mmax;
                tempr = wr * data[ j ] - wi * data[ j + 1 ];
                tempi = wr * data[ j + 1 ] + wi * data[ j ];
                data[ j ] = data[ i ] - tempr;
                data[ j + 1 ] = data[ i + 1 ] - tempi;
                data[ i ] = data[ i ] + tempr;
                data[ i + 1 ] = data[ i + 1 ] + tempi;
                i = i + istep;
            }
            wtemp = wr;
            wr = wtemp * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
            m = m + 2;
        }
        mmax = istep;
    }
    double *dOut = new double [ nn / 2 ];
 
    for( i = 0; i < ( nn / 2 ); i++ )
    {
        dOut[ i ] = sqrt( data[ i * 2 ] * data[ i * 2 ] + data[ i * 2 + 1 ] * data[ i * 2 + 1 ] );
    }
 
    delete []data; 
    return dOut;
}
 
int main()
{
    int *dsin = new int [ 1801 ], *pdsin = dsin;
 
    for( double x = 0.; x <= 10. * M_PI; x += M_PI / 180.0 )
    {
        *pdsin++ =  sin( x ) * 1000. + cos( 0.5 * x ) * 1000.;
    }
 
    double *dfourier = FFT( dsin, 1024 );
    delete []dsin;
 
    for( int i = 0; i < 512; i++ )
    {
        printf( "%d\t%lf\r\n", i, dfourier[ i ] );
    }
 
    delete []dfourier;
 
    return 0;
}


Автор: maxdiver 27.7.2010, 00:29
Да, реализация не самая лучшая (с точки зрения понятности; слишком она наоптимизированная, чтобы такую в вики пихать)...

В их обозначения переменных я не углублялся, но кажется, что здесь может быть единственный вариант:
Код

theta = 2.0 * M_PI / ( isign * mmax );

здесь всё умножить на -1 (или у isign, который может быть для этого и предназначен, изменить значение на противоположное).

И вообще, чтобы из прямого БПФ получить обратное, надо не только заменить экспоненту на противоположную, но и разделить потом каждый элемент результата на n.

Автор: Proger10 11.8.2010, 19:07
Опробовал, но не прокатило :(

В конце кода ФФТ, есть такая штука: dOut[ i ] = sqrt( .... ), случаем не может ли быть из-за неё?

У меня получаются вот какие графики: http://s001.radikal.ru/i196/1008/f5/65fd8b141fac.png
Тут видите, абсолютно какие-то другие числа получаются. Такое впечатление, что это кепстр исходного получается)))

Алгоритм такой:

ффтМассив = ФФТ( исходные_данные, 1024-окно );
ффтОбратноеМассив = рФФТ( ффтМассив, 1024 ); // рФФТ - обратное ффт, в котором заменён isign = -1 на +1.

И чего-то ерунда какая-то получается. Может у кого есть прямое и обратное FFT на C++?

Автор: Pavia 11.8.2010, 19:22
Proger10, 
// http://psi-logic.shadanakar.org/fft/fftf.htm
/// http://student.kuleuven.be/~m0216922/CG/fourier.html#introduction

Автор: Proger10 11.8.2010, 23:32
А покороче ничего нет? smile Чего там кода-то столько?)

Я вот нашёл некоторый код, короткий (он был комплексный, я им пользуюсь как не комплексным):
Код

/*
 This computes an in-place complex-to-complex FFT 
 x and y are the real and imaginary arrays of 2^m points.
 dir =  1 gives forward transform
 dir = -1 gives reverse transform 
 */
//short FFTc( short int dir, long m, double *x, double *y )
short FFTc( short int dir, long m, vector<double> &x, vector<double> &y )
{
    
    long n,i,i1,j,k,i2,l,l1,l2;
    double c1,c2,tx,ty,t1,t2,u1,u2,z;
    
    /* Calculate the number of points */
    n = 1;
    for (i=0;i<m;i++) 
        n *= 2;
    
    /* Do the bit reversal */
    i2 = n >> 1;
    j = 0;
    for (i=0;i<n-1;i++) {
        if (i < j) {
            tx = x[i];
            //ty = y[i];
            x[i] = x[j];
            //y[i] = y[j];
            x[j] = tx;
            //y[j] = ty;
        }
        k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
    
    /* Compute the FFT */
    c1 = -1.0; 
    c2 = 0.0;
    l2 = 1;
    for (l=0;l<m;l++) {
        l1 = l2;
        l2 <<= 1;
        u1 = 1.0; 
        u2 = 0.0;
        for (j=0;j<l1;j++) {
            for (i=j;i<n;i+=l2) {
                i1 = i + l1;
                t1 = u1 * x[i1] - u2 * y[i1];
                t1 = u1 * x[i1];
                t2 = u2 * x[i1];
                t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1; 
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            z =  u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
        c2 = sqrt((1.0 - c1) / 2.0);
        if (dir == 1) 
            c2 = -c2;
        c1 = sqrt((1.0 + c1) / 2.0);
    }
    
    /* Scaling for forward transform */
    if (dir == 1) {
        for (i=0;i<n;i++) {
            x[i] /= n;
            y[i] /= n;
        }
    }
    
    return(TRUE);
    
}



Он преобразование для 1024 точек считает секунд 20 на двуядерном 2.24GHz, что неприемлемо для меня (мне в рил-тайме обрабатывать надо сигнал). Вышеприведённый код считает меньше, чем за пол секунды.

Если я начну реализовывать самостоятельно, он, скорее всего, получится у меня таким же тормозом.. а почему же первый работает так быстро? Накрутили его нехило, зато, чёрт, очень работоспособный.. И то на половину. Обратное сделать в нём не получается.

Автор: Pavia 12.8.2010, 08:34
Proger10, 
2 код тормазит из-за vector. Используй указатель.
Хочешь скорость используй fftw
http://www.fftw.org/

Автор: Proger10 12.8.2010, 10:24
Вряд ли из-за vector он тормозит. Покуда прошлый код у меня тоже работает на vector. (я поменял с указателя на vector). Видимых отличий в скорости не оказалось, а работать удобнее, поэтому и оставил. 

Спасибо за совет. Попробую. Я как-то смотрел уже те мануалы.. первое впечатление - чёрт ногу сломит. Посмотрим ещё раз.

Автор: maxdiver 13.8.2010, 23:50
У меня на http://e-maxx.ru/algo/fft_multiply есть вполне короткий и адекватно быстрый код:

Код

const double PI = 3.1415926535897932384626433832795;
typedef complex<double> base;
 
void fft (vector<base> & a, bool invert) {
    int n = (int) a.size();
    if (n == 1)  return;
 
    vector<base> a0 (n/2),  a1 (n/2);
    for (int i=0, j=0; i<n; i+=2, ++j) {
        a0[j] = a[i];
        a1[j] = a[i+1];
    }
    fft (a0, invert);
    fft (a1, invert);
 
    double ang = 2*PI/n * (invert ? -1 : 1);
    base w (1),  wn (cos(ang), sin(ang));
    for (int i=0; i<n/2; ++i) {
        a[i] = a0[i] + w * a1[i];
        a[i+n/2] = a0[i] - w * a1[i];
        if (invert)
            a[i] /= 2,  a[i+n/2] /= 2;
        w *= wn;
    }
}


Для ста тысяч элементов отрабатывать должен где-то за секунду.

Автор: Proger10 16.8.2010, 20:40
maxdiver
А не в комплексном случае, а вещественном, инициироваться wn должна каким значением? Инициировал косинусом от ang - получилось что-то не очень походящее на спектр Фурье smile Нужно ли делать какую-то пост-обработку получившихся коэффициентов после этой функции? (я ещё прологорифмировал каждый получившийся коэффициент). И тем не менее, что-то не то получается..

Код у меня сейчас вот так модифицирован:
Код
#ifndef M_PI 
#define M_PI 3.14159265358979323846
#endif

void FFTd (vector<double> &a, bool invert) {

    int n = (int) a.size();
    if (n == 1)  return;
    
    vector<double> a0 (n/2),  a1 (n/2);
    for (int i=0, j=0; i<n; i+=2, ++j) {
        a0[j] = a[i];
        a1[j] = a[i+1];
    }
    FFTd (a0, invert);
    FFTd (a1, invert);
    
    double ang = 2 * M_PI / n * ( invert ? -1 : 1 );
    double w = 1.0;
    double wn = cos( ang );
    
    for (int i=0; i<n/2; ++i) {
        
        a[i] = a0[i] + w * a1[i];
        a[i+n/2] = a0[i] - w * a1[i];
        if (invert)
            a[i] /= 2,  a[i+n/2] /= 2;
        w *= wn;
        
    }
    
}


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

Добавлено
Квадраты суммы высчитал - стало уже получше, но всё равно очень уж много шума.. совсем непонять какие там частоты наиболее активны. Сигнал разный, а они практически одинаковые везде smile заметно, что чуток разнятся, но на глаз я бы не сказал ничего по такому спектру.. слишком сложно) а у первого кода - легко всё различимо..

Автор: Proger10 16.8.2010, 21:30
Чтобы было более понятно о какой такой точности двух кодов я говорю, вот картинка:
http://i064.radikal.ru/1008/c7/63134840417c.gif

левая колонка - верхний код этой темы, а правая - последний код (аналогичный сигнал). Вот о чём я и говорю.. блин, код хороший, короткий, но какой-то не очень точный smile как такой поаанлизируешь? smile

Автор: maxdiver 16.8.2010, 22:16
В БПФ числа всегда комплексные. Внутри функции ничего менять не надо. Входные параметры (а они вещественные - например, vector<double>) надо перевести в комплексные, для чего просто сделать вектор комплексных чисел и присвоить ему вектор вещественных (конвертация пройзойдет автоматически). Результат потом обратно надо перевести в действительные, для чего можно взять только действительную часть каждого комплексного числа (ну мнимая часть как бы = 0 в теории), или можно, как в коде в первом посте, взять модули всех чисел.

В коде, приведённом в первом посте, комплексные числа тоже есть, просто они хранятся не в виде отдельной структуры, а как два соседних дабла (опять же, ничего не скажешь, "удобно" для человека, пытающегося разобраться в БПФ и видящего перед собой этот код).

Автор: Proger10 16.8.2010, 23:08
Цитата

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

а где там модули берутся? не могу пока найти... Квадраты суммы квадратов, имеете ввиду?)

Опробовал с комплексными - это уже намного интереснее!! smile 

Автор: maxdiver 18.8.2010, 22:41
Вот где там брался модуль:
Код

dOut[ i ] = sqrt( data[ i * 2 ] * data[ i * 2 ] + data[ i * 2 + 1 ] * data[ i * 2 + 1 ] );


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

Автор: Proger10 19.8.2010, 17:43
А я тоже пока только БПФ использую smile

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)