Справочник функций

Ваш аккаунт

Войти через: 
Забыли пароль?
Регистрация
Информацию о новых материалах можно получать и без регистрации:

Почтовая рассылка

Подписчиков: -1
Последний выпуск: 19.06.2015

FFT. Помогите найти ошибку

75K
25 декабря 2011 года
Teddy_bear
8 / / 29.11.2011
Вместо ожидаемого пика на 1 кГц наблюдаю следующее непойми-что:




Код:
//---------------------------------------------------------------------------
#include <math.h>
#include <vcl.h>
#include "open_wav.cpp"

#pragma hdrstop

#include "fft.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;

#define N 8192  //кол-во элементов в окне

//Функция, инвертирующая биты числа. С ее помощью делается
//двоично-инверсная перестановка
unsigned int BinaryInverse(unsigned int InValue, unsigned int L)
{
static unsigned char reverse256[]= {
    0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0,
    0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
    0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,
    0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
    0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,
    0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
    0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,
    0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
    0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,
    0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
    0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,
    0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
    0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,
    0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
    0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,
    0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
    0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,
    0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
    0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,
    0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
    0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,
    0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
    0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED,
    0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
    0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,
    0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
    0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,
    0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
    0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,
    0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
    0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,
    0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF,
};

 unsigned int OutValue;
 unsigned char *pInValue = (unsigned char*)&InValue;
 unsigned char *pOutValue = (unsigned char*)&OutValue;

 pOutValue[0]=reverse256[pInValue[3]];
 pOutValue[1]=reverse256[pInValue[2]];
 pOutValue[2]=reverse256[pInValue[1]];
 pOutValue[3]=reverse256[pInValue[0]];

 OutValue >>=(32-L);
 return (OutValue);
}


struct complex {double re,im;};
short *in; //указатель на входные данные
complex A[N],B[N]; //массивы для промежуточных результатов

//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}

//Читаем данные из wav-файла в память, потом с "двоично-инверсной" перестановкой
//записываем в массив. Выводим в Chart в "нормальной" последовательности.
void __fastcall TForm1::Button2Click(TObject *Sender)
{
int i;
FILE *inptr;

OpenDialog1->Execute();
char *a=OpenDialog1->FileName.c_str();
inptr = fopen(a,"rb");
ReadWaveHeader(inptr);
in=(short*)malloc(N*2);

for(i=0;i<N;i++)
{
 fread(in+i,sizeof(short),1,inptr);
 A[BinaryInverse(i,13)].re=in;
 A[BinaryInverse(i,13)].im=0;
}
for(i=0;i<N;i++) Chart1->Series[0]->AddXY(i,A[BinaryInverse(i,13)].re,"",clGreen);

}
//---------------------------------------------------------------------------

void __fastcall TForm1::Button1Click(TObject *Sender)
{
int i,j,Nsub=1; //Nsub - длина подпоследовательности в общем массиве
double ampl[N],n,k,W;
bool a=false,b=false;
   
// основной цикл. Данные каждого следующего уровня объединения
// формируются поочередно из массивов А и В, и записываься в другой
// массив, из к-рого формируется следующий уровень и т.д. до тех пор,
// пока Nsub не станет равно N.              
while(1)
 {
  if(!(Nsub==N))Nsub*=2;
  else {a=true; break;}

  for(i=0,n=0;i<=N-Nsub;i+=Nsub,n+=Nsub)
  {
   for(k=i,j=i;j<(Nsub/2)+i;k++,j++)
   {
    W=2*M_PI*(k-n)/Nsub;
    B[j].re=A[j].re+cos(W)*A[j+Nsub/2].re;
    B[j].im=A[j].im-sin(W)*A[j+Nsub/2].im;
   }
   for(k,j;j<Nsub+i;k++,j++)
   {
    W=2*M_PI*(k-n-0.5*Nsub)/Nsub;
    B[j].re=A[j-Nsub/2].re-cos(W)*A[j].re;
    B[j].im=A[j-Nsub/2].im+sin(W)*A[j].im;
   }
  }

  if(!(Nsub==N))Nsub*=2;
  else {b=true; break;}

  for(i=0,n=0;i<=N-Nsub;i+=Nsub,n+=Nsub)
  {
   for(k=i,j=i;j<(Nsub/2)+i;k++,j++)
   {
    W=2*M_PI*(k-n)/Nsub;
    A[j].re=B[j].re+cos(W)*B[j+Nsub/2].re;
    A[j].im=B[j].im-sin(W)*B[j+Nsub/2].im;
   }
   for(k,j;j<Nsub+i;k++,j++)
   {
    W=2*M_PI*(k-n-0.5*Nsub)/Nsub;
    A[j].re=B[j-Nsub/2].re-cos(W)*B[j].re;
    A[j].im=B[j-Nsub/2].im+sin(W)*B[j].im;
   }
  }
 }
 for(i=3, k=3; i<N; i++, k++)
 {
  if(b)
  {
   ampl=sqrt(B.re*B.re + B.im*B.im)/N;
   Chart2->Series[0]->AddXY(k*44100/N,ampl,"",clGreen);
  }
  if(a)
  {
   ampl=sqrt(A.re*A.re + A.im*A.im)/N;
   Chart2->Series[0]->AddXY(k*44100/N,ampl,"",clGreen);
  }
 }
 a=false;
 b=false;
}


