I get different results using Kiss FFT and Octave/Matlab fft.
Here's the input/output of my code:
=========================================================
xtmp[2][4]
0 1 2 3
0 (0.2348,-0.5121) (-0.2745,0.0647)(-0.06092,-0.4704) (-0.04771,0.1112)
1 (-0.14,-0.04515)(-0.5469,0.4785) (-0.3687,-0.2009)(-0.04413,0.1298)
Input: Output:
0.234797 -0.512116 -0.1706950.026731
-0.274483 0.064698 0.128480-0.021765
-0.060917 -0.470399 0.0131070.008202
-0.047714 0.111247 0.012300-0.097780
-0.140031 -0.045153 0.0576210.186990
-0.546944 0.478488 0.0032580.104076
-0.368746 -0.200878 0.088650-0.140648
-0.044125 0.129769 -0.0153230.051592
x_tmp from ifft:
=========================================================
x_tmp[2][4]
0 1 2 3
0 (-0.1707,0.02673)(0.1285,-0.02177) (0.01311,0.008202)(0.0123,-0.09778)
1 (0.05762,0.187) (0.003258,0.1041)(0.08865,-0.1406)(-0.01532,0.05159)
-----------------------------------------------------------------------------
And, this is the correct output from Matlab/Octave:
x_tmp =
0.2348 - 0.5121i -0.2745 + 0.0647i -0.0609 - 0.4704i -0.0477 + 0.1112i
-0.1400 - 0.0452i -0.5469 + 0.4785i -0.3687 - 0.2009i -0.0441 + 0.1298i
x_tmp_after_ifft =
0.0474 - 0.2786i -0.4107 + 0.2716i -0.2148 - 0.3356i -0.0459 + 0.1205i
0.1874 - 0.2335i 0.1362 - 0.2069i 0.1539 - 0.1348i -0.0018 - 0.0093i
----------------------------------------------------------------------------
Here's my code calling fft from Kiss:
<code>
int i;
int isinverse = 1;
// TODO try row wise storage
int rows = x_tmp.size(); int cols=x_tmp[0].size();
double array_r[rows*cols];int idx=0;
double array_i[rows*cols];
double buf[rows*cols]; int size=rows*cols;
double buf_i[rows*cols];
for(int i=0; i<rows;i++)
for(int j=0;j<cols;j++)
{
array_r[idx]=x_tmp[i][j].real();
array_i[idx]=x_tmp[i][j].imag();
idx++;
}
kiss_fft_cpx out_cpx[size],out[size],*cpx_buf;
kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,0,0);
kiss_fftr_cfg ifft = kiss_fftr_alloc(size*2,isinverse,0,0);
cpx_buf = copycpx(array_r, array_i, size);
// kiss_fftr(fft,(kiss_fft_scalar*)cpx_buf, out_cpx);
// out_cpx has our input
//kiss_fftri(ifft,out_cpx,(kiss_fft_scalar*)out );
kiss_fftri(ifft,cpx_buf,(kiss_fft_scalar*)out );
printf("Input:\t\tOutput:\n");
for(i=0;i<size;i++)
{
buf[i] = (out[i].r)/(size*2);
buf_i[i]=(out[i].i)/(size*2);
printf("%f%f\t\t%f%f\n",array_r[i]
, array_i[i],buf[i], buf_i[i]);
}
kiss_fft_cleanup();
free(fft);
free(ifft);
</code>
------------------------------------------------------
Sorry about the formatting.
If Mark Borgerding is around, that would be great.
Thank you very much. I don't know where is my mistake.
Problem Using Kiss FFT
Started by ●July 20, 2014
Reply by ●July 20, 20142014-07-20
On Sunday, July 20, 2014 2:06:34 AM UTC-4, Javier Pajuelo wrote:> I get different results using Kiss FFT and Octave/Matlab fft. > > > > Here's the input/output of my code: > > ========================================================= > > xtmp[2][4] > > 0 1 2 3 > > 0 (0.2348,-0.5121) (-0.2745,0.0647)(-0.06092,-0.4704) (-0.04771,0.1112) > > 1 (-0.14,-0.04515)(-0.5469,0.4785) (-0.3687,-0.2009)(-0.04413,0.1298) > > > > Input: Output: > > 0.234797 -0.512116 -0.1706950.026731 > > -0.274483 0.064698 0.128480-0.021765 > > -0.060917 -0.470399 0.0131070.008202 > > -0.047714 0.111247 0.012300-0.097780 > > -0.140031 -0.045153 0.0576210.186990 > > -0.546944 0.478488 0.0032580.104076 > > -0.368746 -0.200878 0.088650-0.140648 > > -0.044125 0.129769 -0.0153230.051592 > > > > > > x_tmp from ifft: > > > > ========================================================= > > x_tmp[2][4] > > 0 1 2 3 > > 0 (-0.1707,0.02673)(0.1285,-0.02177) (0.01311,0.008202)(0.0123,-0.09778) > > 1 (0.05762,0.187) (0.003258,0.1041)(0.08865,-0.1406)(-0.01532,0.05159) > > > > > > > > ----------------------------------------------------------------------------- > > And, this is the correct output from Matlab/Octave: > > > > x_tmp = > > > > 0.2348 - 0.5121i -0.2745 + 0.0647i -0.0609 - 0.4704i -0.0477 + 0.1112i > > -0.1400 - 0.0452i -0.5469 + 0.4785i -0.3687 - 0.2009i -0.0441 + 0.1298i > > > > > > x_tmp_after_ifft = > > > > 0.0474 - 0.2786i -0.4107 + 0.2716i -0.2148 - 0.3356i -0.0459 + 0.1205i > > 0.1874 - 0.2335i 0.1362 - 0.2069i 0.1539 - 0.1348i -0.0018 - 0.0093i > > > > ---------------------------------------------------------------------------- > > Here's my code calling fft from Kiss: > > > > <code> > > int i; > > int isinverse = 1; > > > > // TODO try row wise storage > > int rows = x_tmp.size(); int cols=x_tmp[0].size(); > > double array_r[rows*cols];int idx=0; > > double array_i[rows*cols]; > > double buf[rows*cols]; int size=rows*cols; > > double buf_i[rows*cols]; > > > > > > for(int i=0; i<rows;i++) > > for(int j=0;j<cols;j++) > > { > > array_r[idx]=x_tmp[i][j].real(); > > array_i[idx]=x_tmp[i][j].imag(); > > idx++; > > } > > > > kiss_fft_cpx out_cpx[size],out[size],*cpx_buf; > > > > kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,0,0); > > kiss_fftr_cfg ifft = kiss_fftr_alloc(size*2,isinverse,0,0); > > > > cpx_buf = copycpx(array_r, array_i, size); > > > > // kiss_fftr(fft,(kiss_fft_scalar*)cpx_buf, out_cpx); > > > > // out_cpx has our input > > //kiss_fftri(ifft,out_cpx,(kiss_fft_scalar*)out ); > > kiss_fftri(ifft,cpx_buf,(kiss_fft_scalar*)out ); > > > > printf("Input:\t\tOutput:\n"); > > for(i=0;i<size;i++) > > { > > buf[i] = (out[i].r)/(size*2); > > buf_i[i]=(out[i].i)/(size*2); > > > > printf("%f%f\t\t%f%f\n",array_r[i] > > , array_i[i],buf[i], buf_i[i]); > > } > > > > kiss_fft_cleanup(); > > free(fft); > > free(ifft); > > > > </code> > > ------------------------------------------------------ > > > > Sorry about the formatting. > > If Mark Borgerding is around, that would be great. > > > > Thank you very much. I don't know where is my mistake.Just a thought. Why not use FFTW instead? Dirk
Reply by ●July 20, 20142014-07-20
> > Just a thought. Why not use FFTW instead? > > > > DirkYeah, I will try it as my last resort. The problem is that it is too heavy and I don't know how to use it. I am still looking for help using the Kiss library, so if someone recognizes the problem please let me know.
Reply by ●July 20, 20142014-07-20
Javier Pajuelo wrote:> >> >> Just a thought. Why not use FFTW instead? >> >> >> >> Dirk > > Yeah, I will try it as my last resort. The problem is that it is too heavy and I don't know how to use it. >I found it easier to use by wrapping it with C++. That's usually a waste of time but with FFTW, not so much. It helps with managing all the state.> I am still looking for help using the Kiss library, so if someone recognizes the problem please let me know. >-- Les Cargill
Reply by ●July 20, 20142014-07-20
Ok, I managed to call FFTW3, but I still get different results.
I am calling the ifft directly from some row major complex vector.
TEST03
Demonstrate FFTW3 on a single vector of complex data.
Transform data to FFT coefficients.
Backtransform FFT coefficients to recover data.
Compare recovered data to original data.
Input Data:
0 0.234797 -0.512116
1 -0.274483 0.064698
2 -0.060917 -0.470399
3 -0.047714 0.111247
4 -0.140031 -0.045153
5 -0.546944 0.478488
6 -0.368746 -0.200878
7 -0.044125 0.129769
Recovered input data:
0 -1.248163 -0.444346
1 1.145235 -0.248510
2 0.222258 -0.615579
3 0.215801 -0.305174
4 0.578368 -2.012748
5 0.143462 -0.069761
6 0.826599 0.843596
7 -0.005187 -1.244409
Recovered input data divided by N:
0 -0.156020 -0.055543
1 0.143154 -0.031064
2 0.027782 -0.076947
3 0.026975 -0.038147
4 0.072296 -0.251594
5 0.017933 -0.008720
6 0.103325 0.105449
7 -0.000648 -0.155551
x_tmp from ifft:
x_tmp[2][4]
0 1 2 3
0 (-0.156,-0.05554) (0.1432,-0.03106) (0.02778,-0.07695) (0.02698,-0.03815)
1 (0.0723,-0.2516) (0.01793,-0.00872) (0.1033,0.1054) (-0.0006484,-0.1556)
-------------
Code is:
04 November 2007
Author:
John Burkardt
*/
{
int i;
fftw_complex *in;
fftw_complex *in2;
int n = 100;
fftw_complex *out;
fftw_plan plan_backward;
fftw_plan plan_forward;
unsigned int seed = 123456789;
printf ( "\n" );
printf ( "TEST01\n" );
printf ( " Demonstrate FFTW3 on a single vector of complex data.\n" );
printf ( "\n" );
printf ( " Transform data to FFT coefficients.\n" );
printf ( " Backtransform FFT coefficients to recover data.\n" );
printf ( " Compare recovered data to original data.\n" );
/*
Create the input array.
*/
in = fftw_malloc ( sizeof ( fftw_complex ) * n );
// srand ( seed );
// filled row major, n is rows*cols
for ( i = 0; i < n; i++ )
{
in[i][0] = frand ( ); // I filled this with my real info
in[i][1] = frand ( ); // same but for complex
}
printf ( "\n" );
printf ( " Input Data:\n" );
printf ( "\n" );
for ( i = 0; i < n; i++ )
{
printf ( " %3d %12f %12f\n", i, in[i][0], in[i][1] );
}
/*
Create the output array.
*/
//out = fftw_malloc ( sizeof ( fftw_complex ) * n );
//plan_forward = fftw_plan_dft_1d ( n, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
/* fftw_execute ( plan_forward );
printf ( "\n" );
printf ( " Output FFT Coefficients:\n" );
printf ( "\n" );
for ( i = 0; i < n; i++ )
{
printf ( " %3d %12f %12f\n", i, out[i][0], out[i][1] );
}*/
/*
Recreate the input array.
*/
in2 = fftw_malloc ( sizeof ( fftw_complex ) * n );
/* notice in is matrix, I want ifft */
plan_backward = fftw_plan_dft_1d ( n, in, in2, FFTW_BACKWARD, FFTW_ESTIMATE );
fftw_execute ( plan_backward );
printf ( "\n" );
printf ( " Recovered input data:\n" );
printf ( "\n" );
for ( i = 0; i < n; i++ )
{
printf ( " %3d %12f %12f\n", i, in2[i][0], in2[i][1] );
}
printf ( "\n" );
printf ( " Recovered input data divided by N:\n" );
printf ( "\n" );
for ( i = 0; i < n; i++ )
{
printf ( " %3d %12f %12f\n", i,
in2[i][0] / ( double ) ( n ), in2[i][1] / ( double ) ( n ) );
}
/*
Free up the allocated memory.
*/
fftw_destroy_plan ( plan_forward );
fftw_destroy_plan ( plan_backward );
fftw_free ( in );
fftw_free ( in2 );
fftw_free ( out );
I assume it is correct. It's described here: http://people.sc.fsu.edu/~jburkardt/c_src/fftw3/fftw3_prb.c
Reply by ●July 20, 20142014-07-20
According to FFTW3's website, the transforms are not normalized on the other hand matlab/octave compute normalized transforms. How do I normalize my matrix to get same outputs as Matlab/Octave?
Reply by ●July 20, 20142014-07-20
Javier Pajuelo wrote:> According to FFTW3's website, the transforms are not normalized on > the other hand matlab/octave compute normalized transforms. > > How do I normalize my matrix to get same outputs as Matlab/Octave? >Divide the vector of results by the number of buckets. -- Les Cargill
Reply by ●July 20, 20142014-07-20
Thank you all. Besides that normalization issue. The problem lied in definition of fft in matlab/Octave. 'fft and ifft' do the FFTs in each column independently. I was calling a FTW3 2d FFT, which was incorrect. Now, I get similar results. Thanks again.
Reply by ●July 21, 20142014-07-21
Just a bit OT, but may somewhat fit in here... I noticed that when I do a LOT of FFTs in scilab 5.5, it keeps allocating memory which is never freed. Scilab uses FFTW as far as I know... Has anybody knows any cure for this? simple sample code: while(1) a = rand(1,1024); b = fft(a); end I came across this when using FFTs heavily for iteratively optimizing FIR filters. Scilab will crash after some 100k loops. best regards, Andre>> >> >> Thank you very much. I don't know where is my mistake. > > Just a thought. Why not use FFTW instead? > > Dirk >
Reply by ●July 21, 20142014-07-21
> On Sunday, July 20, 2014 2:06:34 AM UTC-4, Javier Pajuelo wrote: > I get different results using Kiss FFT and Octave/Matlab fft. > Here's the input/output of my code: > Thank you very much. I don't know where is my mistake.The simplest thing to do is run each set of software through a 2-way translation test A -> B -> A' involving the "A -> B" examples in question and find out which ones are most closely matching A' = A. On Sunday, July 20, 2014 1:31:23 AM UTC-5, belld...@gmail.com wrote:> Just a thought. Why not use FFTW instead?I once wrote a parody of the legend of Phaeton that had the throwaway line "He climbed one of the tall mountains and found a palace on its summit which he immediately knew belonged to Apollo. 'Sheesh', Phaeton thought to himself, 'I knew it would be tough getting here, but who'd have ever thought trying to reach Apollo would become such a huge NASA mission?' I think that's the gist of why not to use FFTW, if you don't have to. I found an interesting alternative a short while back ... and one that's almost as old as Phaeton: AlgLib. It used to be in FORTRAN but was moved over (by automated translation) to C++. Although much of it is buried underneath the author's attempt at hand-coding a multi-threaded runtime system, the code for the key routines is easy to extract and use elsewhere. And it's not too difficult to smoothe over into C if you're not using C++. In particular, the Fast Transforms section: http://www.alglib.net/fasttransforms/ has fast routines for convolution and DFT that are NOT limited to powers of 2 sizes. Small prime factors in the DFT size N are done with codelets (like FFTW or its Fortran ancestor FFTPACK), large prime factors are done with the Bluestein algorithm. The FFT parts of the .cpp and .h files are about 1000 lines each only because the entire manual is practically inserted in the code (a habit inherited from the FORTRAN days). Separate out the manual into its own reference and the code will be small and easy to work with. (And actually: AlgLib already has a web page that contains a copy of the manual).
