Фурье с учетом предыдущих отсчетов, как сделать чтобы целочисленное работало? Опции

[alexPec]

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

Код

// Функция вычисления Фурье-преобразования последовательности {X1, X2,...,Xn} через заранее известное // Фурье последовательности {X0,X1,...,Xn-1}
// Re_prev, Im_prev   ---  действительная и мнимая части  Фурье предыдущего такта
// delta1 = (Re(Xn)-Re(X0)), delta2 = (Im(Xn)-Im(X0)), где Xn - следующий за последним, а X0 - начальный отсчеты предыдущего такта

void fourierShiftReduced( double *Re, double *Im, double delta ) {
        // вычитание спектра дельта-импульса
        for( int k = 0; k < fourierLength/2; k++ )
          Re[k] += delta1;
      Im[k] += delta2;

        // сдвиг спектра на 1/N
        double shift = 2 * PI / fourierLength;

        for( int k = 0; k < fourierLength/2; k++ )
        {
          phaseShift( Re[k], Im[k], shift * k );
        }
}

void phaseShift( double &re, double &im, double shift) {
        double a = cos( shift );
        double b = sin( shift );
        double x = a * re - b * im;
        im = b * re + a * im;
        re = x;
}

Проблема: когда используем арифметику с плавающей точкой, даже single, через 50х1024 входных отсчетов картина модуля спектра (fft_float). Когда все переводим на фиксированную точку, и именно синус, косинус умножаем на 65535 (получаем оцифровку 17 бит) и оставляем только целую часть, входной сигнал оцифровываем до 14 бит, то после 1024 отсчетов картина еще нормальная (fft_1), а через те же 50x1024 отсчета получаем это (fft_int). По ходу вычислений гармоники плывут какие-то вверх, какие-то вниз. С плавающей точкой движения никакого не наблюдается вообще. Похоже что по ходу в буфере накапливается ошибка округления - и результат нехороший. Пробовал синус, косинус оцифровывать до 19 бит - чуть медленее плывут.
Вопрос: можно ли в принципе реализовать таку схему с целочисленной арифметикой, или никак?

Сообщение отредактировал alexPec - Sep 1 2010, 05:58

Эскизы прикрепленных изображений

 

[alex_os]

...
Вопрос: можно ли в принципе реализовать таку схему с целочисленной арифметикой, или никак?

Можно попробовать сделать так, чтобы старые данные (и ошибка с ними вместе) "забывались".
т.е.

Код

      
          Re[k] += delta1;
      Im[k] += delta2;

заменить на

Код

      
          Re[k] = Re[k]  *аlpha + delta1;
      Im[k] = Im[k] *аlpha + delta2;

где alpha - число близкое к 1, но меньшее 1, например 0.9999.

Конечно это будет не совсем честное фурье, а что-то вроде фурье от сигнала на который наложено экспоненциальное окно
w = alpha^n , n=0...1023.