//---------------------------------------------------------------------------
#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;
}
FFT. Помогите найти ошибку
Код:
Вот значения поворачивающих коэф-тов каждого уровня для 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)
Что не так?
Код:
{
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;
}
}
}
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;
}
}
}
Цитата:
#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 кГц
Кстати забыл сказать: очень рекомендую пользоваться библиотекой fftw для подобных задач. Быстрее реализации просто не бывает. Или хотя бы alglib. Fft там реализовано похуже, зато много других вкусных вещей.
спасибо, буду разбираться. Может посоветуете, что почитать по этой теме? Чтоб поменьше формул и не слишком заумно)
Чтобы разобраться с fft советую почитать книгу Ахмед, Рао "Ортогональные преобразования при обработке цифровых сигналов" - описано досконально, хотя книга явно не для новичков. И формул там полно, но именно ффт прям разжевано. А для новичков, на мой взгляд самый легкий и понятный учебник по цос - Сергиенко "Цифровая обработка сигналов". Еще можно почитать мануалы и исходники alglib, спасибо код открыт