Вот значения поворачивающих коэф-тов каждого уровня для N=8:
1
-1
1
-1
1
-1
1
-1
-----------------
1
6,12303176911189E-17
-1
-6,12303176911189E-17
1
6,12303176911189E-17
-1
-6,12303176911189E-17
-----------------
1
0,707106781186548
6,12303176911189E-17
-0,707106781186547
-1
-0,707106781186548
-6,12303176911189E-17
0,707106781186547
-----------------
Начальная перестановка и последующая схема движения эл-тов проверены и работают правильно (как описано тут:http://www.dsplib.ru/content/thintime/thintime.html)
Что не так?
75K
13 января 2012 года
Teddy_bear
8 / / 29.11.2011
Упростил реализацио основного цикла (см. код), но проблема осталась. Подозреваю, что что-то не так с комплексными операциями. Но что именно? Помогите плиз, уже все мозги сломал:confused:

Код:
{
 int Nsub,n,k,c=1;
 long double ampl[N],W,temp,i,j;

 for(Nsub=2;Nsub<=N;Nsub*=2)
 {
  for(n=0,i=0;n<N;n+=Nsub,i+=Nsub)
  for(k=n,j=i;k<(Nsub/2)+n;k++,j++)
  {
   W=(2*M_PI*(i-j)/Nsub);
   temp=A[k].re;

   A[k].re+=cos(W)*A[k+Nsub/2].re;
   A[k+Nsub/2].re=temp-cos(W)*A[k+Nsub/2].re;

   A[k].im+=sin(W)*A[k+Nsub/2].im;
   A[k+Nsub/2].im=temp-sin(W)*A[k+Nsub/2].im;
 
  }
 }

}
79K
15 января 2012 года
vital792
5 / / 15.01.2012
К сожалению билдеров не имею, проверял на gcc. Честно говоря особо не вникал в основной цикл, в чем была ошибка не скажу. Сразу заменил своим и заработало.

Цитата:

#include "math.h"
#include "stdio.h"

#define N 8192

struct complex {double re,im;};

unsigned int BinaryInverse(unsigned int InValue, unsigned int L)
{
static unsigned char reverse256[]= {
0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0,
0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,
0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,
0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,
0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,
0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,
0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,
0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,
0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,
0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,
0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,
0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED,
0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,
0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,
0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,
0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,
0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF,
};

unsigned int OutValue;
unsigned char *pInValue = (unsigned char*)&InValue;
unsigned char *pOutValue = (unsigned char*)&OutValue;

pOutValue[0]=reverse256[pInValue[3]];
pOutValue[1]=reverse256[pInValue[2]];
pOutValue[2]=reverse256[pInValue[1]];
pOutValue[3]=reverse256[pInValue[0]];

OutValue >>= (32-L);
return OutValue;
}

void fft(complex *A)
{
/*
int Nsub,n,k,c=1;
double ampl[N], W, temp, i, j;

for(Nsub=2; Nsub<=N; Nsub*=2)
{
for(n=0, i=0;n<N; n+=Nsub, i+=Nsub)
for(k=n, j=i; k<(Nsub/2)+n; k++, j++)
{
W = (2*M_PI*(i-j)/Nsub);
temp = A[k].re;

A[k].re += cos(W)*A[k+Nsub/2].re;
A[k+Nsub/2].re = temp-cos(W)*A[k+Nsub/2].re;

A[k].im += sin(W)*A[k+Nsub/2].im;
A[k+Nsub/2].im = temp-sin(W)*A[k+Nsub/2].im;
}
}
*/
double *out_Buffer = (double*)A;
double wr, wi, arg, *p1, *p2, temp;
double tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
long i, bitm, j, le, le2, k;
long lim = (long)(log((double)N)/log(2.)+.5);
for (k = 0, le = 2; k < lim; k++)
{
le <<= 1;
le2 = le>>1;
ur = 1.0;
ui = 0.0;
arg = M_PI / (le2>>1);
wr = cos(arg);
wi = -sin(arg);
for (j = 0; j < le2; j += 2)
{
p1r = out_Buffer+j; p1i = p1r+1;
p2r = p1r+le2; p2i = p2r+1;
for (i = j; i < 2*N; i += le)
{
tr = *p2r * ur - *p2i * ui;
ti = *p2r * ui + *p2i * ur;
*p2r = *p1r - tr; *p2i = *p1i - ti;
*p1r += tr; *p1i += ti;
p1r += le; p1i += le;
p2r += le; p2i += le;
}
tr = ur*wr - ui*wi;
ui = ur*wi + ui*wr;
ur = tr;
}
}
}

int main()
{
double *signal = new double[N];
for(int i=0;i<N;i++)
signal = sin(2*M_PI*1000*i/44100); // 1 kHz Fs = 44100

complex A[N];
for(int i=0; i<N; i++)
{
A[BinaryInverse(i,13)].re=signal;
A[BinaryInverse(i,13)].im=0;
}

fft(A);

double *AmpSpec = new double [N];
for(int i=0; i<N; i++)
AmpSpec = sqrt(A.re*A.re+A.im*A.im);

FILE *f = fopen("spec.txt", "w+");
double re, im;
for(unsigned i=0; i<N; i++)
{
re = A.re;
im = A.im;
fprintf(f,"%.7e%c", sqrt(re*re + im*im), ' ');
}
fclose(f);

delete [] signal;
delete [] AmpSpec;
return 0;
}


в функции fft ваш цикл заменен на мой. Результат сохранил в текстовый файлик, смотрел в матлабе - все честно, пик на частоте 1 кГц

79K
15 января 2012 года
vital792
5 / / 15.01.2012
Кстати забыл сказать: очень рекомендую пользоваться библиотекой fftw для подобных задач. Быстрее реализации просто не бывает. Или хотя бы alglib. Fft там реализовано похуже, зато много других вкусных вещей.
75K
16 января 2012 года
Teddy_bear
8 / / 29.11.2011
спасибо, буду разбираться. Может посоветуете, что почитать по этой теме? Чтоб поменьше формул и не слишком заумно)
79K
16 января 2012 года
vital792
5 / / 15.01.2012
Чтобы разобраться с fft советую почитать книгу Ахмед, Рао "Ортогональные преобразования при обработке цифровых сигналов" - описано досконально, хотя книга явно не для новичков. И формул там полно, но именно ффт прям разжевано. А для новичков, на мой взгляд самый легкий и понятный учебник по цос - Сергиенко "Цифровая обработка сигналов". Еще можно почитать мануалы и исходники alglib, спасибо код открыт
Реклама на сайте | Обмен ссылками | Ссылки | Экспорт (RSS) | Контакты
Добавить статью | Добавить исходник | Добавить хостинг-провайдера | Добавить сайт в каталог