1

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?

Paul R
  • 208,748
  • 37
  • 389
  • 560
Matt
  • 5,404
  • 3
  • 27
  • 39
  • If there's any more information I can give that would be helpful to solving the problem, just let me know and I would be happy to post it. – Matt Jul 29 '11 at 17:55
  • I changed my code so that it no longer converts the values into the domain [-1,1] and leaves them [-128,127] ... It still appears to not work. – Matt Jul 29 '11 at 18:54
  • if your sin-wave is phase shifted respect to its corespond frequency window or its period is not an multiply of any frequency window period then your DFT holds not a valid set of data for reconstruction which is basic problem of frequency domain signal reconstruction. To fix that you have to have a big enough number of samples and sometimes even shift signal before DFT/IDFT to synchronize your DFT/IDFT with it... – Spektre Feb 09 '14 at 15:42

2 Answers2

3

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.

Oliver Charlesworth
  • 267,707
  • 33
  • 569
  • 680
  • Wow... okay... first off, thanks... second, how do you break x[n] into x_real[n] and x_imag[n]? Third, x_real[n] is supposed to be a separate array than x_real[k]?or no? It seems recursively defined... – Matt Jul 29 '11 at 16:41
  • 1
    @Matt: Ah, I'm distinguishing between time-domain (`x`) and freq-domain (`X`) (this capitalisation difference is fairly common in signal-processing stuff). If your data is real-only (like an audio recording or something), then `x_real[n] = x[n]` and `x_imag[n] = 0`, so all the terms with `x_imag` disappear. – Oliver Charlesworth Jul 29 '11 at 16:47
  • My data is seemingly random and thus I am unsure if my data is real only. I mean... I only have one array of data so I (think) I could treat it all as real, but I'm unsure if that has complications. Does just having one array mean my data is all-real? – Matt Jul 29 '11 at 16:53
  • Also, I feel like all my math cred went down the drain for derping on Euler's formula. Oh well, at least it's getting fixed, thanks again. – Matt Jul 29 '11 at 16:54
  • Oh, also... my data is not from [-1,1] it's [-128,127] do I need to scale it down? or is it fine as it is? – Matt Jul 29 '11 at 16:56
1

As well as the points that @Oli correctly makes, you also have a fundamental misunderstanding about conversion between time and frequency domains. Your real input signal becomes a complex signal in the frequency domain. You should not be taking the magnitude of this and converting back to the time domain (this will actually give you the time domain autocorrelation if done correctly, but this is not what you want). If you want to be able to reconstruct the time domain signal then you must keep the complex frequency domain signal as it is (i.e. separate real/imaginary components) and do a complex-to-real IDFT to get back to the time domain.

E.g. your forward transform should look something like this:

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].real = real;
        windows[k].imag = image;
}

where windows is defined as an array of complex values.

Paul R
  • 208,748
  • 37
  • 389
  • 560
  • Makes sense. Then how do I reconstructed the signal. The forward transformation is certainly helpful, but it is only half of the problem. I think I can figure it out from @Oli's answer. Thanks for your help though. – Matt Aug 01 '11 at 14:46
  • 1
    The inverse transform is very similar to the forward transform, except that you have complex input data and you will just ignore the imaginary part of the output data (which should be close to zero anyway) and use the real part of the output. The biggest difference will be that you need to perform complex multiplication. – Paul R Aug 01 '11 at 16:27
  • by complex multiplication, you mean what @Oli had in his answer, correct? [following comment] – Matt Aug 01 '11 at 19:29
  • X_real[k] = SUM { x_real[n] * cos(2*pi*n*k/N) + x_imag[n] * sin(2*pi*n*k/N) } – Matt Aug 01 '11 at 19:30
  • 1
    Yes - I think there's also a sign change in at least one of the twiddle factors for the inverse transform, but it should be easy enough to look this up. – Paul R Aug 01 '11 at 20:53