Привет!

Народ, не глянете, что я не так программирую в быстром преобразовании Фурье? Я пишу программный фильтр и пока занимаюсь отсеиванием наводки в 50 герц. Обычное вычисление дискретного преобразования Фурье дает красивую синусоиду, отлично вписывающуюся в сигнал, а вот БПФ дает какую-то прыгающую синусоиду, которая пляшет и по фазе и по амплитуде. В чем может быть дело?

Я использую рекусрсивный алгоритм БПФ с прореживанием по частоте, с 1024 точками. Эти 1024 точки промеряют сигнал таким образом, что на них попадает два полных колебания 50-герцовой наводки.

Я не совсем в нем разобрался, поэтому, где-то наверное ошибся. Помогите найти, вроде не должно быть сложно, кто разбирался!

Комплексной библиотекой я не пользуюсь, поэтому работаю отдельно с косинусоидальными и синусоидальными частями. Сначала я иницализирую таблицу ортогональных функций (синусов и косинусов), длинной в половину, то есть в 512 занчений:

for( i=0; i<SampleCount/2; ++i ) {
UnitCosine[i] = cos( -2 * M_PI * i / SampleCount );
UnitSine[i] = sin( -2 * M_PI * i / SampleCount );
}

Делаю это один раз.

Рекурсивная процедура устроена так:


void dif( int Length, int EvenStart, int UnitIndexStep ) {

// делим полученную длину пополам
int SubLength = Length / 2;

// получаем сдвиг индекса для второй половины
int OddStart = EvenStart + SubLength;

// если рекурсия дошла до длины в одну точку — прекращам
if( Length > 1 ) {

// пробегаем половину массива, а остальную половину берем при помощи сдвига
// для пользования единой таблицей ортогональных функций используем шаг u
for( i=0, u=0; i<SubLength; ++i, u+=UnitIndexStep ) {

// первое крыло "бабочки" — комплексная сумма значений из двух половин
EvenCosine = ValueCosine[EvenStart+i] + ValueCosine[OddStart+i];
EvenSine = ValueSine[EvenStart+i] + ValueSine[OddStart+i];

// второе крыло "бабочки" — комплексная разность значений из двух половин
OddPreCosine = ValueCosine[EvenStart+i] — ValueCosine[OddStart+i];
OddPreSine = ValueSine[EvenStart+i] — ValueSine[OddStart+i];

// комплексно умноженная на ортогональную функцию
OddCosine = OddPreCosine * UnitCosine[u] — OddPreSine * UnitSine[u];
OddSine = OddPreCosine * UnitSine[u] + OddPreSine * UnitCosine[u];

// кладем рассчитаные значения в тот же массив
ValueCosine[EvenStart+i] = EvenCosine;
ValueSine[EvenStart+i] = EvenSine;
ValueCosine[OddStart+i] = OddCosine;
ValueSine[OddStart+i] = OddSine;
}

// повторяем процедуру на следующем уровне рекурсии раздельно по половинам
dif( SubLength, EvenStart, UnitIndexStep * 2 );
dif( SubLength, OddStart, UnitIndexStep * 2 );
}



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

for( i=0; i<SampleCount; ++i ) {
ValueCosine[i] = Sample[i];
ValueSine[i] = 0;
}
dif( SampleCount, 0, 1 );


Результат получается расположенным в реверсивном коде. Для реверсивной индексации использую функцию:

unsigned reverse( unsigned value ) {
int res = 0;
for( int b=0; b<bitsize; ++b ) {
res |= (value & 0x1);
res <<= 1;
value >>= 1;
}
return res;
}

Проверку осуществляю путем расчета одной 50-герцовой гармоники по обратному преобразованию Фурье:

bin = AvgSampleCount; // = 2
for( i=0; i<SampleCount; ++i ) {
Noise[i] =
(ValueCosine[reverse(bin)]*sin(2*M_PI*i*bin/SampleCount) —
ValueSine[reverse(bin)]*cos(2*M_PI*i*bin/SampleCount));

}


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

Помогите, плз, кому не в лом!