Heja, mam problem z transforamta fouriera. Zaimpelementowalem algorytm radix 2 ale nie dziala on do konca dobrze. Tzn oglony charakter jest zachowany ale widac roznice w wartosciach. Ponizej przedstawiam program prawidlowy z matlaba i ten moj z C. Prosze o pomoc w znalezieniu bledu.
Matlab:
N=256;
x=0:N-1;
typBitReverse=2;
typFFT=2;
t=0:0.1:N/10-0.1;
x=sin(t*pi)+0.2*sin(t*pi*7);
xc=x;
a=1;
for b=1:N-1
if (b<a)
T=x(a); x(a)=x(b); x(b)=T;
end
c=N/2;
while (c<a)
a=a-c;
c=c/2;
end
a=a+c;
end
Moje progi w C:
// program realizuje test transforamaty fouriera i zapis do plikow wynikowych
#include <stdio.h>
#include <math.h>
#include "fourier.h"
#define N 256
#define LPI 3.14159265358979323846
int main(void)
{
int i;
double sinus[N],real[N],imag[N];
FILE *plik_dane;
FILE *plik_realis;
FILE *plik_imag;
//--------------------------------------------Utworzenie sygnalu sinusoidalnego
for(i=0;i<N;i++)
{
sinus[i]=sin(LPI*i*0.1)+0.2*(sin(LPI*i*7*0.1));
}
fourier(sinus,imag);
//----------------------------------------zapis danych do pliku z tablic
plik_realis=fopen("C:\\pliki/realis.txt", "w");
for(i=0;i<N;i++)
{
fprintf(plik_realis,"%lf\n",sinus[i] );
}
fclose(plik_realis);
plik_imag=fopen("C:\\pliki/imag.txt", "w");
for(i=0;i<N;i++)
{
fprintf (plik_imag,"%lf\n",imag[i] );
}
fclose(plik_imag);
return 0;
}
fourier.h
//realizacja fouriera dla 256 probek metoda radix-2
#define NN 256 //liczba probek
#include <math.h>
#define LPI 3.14159265358979323846
void fourier(double x[],double y[])
{
int a=1,z,licznik;
int b,c,d,e,g,L,M,mm,i;
double T;
double W_r,W_i,Wi_r,Wi_i,T_r,T_i;
double licz[1024];
licznik=0;
FILE *wyn;
//-----------------------------------------przestawianie elementow
for(b=1;b<NN;b++)
{
if (b<a)
{
T=x[a-1];
x[a-1]=x[b-1];
x[b-1]=T;
}
c=NN/2;
while (c<a)
{
a-=c;
c/=2;
}
a+=c;
}
//--------------------------------------FFT
for(e=0;e<NN;e++)
{
y[e]=0;
}
for(e=1;e<9;e++)
{
L = pow(2,e);
M = L/2;
Wi_r=1;
Wi_i=0;
W_r=cos(2*LPI/L);
W_i=-(sin(2*LPI/L));
for(mm=1;mm<M+1;mm++)
{
for(g=mm;g<NN+1;g+=L)
{
d = g+M;
//T=x(d)*Wi
T_r = x[d-1]*Wi_r-y[d-1]*Wi_i;
T_i = x[d-1]*Wi_i+y[d-1]*Wi_r;
// x(d)=x(g)-T
x[d-1] = x[g-1]-T_r;
y[d-1] = y[g-1]-T_i;
//x(g)=x(g)+T
x[g-1] = x[g-1]+T_r;
y[g-1] = y[g-1]+T_i;
}
Wi_r=Wi_r*W_r-Wi_i*W_i;
Wi_i=Wi_r*W_i+Wi_i*W_r;
}
}
}
Z gory dzieki za pomoc jest mi to naprawde potrzebne, pozdry