How can I transfer a discrete set of data into the frequency domain and back (preferrably losslessly)

Question

I would like to take an array of bytes of roughly size 70-80k and transform them from the time domain to the frequency domain (probably using a DFT). I have been following wiki and gotten this code so far.

for (int k = 0; k < windows.length; k++) {
        double imag = 0.0;
        double real = 0.0;
        for (int n = 0; n < data.length; n++) {
            double val = (data[n])
                    * Math.exp(-2.0 * Math.PI * n * k / data.length)
                    / 128;
            imag += Math.cos(val);
            real += Math.sin(val);
        }
        windows[k] = Math.sqrt(imag * imag + real
                * real);
}

and as far as I know, that finds the magnitude of each frequency window/bin. I then go through the windows and find the one with the highest magnitude. I add a flag to that frequency to be used when reconstructing the signal. I check to see if the reconstructed signal matches my original data set. If it doesn't find the next highest frequency window and flag that to be used when reconstructing the signal.

Here is the code I have for reconstructing the signal which I'm mostly certain is very wrong (it is supposed to perform an IDFT):

for (int n = 0; n < data.length; n++) {
        double imag = 0.0;
        double real = 0.0;
        sinValue[n] = 0;
        for (int k = 0; k < freqUsed.length; k++) {
            if (freqUsed[k]) {
                double val = (windows[k] * Math.exp(2.0 * Math.PI * n
                        * k / data.length));
                imag += Math.cos(val);
                real += Math.sin(val);
            }
        }
        sinValue[n] = imag* imag + real * real;
        sinValue[n] /= data.length;
        newData[n] = (byte) (127 * sinValue[n]);
}

freqUsed is a boolean array used to mark whether or not a frequency window should be used when reconstructing the signal.

Anyway, here are the problems that arise:

1. Even if all of the frequency windows are used, the signal is not reconstructed. This may be due to the fact that ...
2. Sometimes the value of Math.exp() is too high and thus returns infinity. This makes it difficult to get accurate calculations.
3. While I have been following wiki as a guide, it is hard to tell whether or not my data is meaningful. This makes it hard to test and identify problems.

Off hand from the problem:

I am fairly new to this and do not fully understand everything. Thus, any help or insight is appreciated. Thanks for even taking the time to read all of that and thanks ahead of time for any help you can provide. Any help really would be good, even if you think I'm doing this the most worst awful way possible, I'd like to know. Thanks again.

-

EDIT:

So I updated my code to look like:

for (int k = 0; k < windows.length; k++) {
        double imag = 0.0;
        double real = 0.0;
        for (int n = 0; n < data.length; n++) {
            double val = (-2.0 * Math.PI * n * k / data.length);
            imag += data[n]*-Math.sin(val);
            real += data[n]*Math.cos(val);
        }
        windows[k] = Math.sqrt(imag * imag + real
                * real);
}

for the original transform and :

for (int n = 0; n < data.length; n++) {
    double imag = 0.0;
    double real = 0.0;
    sinValue[n] = 0;
    for (int k = 0; k < freqUsed.length; k++) {
        if (freqUsed[k]) {
            double val = (2.0 * Math.PI * n
                    * k / data.length);
            imag += windows[k]*-Math.sin(val);
            real += windows[k]*Math.cos(val);
        }
    }
    sinValue[n] = Math.sqrt(imag* imag + real * real);
    sinValue[n] /= data.length;
    newData[n] = (byte) (Math.floor(sinValue[n]));
}

for the inverse transform. Though I am still concerned that it isn't quite working correctly. I generated an array holding a single sine wave and it can't even decompose/reconstruct that. Any insight as to what I'm missing?

Answer 1

Yes, your code (for both DFT and IDFT) is broken. You are confusing the issue of how to use the exponential. The DFT can be written as:

       N-1
X[k] = SUM { x[n] . exp(-j * 2 * pi * n * k / N) }
       n=0

where j is sqrt(-1). That can be expressed as:

       N-1
X[k] = SUM {   (x_real[n] * cos(2*pi*n*k/N) + x_imag[n] * sin(2*pi*n*k/N))
       n=0  +j.(x_imag[n] * cos(2*pi*n*k/N) - x_real[n] * sin(2*pi*n*k/N)) }

which in turn can be split into:

            N-1
X_real[k] = SUM { x_real[n] * cos(2*pi*n*k/N) + x_imag[n] * sin(2*pi*n*k/N) }
            n=0

            N-1
X_imag[k] = SUM { x_imag[n] * cos(2*pi*n*k/N) - x_real[n] * sin(2*pi*n*k/N) }
            n=0

If your input data is real-only, this simplifies to:

            N-1
X_real[k] = SUM { x[n] * cos(2*pi*n*k/N) }
            n=0

            N-1
X_imag[k] = SUM { x[n] * -sin(2*pi*n*k/N) }
            n=0

So in summary, you don't need both the exp and the cos/sin.