CODEWITHSUNDEEP
c++dft
mr. Vachovsky
mr. Vachovsky

Reputation: 1128

DFT algorithm and convolution. what is wrong?




#include <vector>
using std::vector;

#include <complex>
using std::complex;
using std::polar;

typedef complex<double> Complex;

#define Pi 3.14159265358979323846

// direct Fourier transform
vector<Complex> dF( const vector<Complex>& in )
{
 const int N = in.size();

 vector<Complex> out( N );

 for (int k = 0; k < N; k++)
 {
  out[k] = Complex( 0.0, 0.0 );

  for (int n = 0; n < N; n++)
  {
   out[k] += in[n] * polar<double>( 1.0, - 2 * Pi * k * n / N );
  }
 }

 return out;
}

// inverse Fourier transform
vector<Complex> iF( const vector<Complex>& in )
{
 const int N = in.size();

 vector<Complex> out( N );

 for (int k = 0; k < N; k++)
 {
  out[k] = Complex( 0.0, 0.0 );

  for (int n = 0; n < N; n++)
  {
   out[k] += in[n] * polar<double>(1, 2 * Pi * k * n / N );
  }

  out[k] *= Complex(  1.0 / N , 0.0 );
 }

 return out;
}

Who can say, what is wrong???

Maybe i don't understand details of implementation this algorithm... But i can't find it )))

also, i need to calculate convolution.

But i can't find test example.

UPDATE

// convolution. I suppose that x0.size == x1.size
vector convolution( const vector& x0, const vector& x1) 
{
    const int N = x0.size();



vector<Complex> tmp( N );

for ( int i = 0; i < N; i++ )
{
    tmp[i] = x0[i] * x1[i];
}

return iF( tmp );




}

Upvotes: 1

Views: 2864

Answers (3)

lijie
lijie

Reputation: 4871

The transforms look fine, but there's nothing in the program that is doing convolution.

UPDATE: the convolution code needs to forward transform the inputs first before the element-wise multiplication.

Upvotes: 0

Jason B
Jason B

Reputation: 12975

I really don't know exactly what your asking, but your DFT and IDFT algorithms look correct to me. Convolution can be performed using the DFT and IDFT using the circular convolution theorem which basically states that f**g = IDFT(DFT(f) * DFT(g)) where ** is circular convolution and * is simple multiplication.

To compute linear convolution (non-circular) using the DFT, you must zero-pad each of the inputs so that the circular wrap-around only occurs for zero-valued samples and does not affect the output. Each input sequence needs to be zero padded to a length of N >= L+M-1 where L and M are the lengths of the input sequences. Then you perform circular convolution as shown above and the first L+M-1 samples are the linear convolution output (samples beyond this should be zero).

Note: Performing convolution with the DFT and IDFT algorithms you have shown is much more inefficient than just computing it directly. The advantage only comes when using an FFT and IFFT(O(NlogN)) algorithm in place of the DFT and IDFT (O(N^2)).

Upvotes: 2

Edward83
Edward83

Reputation: 6686

Check FFTW library "for computing the discrete Fourier transform (DFT)" and its C# wrapper;) Maybe this too;)

Good luck!

Upvotes: 1

Related Questions