Народ, не глянете, что я не так программирую в быстром преобразовании Фурье? Я пишу программный фильтр и пока занимаюсь отсеиванием наводки в 50 герц. Обычное вычисление дискретного преобразования Фурье дает красивую синусоиду, отлично вписывающуюся в сигнал, а вот БПФ дает какую-то прыгающую синусоиду, которая пляшет и по фазе и по амплитуде. В чем может быть дело?
Я использую рекусрсивный алгоритм БПФ с прореживанием по частоте, с 1024 точками. Эти 1024 точки промеряют сигнал таким образом, что на них попадает два полных колебания 50-герцовой наводки.
Я не совсем в нем разобрался, поэтому, где-то наверное ошибся. Помогите найти, вроде не должно быть сложно, кто разбирался!
Комплексной библиотекой я не пользуюсь, поэтому работаю отдельно с косинусоидальными и синусоидальными частями. Сначала я иницализирую таблицу ортогональных функций (синусов и косинусов), длинной в половину, то есть в 512 занчений:
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];
Понимаю, что где-то наверное не учел нормировку, но все равно непонятно, почему гармоника получается прыгающей. Пытался подобрать в последней части знаки плюс/минус, функции синусы/косинусы и т.п., но картина качественно не меняется, наверное ошибка в рекурсии или в расчете ортогональных функций.
Помогите, плз, кому не в лом!
Re: Что я делаю не так, реализуя быстрое преобразование Фурь
Здравствуйте, dims12, Вы писали:
D>Результат получается расположенным в реверсивном коде. Для реверсивной индексации использую функцию:
Возможно я и ошибаюсь, поскольку писал FFT давно, да и то итеративный вариант, а не рекурсивный, но
вроде бы не результат будет лежать в обратном битовом порядке, а исходниые данные
должны быть так размещены. Просто слабо верится, что операции рекурсивного FFT и битового
обращения коммутируют...
D>Проверку осуществляю путем расчета одной 50-герцовой гармоники по обратному преобразованию Фурье:
D>bin = AvgSampleCount; // = 2 D>for( i=0; i<SampleCount; ++i ) { D>Noise[i] = D>(ValueCosine[reverse(bin)]*sin(2*M_PI*i*bin/SampleCount) - D>ValueSine[reverse(bin)]*cos(2*M_PI*i*bin/SampleCount)); D>}
Вот это и вызывает сомнения -- обращения вида ValueCosine[reverse(bin)], ValueSine[reverse(bin)].
Re: Что я делаю не так, реализуя быстрое преобразование Фурь
/*============================================================================
fourierd.c - Don Cross <dcross@intersrv.com>
http://www.intersrv.com/~dcross/fft.html
Contains definitions for doing Fourier transforms
and inverse Fourier transforms.
This module performs operations on arrays of 'double'.
Revision history:
1998 September 19 [Don Cross]
Updated coding standards.
Improved efficiency of trig calculations.
============================================================================*/#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <fourier.h>
#define DDC_PI (3.14159265358979323846)
#define CHECKPOINTER(p) CheckPointer(p,#p)
static void CheckPointer ( void *p, char *name )
{
if ( p == NULL )
{
fprintf ( stderr, "Error in fft_double(): %s == NULL\n", name );
exit(1);
}
}
void fft_double (
unsigned NumSamples,
int InverseTransform,
double *RealIn,
double *ImagIn,
double *RealOut,
double *ImagOut )
{
unsigned NumBits; /* Number of bits needed to store indices */unsigned i, j, k, n;
unsigned BlockSize, BlockEnd;
double angle_numerator = 2.0 * DDC_PI;
double tr, ti; /* temp real, temp imaginary */if ( !IsPowerOfTwo(NumSamples) )
{
fprintf (
stderr,
"Error in fft(): NumSamples=%u is not power of two\n",
NumSamples );
exit(1);
}
if ( InverseTransform )
angle_numerator = -angle_numerator;
CHECKPOINTER ( RealIn );
CHECKPOINTER ( RealOut );
CHECKPOINTER ( ImagOut );
NumBits = NumberOfBitsNeeded ( NumSamples );
/*
** Do simultaneous data copy and bit-reversal ordering into outputs...
*/for ( i=0; i < NumSamples; i++ )
{
j = ReverseBits ( i, NumBits );
RealOut[j] = RealIn[i];
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
}
/*
** Do the FFT itself...
*/
BlockEnd = 1;
for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
{
double delta_angle = angle_numerator / (double)BlockSize;
double sm2 = sin ( -2 * delta_angle );
double sm1 = sin ( -delta_angle );
double cm2 = cos ( -2 * delta_angle );
double cm1 = cos ( -delta_angle );
double w = 2 * cm1;
double ar[3], ai[3];
double temp;
for ( i=0; i < NumSamples; i += BlockSize )
{
ar[2] = cm2;
ar[1] = cm1;
ai[2] = sm2;
ai[1] = sm1;
for ( j=i, n=0; n < BlockEnd; j++, n++ )
{
ar[0] = w*ar[1] - ar[2];
ar[2] = ar[1];
ar[1] = ar[0];
ai[0] = w*ai[1] - ai[2];
ai[2] = ai[1];
ai[1] = ai[0];
k = j + BlockEnd;
tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];
ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];
RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;
RealOut[j] += tr;
ImagOut[j] += ti;
}
}
BlockEnd = BlockSize;
}
/*
** Need to normalize if inverse transform...
*/if ( InverseTransform )
{
double denom = (double)NumSamples;
for ( i=0; i < NumSamples; i++ )
{
RealOut[i] /= denom;
ImagOut[i] /= denom;
}
}
}
/*--- end of file fourierd.c ---*/
Прочее:
/*============================================================================
fftmisc.c - Don Cross <dcross@intersrv.com>
http://www.intersrv.com/~dcross/fft.html
Helper routines for Fast Fourier Transform implementation.
Contains common code for fft_float() and fft_double().
See also:
fourierf.c
fourierd.c
..\include\fourier.h
Revision history:
1998 September 19 [Don Cross]
Improved the efficiency of IsPowerOfTwo().
Updated coding standards.
============================================================================*/#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <fourier.h>
#define TRUE 1
#define FALSE 0
#define BITS_PER_WORD (sizeof(unsigned) * 8)
int IsPowerOfTwo ( unsigned x )
{
if ( x < 2 )
return FALSE;
if ( x & (x-1) ) // Thanks to 'byang' for this cute trick!return FALSE;
return TRUE;
}
unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo )
{
unsigned i;
if ( PowerOfTwo < 2 )
{
fprintf (
stderr,
">>> Error in fftmisc.c: argument %d to NumberOfBitsNeeded is too small.\n",
PowerOfTwo );
exit(1);
}
for ( i=0; ; i++ )
{
if ( PowerOfTwo & (1 << i) )
return i;
}
}
unsigned ReverseBits ( unsigned index, unsigned NumBits )
{
unsigned i, rev;
for ( i=rev=0; i < NumBits; i++ )
{
rev = (rev << 1) | (index & 1);
index >>= 1;
}
return rev;
}
double Index_to_frequency ( unsigned NumSamples, unsigned Index )
{
if ( Index >= NumSamples )
return 0.0;
else if ( Index <= NumSamples/2 )
return (double)Index / (double)NumSamples;
return -(double)(NumSamples-Index) / (double)NumSamples;
}
/*--- end of file fftmisc.c---*/
fourier.h
/*============================================================================
fourier.h - Don Cross <dcross@intersrv.com>
http://www.intersrv.com/~dcross/fft.html
Contains definitions for doing Fourier transforms
and inverse Fourier transforms.
============================================================================*/#ifdef __cplusplus
extern"C" {
#endif/*
** fft() computes the Fourier transform or inverse transform
** of the complex inputs to produce the complex outputs.
** The number of samples must be a power of two to do the
** recursive decomposition of the FFT algorithm.
** See Chapter 12 of "Numerical Recipes in FORTRAN" by
** Press, Teukolsky, Vetterling, and Flannery,
** Cambridge University Press.
**
** Notes: If you pass ImaginaryIn = NULL, this function will "pretend"
** that it is an array of all zeroes. This is convenient for
** transforming digital samples of real number data without
** wasting memory.
*/void fft_double (
unsigned NumSamples, /* must be a power of 2 */int InverseTransform, /* 0=forward FFT, 1=inverse FFT */double *RealIn, /* array of input's real samples */double *ImaginaryIn, /* array of input's imag samples */double *RealOut, /* array of output's reals */double *ImaginaryOut ); /* array of output's imaginaries */void fft_float (
unsigned NumSamples, /* must be a power of 2 */int InverseTransform, /* 0=forward FFT, 1=inverse FFT */float *RealIn, /* array of input's real samples */float *ImaginaryIn, /* array of input's imag samples */float *RealOut, /* array of output's reals */float *ImaginaryOut ); /* array of output's imaginaries */int IsPowerOfTwo ( unsigned x );
unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo );
unsigned ReverseBits ( unsigned index, unsigned NumBits );
/*
** The following function returns an "abstract frequency" of a
** given index into a buffer with a given number of frequency samples.
** Multiply return value by sampling rate to get frequency expressed in Hz.
*/double Index_to_frequency ( unsigned NumSamples, unsigned Index );
#ifdef __cplusplus
}
#endif/*--- end of file fourier.h ---*/
/bur
Re[2]: Что я делаю не так, реализуя быстрое преобразование Ф
Здравствуйте, mab, Вы писали:
mab>должны быть так размещены. Просто слабо верится, что операции рекурсивного FFT и битового mab>обращения коммутируют...
Конечно, коммутируют. Бит реверсивный код это не что иное, как выделение четного и нечетного, потом четного из четного, нечетного из четного, четного из нечетного и т.д. и т.п. — именно так ведут себя биты, если смотреть с младшего. Например, самый младший отмечает обычный чет-нечет.
Я об этом прочитал в документе The FFT demystified. Все очень понятно разжевано. Сейчас мне пришла в голову мысль, что ошибка действительно где-то в индексации. Если я беру веса какой-то неправильной компоненты спектра и рисую с ее помощи 50-герцовую синусоиду — она и впрямь должна прыгать...