Hi everybody again!! and thanx! everyone who have the time to help me..
Tim Wescott: The first thing you suspect may be my problem... but i dont
know how to determine what scale use Agilent VEE ??..
by the second thought.. I don´t know how to do it with the fftw library.
May be this is my problem... but i don't know how to do it with fftw
library... which method i have to use.. to calculate fft with complex input
and expect real output?. If you know!! PLEASE HELP!!
eric: thanks for your post!! you did well!! but... i cant get the result
that i expect.. :S
How do you get this results?
0.00000000 : 0.00000000
0.00000000 : 0.00000000
0.00000000 : 0.00000000
0.00000000 : 0.00000000
-0.14420960 : 0.23532869
0.06358576 : -0.19604607
0.00000000 : 0.00000000
0.00000000 : 0.00000000
0.00000000 : 0.00000000
0.00000000 : 0.00000000
by the way... i use those results in my ifft program and... i have
different output than yours. So may be my code wrong?? Here is the code...
And in the bottom i put the link to the files. CH00 and CH01 are the input
to my program.
IN IFFT(agilent).txt and IFFT out (agilent).txt: are the input to ifft in
agilent and the output from it.
Note: Agilent make the same calculations that my program... and have the
same result until ifft. After ifft, everything goes to hell. :(
CODE:
#include <fftw3.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <conio.h>
#define SIZE 2048
#define span 4
#define fcorte 1.4
using namespace std;
void mag (fftw_complex *result, long asize);
void phase ( fftw_complex *result, long asize);
void magpha ( fftw_complex *result, long asize);
long mapi ( double vspan, double asize);
void dco (fftw_complex *result, fftw_complex *result2, long vmapi);
void alloc ( fftw_complex *result, long asize);
void fsave (fftw_complex *result, long asize);
double pha[1+(SIZE/2)];
double magn[1+(SIZE/2)];
long hifrec, lowfrec;
int main( int argc, char** argv )
{
double asize,vspan,vmapi;
char tmp [100];
fftw_complex *data, *fft_result, *ifft_result,*data2, *fft_result2,
*ifft_result2;
fftw_plan plan_forward, plan_backward, plan_forward2,
plan_backward2;
int i,n;
/* IN/OUT */
data=( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * SIZE );
data2=( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * SIZE );
fft_result=( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * SIZE
);
fft_result2=( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * SIZE
);
ifft_result=( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * SIZE
);
ifft_result2=( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * SIZE
);
/* make plans
*/
plan_forward = fftw_plan_dft_1d( SIZE, data, fft_result,
FFTW_FORWARD, FFTW_ESTIMATE );
plan_forward2 = fftw_plan_dft_1d( SIZE, data2, fft_result2,
FFTW_FORWARD, FFTW_ESTIMATE );
plan_backward = fftw_plan_dft_1d( SIZE, fft_result, ifft_result,
FFTW_BACKWARD, FFTW_ESTIMATE );
plan_backward2= fftw_plan_dft_1d( SIZE, fft_result2, ifft_result2,
FFTW_BACKWARD, FFTW_ESTIMATE );
/* Read CH1 & CH2*/
ifstream archivo;
archivo.open ("ch00.txt", ifstream::in);
if (archivo.is_open())
{
n = 0;
while (SIZE>=n)
//while (!archivo.eof()
{
archivo.getline (tmp,100);
data[n][0] = atof (tmp);
data[n][1] = 0.0;
n++;
}
}
else
{
cout << "Error al abrir el archivo";
}
archivo.close ();
ifstream archiv;
archiv.open ("ch01.txt", ifstream::in);
if (archiv.is_open())
{
n = 0;
while (SIZE>=n)
//while (!archivo.eof()
{
archiv.getline (tmp,100);
data2[n][0] = atof (tmp);
data2[n][1] = 0.0;
n++;
}
}
else
{
cout << "Error al abrir el archivo";
}
archiv.close ();
/* print initial data
*/
for( i = 0 ; i < (SIZE) ; i++ ) {
cout << i << ": " << data[i][0] << " : "<< data[i][1] <<endl;
}
getch();
for( i = 0 ; i < (SIZE) ; i++ ) {
cout << i << ": " << data2[i][0] << " : "<< data2[i][1] <<endl;
}
getch();
/*execute ffts*/
fftw_execute( plan_forward ); //Ejecuta el plan A
fftw_execute( plan_forward2); //Ejecuta el plan A.1
/* print fft result
*/
for( i = 0 ; i < SIZE/2+1 ; i++ ) {
cout << i << ": " << fft_result[i][0] << " : " << fft_result[i][1]
<< endl;
}
getch();
for( i = 0 ; i < SIZE/2+1 ; i++ ) {
cout << i << ": " << fft_result2[i][0] << " : " << fft_result2[i][1]
<< endl;
}
getch();
asize=SIZE/2+1;
mag (fft_result,asize);
phase(fft_result,asize);
magpha(fft_result,asize);
mag (fft_result2,asize);
phase(fft_result2,asize);
magpha(fft_result2,asize);
for( i = 0 ; i < asize ; i++ ) {
cout << i << ": " << fft_result[i][0] << " : " << fft_result[i][1]
<< endl;
}
getch();
for( i = 0 ; i < asize ; i++ ) {
cout << i << ": " << fft_result2[i][0] << " : " << fft_result2[i][1]
<< endl;
}
getch();
vspan = asize/span; /*Calcula el Stop del Star/Stop del agilent*/
vmapi = mapi(vspan,asize);
dco (fft_result,fft_result2,vmapi);
alloc(fft_result,asize);
alloc(fft_result2,asize);
fsave(fft_result,asize);
for( i = 0 ; i < asize ; i++ ) {
cout << i << ": " << fft_result[i][0] << " : " << fft_result[i][1]
<< endl;
}
getch();
for( i = 0 ; i < asize ; i++ ) {
cout << i << ": " << fft_result2[i][0] << " : " << fft_result2[i][1]
<< endl;
}
getch();
fftw_execute( plan_backward ); //Ejecuta el plan B
fftw_execute( plan_backward2 ); //Ejecuta el plan B.1
/* print ifft result
*/
for( i = 0 ; i < SIZE ; i++ ) {
cout << i << ": " << ifft_result[i][0] / SIZE << " : " <<
fft_result[i][1] << endl;
}
getch();
for( i = 0 ; i < SIZE ; i++ ) {
cout << i << ": " << ifft_result2[i][0] / SIZE << " : " <<
fft_result2[i][1] << endl;
}
getch();
/* free memory
*/
fftw_destroy_plan( plan_forward );
fftw_destroy_plan( plan_backward );
fftw_destroy_plan( plan_forward2);
fftw_destroy_plan( plan_backward2);
fftw_free( data );
fftw_free( fft_result );
fftw_free( ifft_result );
fftw_free( data2 );
fftw_free( fft_result2 );
fftw_free( ifft_result2 );
return 0;
}
//////////////////////////////////////////////////PHASE////////////////////////
/* this function calculate phase.*/
void phase ( fftw_complex *result, long asize)
{
long n;
double atan;
for ( n = 0; n < asize ; n++ )
{
atan=atan2(result[n][1], result[n][0]);
pha[n]=atan*(180/3.14159265);
}
}
//////////////////////////////////////////////////MAGNITUD////////////////
/* this function calculate sqrt(a^2+b^2).*/
void mag (fftw_complex *result, long asize)
{
long n;
for ( n = 0; n < asize ; n++ )
{
magn[n]=sqrt(pow(result[n][0],2) + pow(result[n][1],2));
}
}
//////////////////////////////////////////////MAG-PHASE///////////////////
/*this function use mag and phase to make a new array with magnitude and
phase data*/
void magpha ( fftw_complex *result, long asize)
{
long n;
double atan;
for ( n = 0; n < asize ; n++ )
{
result[n][0] = magn[n];
result[n][1] = pha[n];
}
}
///////////////////////////////MAP A INDEX//////////////////////////
long mapi ( double vspan, double asize)
{
double resultado, intpart, fractpart;
resultado = fcorte*(asize/vspan);
fractpart = modf(resultado,&intpart);
if (fractpart>0.5)
{
intpart = intpart + 1;
}
return intpart;
}
/////////////////////////////////////////CUT///////////////////////////////////
/*this function and "map a index" finds hifrec and lowfrec. With this
values i set up the new array*/
void dco (fftw_complex *result, fftw_complex *result2, long vmapi)
{
long a, pico;
int hifrec1, hifrec2, lowfrec1, lowfrec2, n, flag;
//maxIndex([x[a-1],x[a],x[a+1]])+a-1
if (result[vmapi-1][0]>=result[vmapi][0]) a=0;
else a=1;
if (result[vmapi][0]>=result[vmapi+1][0]) a=1;
else a=2;
pico = a + vmapi - 1;
hifrec1 = pico;
hifrec2 = pico;
lowfrec1 = pico;
lowfrec2 = pico;
n = pico;
flag = 0;
while (flag==0)
{
if (result[n+1][0] <= result[n][0]) n=n+1;
else
{
lowfrec1=n;
flag=1;
}
}
flag=0;
n = pico;
while (flag==0)
{
if (result2[n+1][0] <= result2[n][0]) n=n+1;
else
{
lowfrec2=n;
flag=1;
}
}
flag = 0;
n = pico;
while (flag==0)
{
if (result[n-1][0] <= result[n][0])
{
n=n-1;
if (n==0)
{
hifrec1=n;
flag=1;
}
}
else
{
hifrec1=n;
flag=1;
}
}
flag = 0;
n = pico;
while (flag==0)
{
if (result2[n-1][0] <= result2[n][0])
{
n=n-1;
if (n==0)
{
hifrec2=n;
flag=1;
}
}
else
{
hifrec2=n;
flag=1;
}
}
if (lowfrec1>lowfrec2) lowfrec=lowfrec2;
else lowfrec=lowfrec1;
if (hifrec1>hifrec2) hifrec=hifrec1;
else hifrec=hifrec2;
}
//////////////////////////////ALLOC////////////////////////////////
/*this function make a new array with all values at 0.0 and copy the
hifrec and lowfrec from the original fft array (magnitude,phase)*/
void alloc ( fftw_complex *result, long asize)
{
int n=0;
double tmp1,tmp2;
while (asize>=n)
{
if (n==hifrec){}
else
{
if (n==lowfrec){}
else
{
result[n][0] = 0.0;
result[n][1] = 0.0;
}
}
n=n+1;
}
/*
tmp1=result[hifrec][0];
tmp2=result[hifrec][1];
result[hifrec][0]=result[lowfrec][0];
result[hifrec][1]=result[lowfrec][1];
result[lowfrec][0]=tmp1;
result[lowfrec][1]=tmp2;*/
}
//////////////////////////////SAVE////////////////////////
void fsave (fftw_complex *result, long asize)
{
int n;
ofstream fsalida;
fsalida.open("salida.TXT");
for ( n = 0; n < asize ; n++ )
{
fsalida << result[n][0] << "," << result[n][1] << endl;
}
fsalida.close();
}
Files:
IN IFFT(agilent).txt:
http://www.4shared.com/file/55678120/d91bbce1/IN_IFFT_agilent_.html
IFFT out (agilent).txt:
http://www.4shared.com/file/55678126/307819d4/IFFT_out__agilent_.html
CHO1:
http://www.4shared.com/file/55678122/3715ddcd/CH01.html
CH00:
http://www.4shared.com/file/55678124/de7678f8/CH00.html
Hope you understand my English and THANKS FOR EVERY HELP... I NEED IT...!
Have a nice day!!