Why doesn't my Inverse DFT work correctly?

I'm hand waiving over a lot of theory here, but essentially it boils down to: your input function has an infinite frequency, so you'd need to use an infinite number of frequency samples to perfectly reconstruct it. Because your maximum frequency of your input is greater than the Nyquist frequency (look it up if you don't know it), you lose some information about the original signal.

I plotted your result and it looks like what I'd expect to see given the parameters you've used. Original signal is in red, the result is in blue.

[attachment=2323:dft.jpg]

There's nothing wrong with your code, you just need to brush up on theory a little bit.

I'm hand waiving over a lot of theory here, but essentially it boils down to: your input function has an infinite frequency, so you'd need to use an infinite number of frequency samples to perfectly reconstruct it. Because your maximum frequency of your input is greater than the Nyquist frequency (look it up if you don't know it), you lose some information about the original signal.

That is only true if you sample the signal. His signal is already sampled, so any data loss is already present in the input signal. The forward and backward Fourier transform is a unique one-to-one mapping, meaning the inverse FFT will exactly reconstruct the input.

From Matlab using his input sequence.

>> x = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
>> X = fft(x);
>> y = ifft(X);
y =
Columns 1 through 11
0 0 0 -0.0000 0 0 0 0 0 0 0
Columns 12 through 22
-0.0000 0 0 0 1.0000 0 0 0 -0.0000 0 0
Columns 23 through 32
0 0 0 0 0 -0.0000 0 0 0 0

The sequence is, subject to rounding errors, exactly reconstructed.

I'm just testing something, so no need to lecture me on the algorithm please, thank you. I'm running a DFT on an impulse function. Then printing out the data, then running inverse DFT on the frequency domain values.
I'm not sure why the inverse DFT doesn't return back the impulse function. The algorithm is based straight from the definition. Here is the code :


#include <iostream>
#include <ctime>
#include <set>
#include <algorithm>
#include <vector>
#include <cmath>

using namespace std;

typedef float Type;
typedef std::vector<Type> Vector;
const Type TWO_PI = Type(3.141519 * 2.0);

struct FrequencyDomain{
Vector realPart;
Vector imaginaryPart;
};
FrequencyDomain calculateDFT(const Vector& sampleData){
FrequencyDomain freqDomain;

for(unsigned k = 0; k < sampleData.size()/2; ++k){
Type realVal = 0;
Type imVal = 0;
for(unsigned i = 0; i < sampleData.size(); ++i){
float theta = TWO_PI * k * i / float(sampleData.size());
realVal += sampleData * cos(theta);
imVal += sampleData * sin( theta );
}
freqDomain.realPart.push_back(realVal);
freqDomain.imaginaryPart.push_back(-imVal);
}
return freqDomain;
}
Vector calculateInverseDFT(const FrequencyDomain& d){
Vector v;
const int SAMPLE_SIZE = d.realPart.size() * 2;
FrequencyDomain f;
//adjust proper sizing
f.realPart = Vector(d.realPart.size());
f.imaginaryPart = Vector(d.imaginaryPart.size());

//calculate the special case
f.realPart[0] = d.realPart[0] / SAMPLE_SIZE;
f.realPart[f.realPart.size()-1] = d.realPart.back() / SAMPLE_SIZE;

//convert to cos amplitudes
for(unsigned i = 1; i < d.realPart.size() - 1; ++i){
f.realPart = d.realPart / (SAMPLE_SIZE/2.0f);
}
//convert to sin amplitudes
for(unsigned i = 0; i < d.imaginaryPart.size(); ++i){
f.imaginaryPart = -d.imaginaryPart / (SAMPLE_SIZE/2.0f);
}

//run calculation
for(int i = 0; i < SAMPLE_SIZE; ++i){
Type realVal = 0;
Type imVal = 0;
for(int k = 0; k < SAMPLE_SIZE/2; ++k){
float theta = TWO_PI * k * i/ float(SAMPLE_SIZE);
realVal += f.realPart[k] * cos( theta );
imVal += f.imaginaryPart[k] * sin( theta );
}
v.push_back(realVal + imVal );
}
return v;
}
template<int MIN, int MAX>
struct RandRange{
int operator()(){ return rand() % (MAX-MIN) + MIN; }
};
template<typename Itr>
void print(Itr begin, Itr end, const std::string& delim = " "){
cout << "[ ";
while(begin != end){ cout << *begin << delim; ++begin; }
cout << "]";
}
int main() {
srand( time(0) );
const int GEN_SIZE = 32;
Vector sampleData(GEN_SIZE,0);
sampleData[15] = 1;
//std::generate_n(sampleData.begin(),GEN_SIZE,RandRange<1,200>());

cout << "Sampled data: ";
print( sampleData.begin(), sampleData.end());
cout << endl;

FrequencyDomain f = calculateDFT(sampleData);
cout << "Real Part: ";
print(f.realPart.begin(), f.realPart.end());
cout << "\nImaginary Part: ";
print(f.imaginaryPart.begin(),f.imaginaryPart.end());
cout << "\n\n";

Vector v = calculateInverseDFT(f);
cout << "Time Domain after inverse: ";
print(v.begin(),v.end());
cout << "\n";

return 0;
}



here is the output:

Sampled data: [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Real Part: [ 1 -0.980772 0.923827 -0.831354 0.706912 -0.555284 0.3823 -0.194616 -0.000551164 0.195699 -0.38332 0.556201 -0.707694 0.831966 -0.924249 0.980987 ]

Imaginary Part: [ -0 -0.195158 0.382811 -0.555743 0.707302 -0.831661 0.924038 -0.980879 1 -0.980664 0.923616 -0.831048 0.706519 -0.554826 0.38179 -0.194076 ]




Time Domain after inverse: [ 0.000597399 -0.00118974 0.0029381 -0.00577524 0.00959269 -0.014243 0.0195472 -0.0253024 0.0312866 -0.0372698 0.0430225 -0.0483229 0.0529676 -0.0567778 0.0596068 0.938677 0.0619292 -0.061333 0.0595804 -0.0567395 0.0529184 -0.0482652 0.042958 -0.0372013 0.0312163 -0.0252324 0.0194809 -0.014182 0.00954051 -0.00573279 0.00290673 -0.0011705 ]




can anyone see what i'm doing wrong? BTW, i'm reading from this site. Thanks.

regards, Dibash Chhetri

Not only is the inverse FFT wrong, but the forward FFT is possibly wrong as well.

You're only calculating half the FFT, since the output sequence is 16 samples but the input is 32. The FFT for real input data is indeed complex conjugate symmetric and half the FFT is indeed redundant in itselt. But that doesn't mean you can just discard it and go ahead and calculate the IFFT the usual way without compensating for it. It is complex conjugate symmetric, but it is there for a reason. There are two usual ways to handle the symmetry.

1: Don't care about the symmetry and calculate the full FFT.
2: Calculate half the FFT, but then you also have to modify the IFFT for this. Since all frequency bins except first bin and and the Nyquist bin are effectively present twice due to the symmetry, you need to sum those twice. The first bin and the Nyquist bin are present only once, so you need to sum them only once.

But let's say you're going for option 1 and did indeed calculate the full 32 complex conjugate symmetric set of frequency bins, your IFFT is wrong. Let's focus on the inner loop.
[source]
for(int i = 0; i < SAMPLE_SIZE; ++i){
Type realVal = 0;
Type imVal = 0;
for(int k = 0; k < SAMPLE_SIZE/2; ++k){
float theta = TWO_PI * k * i/ float(SAMPLE_SIZE);
realVal += f.realPart[k] * cos( theta );
imVal += f.imaginaryPart[k] * sin( theta );
}
v.push_back(realVal + imVal );
}
[/source]
First, k should loop from 0 to SAMPLE_SIZE-1, since you need the full FFT, not half of it. But that was what I assumed above. But most importantly, your complex multiplication is wrong. The cos and sin in your code are not just cos and sin factors for the real and imaginary part; they are the two components of the complex exponential correlator. Likewise, the real and imaginary components of the signal are the two parts of a complex value. Thus, you need to perform proper complex multiplication.
[source]
realVal += f.realPart[k] * cos( theta ) + f.imaginaryPart[k] * sin( theta );
imVal += f.imaginaryPart[k] * cos( theta ) - f.realPart[k] * sin( theta );
[/source]
Then, when you push the final calculate value to the vector, the value is not realVal+imVal, but the complex value formed by the two variables representing the real and the imaginary part. You should technically push the complex value (realVal, imVal), but due to the complex conjugate symmetry of the FFT, imVal should be very close to zero since the output must be real valued. You can, effectively, assume that imVal is zero and push only realVal (which is what you end up doing, but conceptually it's wrong to push the sum of the two, since it suggests something you shouldn't do if it's properly implemented).

I'm hand waiving over a lot of theory here, but essentially it boils down to: your input function has an infinite frequency, so you'd need to use an infinite number of frequency samples to perfectly reconstruct it. Because your maximum frequency of your input is greater than the Nyquist frequency (look it up if you don't know it), you lose some information about the original signal.

This is highly incorrect: the input is not a continuous signal (which, indeed, could possibly exceed the Nyquist limit), but 32 samples. Whatever (complex) value the 32 samples have, the DFT transforms them into 32 frequency-domain samples, and the inverse DFT is supposed to reconstruct the input exactly (up to small differences due to floating point roundoff, underflow, etc.); errors exceeding 0.05 are bad enough to investigate.

This is highly incorrect: the input is not a continuous signal (which, indeed, could possibly exceed the Nyquist limit), but 32 samples. Whatever (complex) value the 32 samples have, the DFT transforms them into 32 frequency-domain samples, and the inverse DFT is supposed to reconstruct the input exactly (up to small differences due to floating point roundoff, underflow, etc.); errors exceeding 0.05 are bad enough to investigate.

Thank you, Lorenzo. The DFT is just a big rotation matrix. [Well, almost. Depending on the sample count (even/odd) there can also be a reflection in there, so really it's just an isometry. And it's only really an isometry if you choose the "symmetric" scaling factors for the DFT/IFT pair. But modulo those details, it's pretty much a rotation. ;-) ]

I'm not sure what you did wrong (other than using a pretty bad value for pi, but that doesn't explain the whole problem).

This is a totally naive implementation that does work, so perhaps you can compare intermediate results:
#include <iostream>
#include <complex>

typedef double R;
typedef std::complex<R> C;
C I(0,1);
R TWO_PI = std::atan(1)*8;

int main() {
R original[32]={0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,1,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0};

std::cout << "original: ";
for (int n=0; n<32; ++n)
std::cout << original[n] << ' ';
std::cout << '\n';

C dft[32];
for (int n=0; n<32; ++n) {
C s = 0;
for (int m=0; m<32; ++m)
s += std::exp(-TWO_PI*n*m/R(32)*I)*original[m];
dft[n] = s;
}

std::cout << "dft: ";
for (int n=0; n<32; ++n)
std::cout << dft[n] << ' ';
std::cout << '\n';

C reconstruction[32];
for (int n=0; n<32; ++n) {
C s = 0;
for (int m=0; m<32; ++m)
s += std::exp(-TWO_PI*n*m/R(32)*I)*dft[m];
reconstruction[n] = s/R(32);
}

std::cout << "reconstruction: ";
for (int n=0; n<32; ++n)
std::cout << reconstruction[n] << ' ';
std::cout << '\n';
}



[quote name='karwosts' timestamp='1305608072' post='4811775']
I'm hand waiving over a lot of theory here, but essentially it boils down to: your input function has an infinite frequency, so you'd need to use an infinite number of frequency samples to perfectly reconstruct it. Because your maximum frequency of your input is greater than the Nyquist frequency (look it up if you don't know it), you lose some information about the original signal.

That is only true if you sample the signal. His signal is already sampled, so any data loss is already present in the input signal. The forward and backward Fourier transform is a unique one-to-one mapping, meaning the inverse FFT will exactly reconstruct the input.

From Matlab using his input sequence.

>> x = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
>> X = fft(x);
>> y = ifft(X);
y =
Columns 1 through 11
0 0 0 -0.0000 0 0 0 0 0 0 0
Columns 12 through 22
-0.0000 0 0 0 1.0000 0 0 0 -0.0000 0 0
Columns 23 through 32
0 0 0 0 0 -0.0000 0 0 0 0

The sequence is, subject to rounding errors, exactly reconstructed.

I'm just testing something, so no need to lecture me on the algorithm please, thank you. I'm running a DFT on an impulse function. Then printing out the data, then running inverse DFT on the frequency domain values.
I'm not sure why the inverse DFT doesn't return back the impulse function. The algorithm is based straight from the definition. Here is the code :


#include <iostream>
#include <ctime>
#include <set>
#include <algorithm>
#include <vector>
#include <cmath>

using namespace std;

typedef float Type;
typedef std::vector<Type> Vector;
const Type TWO_PI = Type(3.141519 * 2.0);

struct FrequencyDomain{
Vector realPart;
Vector imaginaryPart;
};
FrequencyDomain calculateDFT(const Vector& sampleData){
FrequencyDomain freqDomain;

for(unsigned k = 0; k < sampleData.size()/2; ++k){
Type realVal = 0;
Type imVal = 0;
for(unsigned i = 0; i < sampleData.size(); ++i){
float theta = TWO_PI * k * i / float(sampleData.size());
realVal += sampleData * cos(theta);
imVal += sampleData * sin( theta );
}
freqDomain.realPart.push_back(realVal);
freqDomain.imaginaryPart.push_back(-imVal);
}
return freqDomain;
}
Vector calculateInverseDFT(const FrequencyDomain& d){
Vector v;
const int SAMPLE_SIZE = d.realPart.size() * 2;
FrequencyDomain f;
//adjust proper sizing
f.realPart = Vector(d.realPart.size());
f.imaginaryPart = Vector(d.imaginaryPart.size());

//calculate the special case
f.realPart[0] = d.realPart[0] / SAMPLE_SIZE;
f.realPart[f.realPart.size()-1] = d.realPart.back() / SAMPLE_SIZE;

//convert to cos amplitudes
for(unsigned i = 1; i < d.realPart.size() - 1; ++i){
f.realPart = d.realPart / (SAMPLE_SIZE/2.0f);
}
//convert to sin amplitudes
for(unsigned i = 0; i < d.imaginaryPart.size(); ++i){
f.imaginaryPart = -d.imaginaryPart / (SAMPLE_SIZE/2.0f);
}

//run calculation
for(int i = 0; i < SAMPLE_SIZE; ++i){
Type realVal = 0;
Type imVal = 0;
for(int k = 0; k < SAMPLE_SIZE/2; ++k){
float theta = TWO_PI * k * i/ float(SAMPLE_SIZE);
realVal += f.realPart[k] * cos( theta );
imVal += f.imaginaryPart[k] * sin( theta );
}
v.push_back(realVal + imVal );
}
return v;
}
template<int MIN, int MAX>
struct RandRange{
int operator()(){ return rand() % (MAX-MIN) + MIN; }
};
template<typename Itr>
void print(Itr begin, Itr end, const std::string& delim = " "){
cout << "[ ";
while(begin != end){ cout << *begin << delim; ++begin; }
cout << "]";
}
int main() {
srand( time(0) );
const int GEN_SIZE = 32;
Vector sampleData(GEN_SIZE,0);
sampleData[15] = 1;
//std::generate_n(sampleData.begin(),GEN_SIZE,RandRange<1,200>());

cout << "Sampled data: ";
print( sampleData.begin(), sampleData.end());
cout << endl;

FrequencyDomain f = calculateDFT(sampleData);
cout << "Real Part: ";
print(f.realPart.begin(), f.realPart.end());
cout << "\nImaginary Part: ";
print(f.imaginaryPart.begin(),f.imaginaryPart.end());
cout << "\n\n";

Vector v = calculateInverseDFT(f);
cout << "Time Domain after inverse: ";
print(v.begin(),v.end());
cout << "\n";

return 0;
}



here is the output:

Sampled data: [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]

Real Part: [ 1 -0.980772 0.923827 -0.831354 0.706912 -0.555284 0.3823 -0.194616 -0.000551164 0.195699 -0.38332 0.556201 -0.707694 0.831966 -0.924249 0.980987 ]

Imaginary Part: [ -0 -0.195158 0.382811 -0.555743 0.707302 -0.831661 0.924038 -0.980879 1 -0.980664 0.923616 -0.831048 0.706519 -0.554826 0.38179 -0.194076 ]




Time Domain after inverse: [ 0.000597399 -0.00118974 0.0029381 -0.00577524 0.00959269 -0.014243 0.0195472 -0.0253024 0.0312866 -0.0372698 0.0430225 -0.0483229 0.0529676 -0.0567778 0.0596068 0.938677 0.0619292 -0.061333 0.0595804 -0.0567395 0.0529184 -0.0482652 0.042958 -0.0372013 0.0312163 -0.0252324 0.0194809 -0.014182 0.00954051 -0.00573279 0.00290673 -0.0011705 ]




can anyone see what i'm doing wrong? BTW, i'm reading from this site. Thanks.

regards, Dibash Chhetri

Not only is the inverse FFT wrong, but the forward FFT is possibly wrong as well.

You're only calculating half the FFT, since the output sequence is 16 samples but the input is 32. The FFT for real input data is indeed complex conjugate symmetric and half the FFT is indeed redundant in itselt. But that doesn't mean you can just discard it and go ahead and calculate the IFFT the usual way without compensating for it. It is complex conjugate symmetric, but it is there for a reason. There are two usual ways to handle the symmetry.

1: Don't care about the symmetry and calculate the full FFT.
2: Calculate half the FFT, but then you also have to modify the IFFT for this. Since all frequency bins except first bin and and the Nyquist bin are effectively present twice due to the symmetry, you need to sum those twice. The first bin and the Nyquist bin are present only once, so you need to sum them only once.

But let's say you're going for option 1 and did indeed calculate the full 32 complex conjugate symmetric set of frequency bins, your IFFT is wrong. Let's focus on the inner loop.
[source]
for(int i = 0; i < SAMPLE_SIZE; ++i){
Type realVal = 0;
Type imVal = 0;
for(int k = 0; k < SAMPLE_SIZE/2; ++k){
float theta = TWO_PI * k * i/ float(SAMPLE_SIZE);
realVal += f.realPart[k] * cos( theta );
imVal += f.imaginaryPart[k] * sin( theta );
}
v.push_back(realVal + imVal );
}
[/source]
First, k should loop from 0 to SAMPLE_SIZE-1, since you need the full FFT, not half of it. But that was what I assumed above. But most importantly, your complex multiplication is wrong. The cos and sin in your code are not just cos and sin factors for the real and imaginary part; they are the two components of the complex exponential correlator. Likewise, the real and imaginary components of the signal are the two parts of a complex value. Thus, you need to perform proper complex multiplication.
[source]
realVal += f.realPart[k] * cos( theta ) + f.imaginaryPart[k] * sin( theta );
imVal += f.imaginaryPart[k] * cos( theta ) - f.realPart[k] * sin( theta );
[/source]
Then, when you push the final calculate value to the vector, the value is not realVal+imVal, but the complex value formed by the two variables representing the real and the imaginary part. You should technically push the complex value (realVal, imVal), but due to the complex conjugate symmetry of the FFT, imVal should be very close to zero since the output must be real valued. You can, effectively, assume that imVal is zero and push only realVal (which is what you end up doing, but conceptually it's wrong to push the sum of the two, since it suggests something you shouldn't do if it's properly implemented).
[/quote]

Thanks for the advice. Taking your advice, I changed the FFT so that the FreqDomain values range from 0-SAMPLE_SIZE - 1 instead of 0 - SAMPLE_SIZE/2. Then made the proper changes in the IFFT, that is change the range of k to 0-SAMPLE_SIZE - 1, and change to the proper complex multiplication.

The results I got from this was:
[source]
Sampled data: [ 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,]

Real Part: [ 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 0.999999 -0.999999 ]

Imaginary Part: [ -0 -7.35839e-05 0.000147168 -0.000221229 0.000294336 -0.000367443 0.000442457 -0.000515564 0.000588671 -0.000661778 0.000734885 -0.0008099 0.000884914 -0.000956114 0.00103113 -0.00110233 ]




Time Domain after inverse: [ 1.86265e-08,-6.28084e-06,-1.12057e-05,-1.53407e-05,-1.84365e-05,-2.05617e-05,-2.22996e-05,-2.31415e-05,1,-2.31378e-05,-2.22139e-05,-2.04779e-05,-1.81757e-05,-1.52197e-05,-1.14515e-05,-5.96419e-06,]


[/source]
which looks very good, and after transforming the data by applying a rounding algorithm, it has perfect integer reconstruction.

Can you explain how you got the proper complex correlator? I understand using euler equation I could reform DFT into cos and sin, but I don't understand why it isn't

realVal += f.realPart[k] * cos( theta ) + f.imaginaryPart[k] * sin( theta );
imVal += f.realPart[k] * cos( theta ) - f.imaginaryPart[k] * sin( theta );

because complex conjugate essentially means we change the axis 'i' lies in, right?

Can you explain how you got the proper complex correlator? I understand using euler equation I could reform DFT into cos and sin, but I don't understand why it isn't

realVal += f.realPart[k] * cos( theta ) + f.imaginaryPart[k] * sin( theta );
imVal += f.realPart[k] * cos( theta ) - f.imaginaryPart[k] * sin( theta );

because complex conjugate essentially means we change the axis 'i' lies in, right?

Take a look at how complex multiplication is performed; here for example. You want to multiply the complex number (f.realPart[k] + i*f.imaginaryPart[k]) with the complex number (cos(theta + i*sin(theta)). The first is your FFT output, and the second is the complex exponential correlator.

However, the formula I gave you had some signed swapped, but it turned out correct because your theta was also incorrect as well, and they happen to cancel eachothers effect. The correct formula is
[source]
realVal += f.realPart[k] * cos( theta ) - f.imaginaryPart[k] * sin( theta );
imVal += f.imaginaryPart[k] * cos( theta ) + f.realPart[k] * sin( theta );
[/source]
and
[source]
float theta = -TWO_PI * k * i/ float(SAMPLE_SIZE);
[/source]
Note the sign change in both codes. The calculation of theta with a minus sign is for the inverse DFT only. The frequencies of the correlators in the forward and inverse transforms are opposite.