Новичок
Профиль
Группа: Участник
Сообщений: 1
Регистрация: 17.6.2008
Репутация: нет
Всего: нет
|
Уважаемые формучане, пишу к вам с просьбой о помощи. Сам в институте еще не учился, поэтому про сабж мало чего знаю. Но в недавнее время очень понадобилось это самое преобразование. На вход ему посылается массив, заполненый синусойдной функцией, от 0 до 255. после преобразование получается массив, заполненый уж больно большими числами. При обратном преобразовании получается либо вообще непонятно что, либо исходный массив, но с некоторыми помехами. Насколько понимаю - при обратном преобразовании должен получатся совершенно точный начальный сигнал. вот код найденый гуглем, и хоть както трансформирующий туда-сюда массив. Нужно избавиться от помех, и сделать чтобы преобразование было правильным. | Код |
using namespace std;
#define CACHE_HALF 65536
#define CONST_PI 3.1415926535897932384626433832
#define CONST_SQRT_2 0.7071067811865475244008443621
#define CONST_SQRT2 1.4142135623730950488016887242
#define max(x,y) ((x) > (y) ? (x) : (y))
#define min(x,y) ((x) < (y) ? (x) : (y))
typedef double real;
typedef unsigned long ulong;
typedef unsigned short ushort;
class Complex {
public:
real r, i;
Complex(void) { }
Complex(real a, real b) { r=a; i=b; }
inline const Complex operator+(const Complex &c) const {
return Complex( r + c.r, i + c.i);
}
inline const Complex operator-(const Complex &c) const {
return Complex( r - c.r, i - c.i);
}
inline const Complex operator*(const Complex &c) const {
return Complex( r*c.r - i*c.i, r*c.i + i*c.r);
}
inline const Complex operator/(const real &divisor) const {
return Complex( r/divisor, i/divisor);
}
};
inline const Complex conj(const Complex &c) {
return Complex( c.r, -c.i);
}
real SineTable[80];
int Log2(ulong Num) {
int x=-1;
if (Num==0) return 0;
while (Num) {x++;Num/=2;}
return x;
}
#define INIT_TRIG(LENGTH) \
ulong x=Log2(LENGTH); \
real Sin0=SineTable[x]; \
real Cos0=SineTable[x+1]; \
Cos0=-2.0*Cos0*Cos0; \
real Sin=Sin0,Cos=1.0+Cos0;
#define NEXT_TRIG_POW { \
real temp=Cos; \
Cos = Cos*Cos0 - Sin*Sin0 + Cos; \
Sin = Sin*Cos0 + temp*Sin0 + Sin; \
}
inline ulong rev_next(ulong r, ulong n) {
do {
n = n >> 1;
r = r^n;
} while ( (r&n) == 0);
return r;
}
// FFTReOrder для действительных векторов
void FHTReOrder(real *Data, ulong Len) {
real temp;
if (Len <= 2) return;
ulong r=0;
for ( ulong x=1; x<Len; x++) {
r = rev_next(r, Len);
if (r>x) { temp=Data[x]; Data[x]=Data[r]; Data[r]=temp; }
}
}
void CreateSineTable(ulong Len) {
int x=0; ulong P=1;
while (P<=Len*4) {
SineTable[x]=sin(CONST_PI/P);
P*=2;
x++;
}
}
#define TRIG_VARS \
ulong TLen,TNdx;int TDir; \
Complex PRoot,Root;
#define INIT_TRIG(LENGTH,DIR) \
TNdx=0;TLen=(LENGTH);TDir=(DIR); \
PRoot.r=1.0;PRoot.i=0.0; \
Root.r=sin(CONST_PI/((LENGTH)*2.0));\
Root.r=-2.0*Root.r*Root.r; \
Root.i=sin(CONST_PI/(LENGTH))*(DIR);
#define NEXT_TRIG_POW \
if (((++TNdx)&15)==0) { \
real Angle=(CONST_PI*(TNdx))/TLen; \
PRoot.r=sin(Angle*0.5); \
PRoot.r=1.0-2.0*PRoot.r*PRoot.r; \
PRoot.i=sin(Angle)*(TDir); \
} else { \
Complex Temp; \
Temp=PRoot; \
PRoot = PRoot*Root; \
PRoot = PRoot+Temp; \
}
void FFTReOrder(Complex *Data, ulong Len) {
Complex temp;
if (Len <= 2) return;
ulong r=0;
for ( ulong x=1; x<Len; x++) {
r = rev_next(r, Len);
if (r>x) { temp=Data[x]; Data[x]=Data[r]; Data[r]=temp; }
}
}
void FFT_T(Complex *Data, ulong Len, int Dir) {
ulong k;
TRIG_VARS;
if (Len <= (CACHE_HALF/sizeof(Complex)) ) {
IFFT_T(Data, Len,Dir);
return;
}
Len /= 2;
INIT_TRIG(Len, Dir);
FFT_T(Data, Len,Dir);
FFT_T(Data+Len,Len,Dir);
for (k=0; k<Len; k++) {
Complex b,c;
b=Data[k];
c = Data[k+Len] * PRoot;
Data[k] = b + c;
Data[k+Len] = b - c;
NEXT_TRIG_POW;
}
}
float * RealFFT(float *ddata, ulong Len, int Dir) {
ulong i, j;
Complex *Data=(Complex*)ddata;
TRIG_VARS;
Len /= 2;
if (Dir > 0) {
FFTReOrder(Data,Len);
FFT_T(Data,Len,1);
}
INIT_TRIG(Len,Dir);
NEXT_TRIG_POW;
for (i = 1, j = Len - i; i < Len/2; i++, j--) {
Complex p1,p2,t;
t = conj(Data[j]);
p1 = Data[i] + t;
p2 = Data[i] - t;
p2 = p2 * PRoot;
t = Complex(-Dir*p2.i,Dir*p2.r);
Data[i] = p1 - t;
Data[j] = p1 + t;
Data[j] = conj(Data[j]);
Data[i] = Data[i]/2;
Data[j] = Data[j]/2;
NEXT_TRIG_POW;
}
{ real r,i;
r=Data[0].r;i=Data[0].i;
Data[0] = Complex(r+i,r-i);
}
if (Dir < 0) {
Data[0] = Data[0]/2.0;
FFTReOrder(Data,Len);
FFT_T(Data,Len,-1);
}
return (float *)Data;
}
................
float *x = new float[P_SIZE];
x = RealFFT(c,P_SIZE,1); //с - массив для трансформации
x = RealFFT(x,P_SIZE,-1);//при обратном преобразовании получаются гадости
|
Заранее спасибо. |
C++ и преобразование фурье, хелп 
Загрузка ...
