Код
// Функция вычисления Фурье-преобразования последовательности {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;
}
// 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 бит - чуть медленее плывут.
Вопрос: можно ли в принципе реализовать таку схему с целочисленной арифметикой, или никак?
