Hi everybody... i am new in signal processing.
something about me... i am a student of software Engineering, My name is
Juan and i live in Argentina.
I am migrating one application developed in VEE (Agilent) to c++. This is
my new task in my work. This application analyzes signals and then... use
that info to set some alarms.. Nothing complicated... so far..
everything goes fine.. until i have to use ifft. I did make a fft
algorithm to do everything... but when i do the ifft... the output data
from VEE (agilent) was not the same as my algorithm. So... I decide to
implement FFTW. And i have the same problem...
i have this input data. (size = 2048)
Agilent(VEE) c++
magnitude phase
0 , @0 0 , 0
0 , @0 0 , 0
0 , @0 0 , 0
0 , @0 0 , 0
0.276 , @121.5 0.276 , 121.5
0.2061, @-72.03 0.2061, -72.03
0 , @0 0 , 0
0 , @0 0 , 0
0 , @0 0 , 0
0 , @0 0 , 0
... ... ... ...
and the output data from the VEE (agilent)
-7.860034643685586E-005
-7.871299096183408E-005
-7.88070606286872E-005
-7.888302098921511E-005
-7.894133976812804E-005
-7.898248647698485E-005
-7.900693202945827E-005
-7.901514835816201E-005
-7.900760803327236E-005
-7.898478388317609E-005
-7.894714861737379E-005
-7.889517445186702E-005
-7.882933273725408E-005
-7.875009358975982E-005
-7.865792552541923E-005
-7.855329509763516E-005
-7.843666653832641E-005
..
and this is the output from c++ (FFTW)
0: 0.000235391 : 0.0241332
1: -0.000100984 : 0.0241359
2: -0.000437429 : 0.0241348
3: -0.000773924 : 0.0241297
4: -0.00111045 : 0.0241208
5: -0.00144697 : 0.024108
6: -0.00178348 : 0.0240913
..
I think that could be the type of the input data... because in agilent is
a waveform. And in my c++ application is real data..
In/out (FFTW)
data=(fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * (SIZE/2+1) );
ifft_result=(fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * (SIZE)
);
I am relative new in the Signal Processing... so i would need some help...
and if you know some books.. And excuse me for my English :D
have a nice day...:cool:
Problem with IFFT, any idea?
Started by ●July 17, 2008
Reply by ●July 17, 20082008-07-17
rapidiablo wrote:> Hi everybody... i am new in signal processing. > > something about me... i am a student of software Engineering, My name is > Juan and i live in Argentina. > > I am migrating one application developed in VEE (Agilent) to c++. This is > my new task in my work. This application analyzes signals and then... use > that info to set some alarms.. Nothing complicated... so far.. > > everything goes fine.. until i have to use ifft. I did make a fft > algorithm to do everything... but when i do the ifft... the output data > from VEE (agilent) was not the same as my algorithm. So... I decide to > implement FFTW. And i have the same problem... > > i have this input data. (size = 2048) > > Agilent(VEE) c++ > > magnitude phase > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0.276 , @121.5 0.276 , 121.5 > 0.2061, @-72.03 0.2061, -72.03 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > ... ... ... ... > > and the output data from the VEE (agilent) > > -7.860034643685586E-005 > -7.871299096183408E-005 > -7.88070606286872E-005 > -7.888302098921511E-005 > -7.894133976812804E-005 > -7.898248647698485E-005 > -7.900693202945827E-005 > -7.901514835816201E-005 > -7.900760803327236E-005 > -7.898478388317609E-005 > -7.894714861737379E-005 > -7.889517445186702E-005 > -7.882933273725408E-005 > -7.875009358975982E-005 > -7.865792552541923E-005 > -7.855329509763516E-005 > -7.843666653832641E-005 > .. > > and this is the output from c++ (FFTW) > > 0: 0.000235391 : 0.0241332 > 1: -0.000100984 : 0.0241359 > 2: -0.000437429 : 0.0241348 > 3: -0.000773924 : 0.0241297 > 4: -0.00111045 : 0.0241208 > 5: -0.00144697 : 0.024108 > 6: -0.00178348 : 0.0240913 > .. > > I think that could be the type of the input data... because in agilent is > a waveform. And in my c++ application is real data.. > > In/out (FFTW) > data=(fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * (SIZE/2+1) ); > > ifft_result=(fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * (SIZE) > ); > > I am relative new in the Signal Processing... so i would need some help... > and if you know some books.. And excuse me for my English :D > > have a nice day...:cool: > >I suspect two things: First, everyone likes to scale their FFT and IFFT algorithms differently. Each one is the best -- just ask the author. So the fact that they're scaled differently is no big surprise. Check to see if they're off by 2048, 2048^2, or sqrt(2048), or any of those three with a factor of pi thrown in. Second, the fact that your Agilent program gives back a single column of numbers implies that it's a real FFT. If you're sticking your data into a complex FFT then you'd get back complex results, which isn't what you're looking for. See if taking these notions into consideration helps. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" gives you just what it says. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●July 18, 20082008-07-18
On Jul 17, 2:55 pm, "rapidiablo" <juan0...@hotmail.com> wrote:> i have this input data. (size = 2048) > > Agilent(VEE) c++ > > magnitude phase > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0.276 , @121.5 0.276 , 121.5 > 0.2061, @-72.03 0.2061, -72.03 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > 0 , @0 0 , 0 > ... ... ... ...To be clear, the above is the data you need to ifft, so you have a complex array in C++ that looks like data = 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 ... If the array is conjugate symmetric it will finish like this: ... 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.06358576 : 0.19604607 -0.14420960 : -0.23532869 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.00000000 : 0.00000000 With the conjugate symmetry shown above, a 2048-point ifft in scipy (a python package that uses either FFTPACK or FFTW) outputs the following: -7.87342198e-05 : -1.43836172e-12 -7.86143723e-05 : -1.44809333e-12 -7.84881940e-05 : -1.45767790e-12 ... -7.90531139e-05 : -1.40828912e-12 -7.89538571e-05 : -1.41845944e-12 -7.88474702e-05 : -1.42848388e-12 which agrees close enough with your VEE result for my money.> and this is the output from c++ (FFTW) > > 0: 0.000235391 : 0.0241332 > 1: -0.000100984 : 0.0241359 > 2: -0.000437429 : 0.0241348 > 3: -0.000773924 : 0.0241297 > 4: -0.00111045 : 0.0241208 > 5: -0.00144697 : 0.024108 > 6: -0.00178348 : 0.0240913 > .. > > I think that could be the type of the input data... because in agilent is > a waveform. And in my c++ application is real data.. > > In/out (FFTW) > data=(fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * (SIZE/2+1) ); > > ifft_result=(fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * (SIZE) > );Are you sure you stored your data properly in the array and are setting up your fftw_plan properly? Why is one array SIZE/2+1 in length while the other is SIZE? They should both have 2048 elements, no? (well, 4096 doubles, counting the real/imaginary components separately) A quick Google search pulled up this tutorial on FFTW: http://www-cecpv.u-strasbg.fr/Documentations/FFTW/v3.0.1/Complex-One-Dimensional-DFTs.html#Complex%20One-Dimensional%20DFTs By the way, you should be able to reinterpret_cast a C++ complex<double> * rather than use fftw_malloc. Both should just be a 2048x2 array of doubles. For example: const int SIZE = 2048; complex<double> *data = new complex<double>[SIZE]; complex<double> *ifft_result = new complex<double>[SIZE]; // ... fftw_plan p; p = fftw_plan_dft_1d( SIZE, reinterpret_cast<fftw_complex*>(data), reinterpret_cast<fftw_complex*>(ifft_result), FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute(p); // ... delete[] data; delete[] ifft_result;> I am relative new in the Signal Processing... so i would need some help... > and if you know some books.. And excuse me for my English :DThe Scientist and Engineer's Guide to DSP, a free online book, might be just what you need: http://www.dspguide.com/pdfbook.htm There's also the dspGuru site: http://www.dspguru.com -Eric
Reply by ●July 18, 20082008-07-18
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!!
Reply by ●July 18, 20082008-07-18
eric!! now i know how 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 You FFT them.. :) But.. i cant get the results form agilen... :( HELP!!!! thanxs!
Reply by ●July 18, 20082008-07-18
I forgot some info about agilent vee... This info about the input, could be my problem?? IN IFFT(agilent).txt: Type: Spectrum shape: Array 1D Size: 1025 mapping: LIN:0 ... 256.25 Thanks again and again..
Reply by ●July 18, 20082008-07-18
On Jul 18, 11:08�am, "rapidiablo" <juan0...@hotmail.com> wrote:> eric!! now i know how 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 > > You FFT them.. :) > > But.. i cant get the results form agilen... :( > > HELP!!!! > > thanxs!Those were your numbers. I simply took the amplitude/phase data you gave me and recreated the array so I could test it using the ifft in scipy. Now I have the precise numbers from IN IFFT(agilent).txt: fft_result = 0: 0.00000000 : 0.00000000 1: 0.00000000 : 0.00000000 2: 0.00000000 : 0.00000000 3: 0.00000000 : 0.00000000 4: 0.00000000 : 0.00000000 5: 0.00000000 : 0.00000000 6:-0.14404949: 0.23546151 7: 0.06356274 : -0.19600286 8: 0.00000000 : 0.00000000 ... 2040: 0.00000000 : 0.00000000 2041: 0.06356274 : 0.19600286 2042:-0.14404949: -0.23546151 2043: 0.00000000 : 0.00000000 2044: 0.00000000 : 0.00000000 2045: 0.00000000 : 0.00000000 2046: 0.00000000 : 0.00000000 2047: 0.00000000 : 0.00000000 ifft_result: -7.86003464e-05 : 0.00000000e+00 -7.87129910e-05 : -4.06575815e-20 -7.88070606e-05 : -1.35525272e-20 ... -7.81463489e-05 : -3.72112162e-19 -7.83174838e-05 : -3.86816124e-19 -7.84686641e-05 : -4.15979204e-19 For reference, these are the corresponding values from IFFT out (agilent).txt : -7.860034643685586E-005 -7.871299096183408E-005 -7.88070606286872E-005 ... -7.814634892479629E-005 -7.831748379365282E-005 -7.846866406201685E-005 Unsurprisingly, they agree, which means the problem is in the C code.
Reply by ●July 18, 20082008-07-18
ERIC!!! YOU ARE A GENIUS!!! I DON´T KNOW HOW TO THANKS YOU!!! ... The answer of my problem was in front of me !!! and i dint see it!! when you say: To be clear, the above is the data you need to ifft, so you have a complex array in C++ that looks like data = 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 .. If the array is conjugate symmetric it will finish like this: .. 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.06358576 : 0.19604607 -0.14420960 : -0.23532869 0.00000000 : 0.00000000 0.00000000 : 0.00000000 0.00000000 : 0.00000000 With the conjugate symmetry shown above, a 2048-point ifft in scipy (a python package that uses either FFTPACK or FFTW) outputs the following: .. my problem was that i just have an array of the half of the size and i never conjugate the values!!! so i have a 1025 arry in!! with the half of the data!!! in your last post!! i realized... what you told me in the first... (in the beginning i belive you just turn the values... for alien reason.. jeje) because you put the numbers of the elements (index). 0: 0.00000000 : 0.00000000 1: 0.00000000 : 0.00000000 2: 0.00000000 : 0.00000000 3: 0.00000000 : 0.00000000 4: 0.00000000 : 0.00000000 5: 0.00000000 : 0.00000000 6:-0.14404949: 0.23546151 7: 0.06356274 : -0.19600286 8: 0.00000000 : 0.00000000 =2E.. 2040: 0.00000000 : 0.00000000 2041: 0.06356274 : 0.19600286 2042:-0.14404949: -0.23546151 2043: 0.00000000 : 0.00000000 2044: 0.00000000 : 0.00000000 2045: 0.00000000 : 0.00000000 2046: 0.00000000 : 0.00000000 2047: 0.00000000 : 0.00000000 So... that was the problem!!! THANKS MEN!! ... Sorry about my English again... i do the best i can!! jeje.. Thxk every body!! :) xD. (im happy)
Reply by ●July 18, 20082008-07-18
On Jul 18, 10:47�am, "rapidiablo" <juan0...@hotmail.com> wrote:> � � vspan = asize/span; �/*Calcula el Stop del Star/Stop del agilent*/ > � � vmapi = mapi(vspan,asize); > � � dco (fft_result,fft_result2,vmapi);In the above code you're calculating the global variables hifrec and lowfrec, right?> � � alloc(fft_result,asize); > � � alloc(fft_result2,asize);Then you zero out everything but those two frequencies in the alloc function, right? But you only do that to elements 0 through 1025, leaving 1026 to 2047 unmodified. As an aside, I think your while loop exit condition in alloc() should have been (asize>n) rather than (asize>=n) since asize equals 1025. Element 0 is DC (0 Hz), and element 1024 is half the sampling frequency (digital frequency pi/2). For a real signal with 2048 sample points, elements 1025-2047 of the DFT (fft) are the conjugated mirror of elements 1-1023, so they don't provide any extra information, which is why they're often ignored.> � � fftw_execute( plan_backward ); �//Ejecuta el plan B > � � fftw_execute( plan_backward2 ); //Ejecuta el plan B.1Then you ifft your modified results. I think I see your problem. You need to zero out the conjugate symmetric frequencies in the 1025-2047 range, except the two at (SIZE - hifrec) and (SIZE - lowfrec). If you don't your result will not match VEE and will be complex valued. I hope this solves your problem. -Eric
Reply by ●July 18, 20082008-07-18
I already did that... everything works just fine!!! thanks eric!!.. you are a pro!! jeje!>On Jul 18, 10:47=A0am, "rapidiablo" <juan0...@hotmail.com> wrote: > >> =A0 =A0 vspan =3D asize/span; =A0/*Calcula el Stop del Star/Stop delagil=>ent*/ >> =A0 =A0 vmapi =3D mapi(vspan,asize); >> =A0 =A0 dco (fft_result,fft_result2,vmapi); > >In the above code you're calculating the global variables hifrec and >lowfrec, right? > >> =A0 =A0 alloc(fft_result,asize); >> =A0 =A0 alloc(fft_result2,asize); > >Then you zero out everything but those two frequencies in the alloc >function, right? But you only do that to elements 0 through 1025, >leaving 1026 to 2047 unmodified. > >As an aside, I think your while loop exit condition in alloc() should >have been (asize>n) rather than (asize>=3Dn) since asize equals 1025. >Element 0 is DC (0 Hz), and element 1024 is half the sampling >frequency (digital frequency pi/2). For a real signal with 2048 sample >points, elements 1025-2047 of the DFT (fft) are the conjugated mirror >of elements 1-1023, so they don't provide any extra information, which >is why they're often ignored. > >> =A0 =A0 fftw_execute( plan_backward ); =A0//Ejecuta el plan B >> =A0 =A0 fftw_execute( plan_backward2 ); //Ejecuta el plan B.1 > >Then you ifft your modified results. > >I think I see your problem. You need to zero out the conjugate >symmetric frequencies in the 1025-2047 range, except the two at (SIZE >- hifrec) and (SIZE - lowfrec). If you don't your result will not >match VEE and will be complex valued. I hope this solves your >problem. > >-Eric >
