DST question

Hi,
with the help of the english wikipedia entry[1] about the
discrete sine transform, I made a little program[2] to
calculate the DST, but it doesn't quite work.
(The inverse DST of the DST doesn't yield the original
data)

The DST-routine is quite straightforward: (it saves the discrete
sine transform in the DST array)

void discrete_sine_transform(double *DST, double *samples, int N) {
    int bin, k;    /* bin means "container", "bucket" */
    double arg;
    for (bin=0; bin<N; bin++) {
        DST[bin] = 0;
        for (k=0; k<N; k++) {
            arg = M_PI * (bin+1)*(k+1)/(N+1);
            DST[bin] += samples[k] * sin(arg);
        }
    }
}

And according to the wikipedia article, the inverse
of this transform is simply the transform times the
factor 2/(N+1). So my inverse function is: (it tries to
recreate the samples array by means of the DST array)

void inverse_discrete_sine_transform(double *DST, double *samples, int N) {
    int k;
    discrete_sine_transform(samples, DST, N);
    for (k=0; k<N; k++) {
        samples[k] *= 2.0/(N+1);
    }
}

The program takes the input array as parameters on the command line and
continues to output the calculated DST and tries to invert the operation
by reconstructing the input via the calculated DST.

Now, shouldn't DST*DST_inverse yield the identity? And thus shouldn't the
"reconstructed" data be equivalent to the original data?

Because it isn't. Here are a few example outputs:

$ ./dst 1 2 3
original data:      { 1.00 2.00 3.00 }
the DST:            { 4.62 -3.62 1.97 }
reconstructed data: { 0.52 1.32 2.51 }

$ ./dst 1 0 0
original data:      { 1.00 0.00 0.00 }
the DST:            { 1.21 -0.21 1.56 }
reconstructed data: { 0.88 -0.18 0.66 }

Curiously, for an input-length of 2 it seems to work:

$ ./dst 1 2
original data:      { 1.00 2.00 }
the DST:            { 2.60 -0.87 }
reconstructed data: { 1.00 2.00 }

$ ./dst 17 23
original data:      { 17.00 23.00 }
the DST:            { 34.64 -5.20 }
reconstructed data: { 17.00 23.00 }


So, is the wiki-article wrong? Or my assumption that the inverse DST of
the DST should yield the identity?
(Or is there a bug in my program? Which I don't think,
because I basically just copied the formula from wikipedia.)


[1] http://en.wikipedia.org/wiki/Discrete_sine_transform (I used the
first definition)
[2] http://nopaste.info/2f00cb7f9b.html

On 12 Mar, 09:47, Newbie <no...@nospam.invalid> wrote:

> So, is the wiki-article wrong? Or my assumption that the inverse DST of
> the DST should yield the identity?
> (Or is there a bug in my program? Which I don't think,
> because I basically just copied the formula from wikipedia.)

The fact that your program works for the N=2 case indicates
that the formulas that *are* on the wiki page are correct.
It does *not* mean these formulas are *complete*. There are
likely technical details missing that are necessary to the
algorithm work for N odd.

I don't know anything about the DST or DCT, but FFTs for N
a power of 2 is standard material in DSP entry / medium level
class. Generalizing the FFT to general N is a highly specialized
skill that very few people do, and taht even fewer want to do.
Or, for that matter, have the skills to pull off successfully.

Rune

Newbie <noone@nospam.invalid> wrote:

> with the help of the english wikipedia entry[1] about the
> discrete sine transform, I made a little program[2] to
> calculate the DST, but it doesn't quite work.
> (The inverse DST of the DST doesn't yield the original
> data)
 
> The DST-routine is quite straightforward: (it saves the discrete
> sine transform in the DST array)

There is something wrong somewhere else.  Show the whole
program, including the program that calls the DST routine.

This one, close to yours, except for the correction factor, works.
I believe that you have something wrong in the calling routine
such that you aren't printing what you think you are.

Note that the DST-I of {1,0,0} is {sqrt(0.5),1.0,sqrt(0.5)}
unlike the one that you show.

#include <stdio.h>
#include <math.h>

int main() {
   int i,j,k;
   double x[3]={1,2,3},y[3];
   for(i=0;i<3;i++) printf("%g ",x[i]);
   printf("\n");
   dst(y,x,3);
   for(i=0;i<3;i++) printf("%g ",y[i]);
   printf("\n");
   dst(x,y,3);
   for(i=0;i<3;i++) printf("%g ",x[i]);
   printf("\n");
   }

dst(double *out, double *in, int N) {
   int bin,k;
   double arg;
   for(bin=0;bin<N;bin++) {
     out[bin]=0;
     for(k=0;k<N;k++) out[bin] += in[k]*sin(M_PI*(bin+1)*(k+1)/(N+1));
     }

-- glen 

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> There is something wrong somewhere else.  Show the whole
> program, including the program that calls the DST routine.

I did, see the second link at the end of the initial post.

> I believe that you have something wrong in the calling routine
> such that you aren't printing what you think you are.

Indeed. In my naivety I treated doubles and floats interchangably.
Specificly, I stored the return value of atof() into a double and
later displayed that double with a "%f" format specifier.
So it was a bug in my program after all.

It now works, though with new "normalizing factors", which I got from
here[1].

(The "forward" transform is multiplied by 2, the inverse by 1/(2N+2))

Thanks for the responses guys.

[1] http://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html

Newbie <noone@nospam.invalid> wrote:
(snip)
 
> Indeed. In my naivety I treated doubles and floats interchangably.
> Specificly, I stored the return value of atof() into a double and
> later displayed that double with a "%f" format specifier.
> So it was a bug in my program after all.

Both of those should work fine.  If you pass (float) to printf,
it will be converted, as all (float) arguments to varargs functions
are, to (double).  The problem should be somewhere else.

-- glen

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> The problem should be somewhere else.

Well, perhaps it's the sin() function or the M_PI from
math.h which causes this behaviour.

All I know is that if you take my initial
program (the one on nopaste) and simply substitute every
"double" with "float", it magically starts to work
as expected.

So, I don't know, somewhere there's a problem with
the doubles.

Greetings

On 13 Mar, 08:21, Newbie <no...@nospam.invalid> wrote:

> All I know is that if you take my initial
> program (the one on nopaste) and simply substitute every
> "double" with "float", it magically starts to work
> as expected.
>
> So, I don't know, somewhere there's a problem with
> the doubles.

Check the type declarations in the calling routine.
Do you call your routine with float arrays instead
of double arrays?

It's been a long time since I programmed C, so I don't
remember if the compiler is supposed to complain if
you call the routine with the wrong types of arguments.

Rune

Newbie <noone@nospam.invalid> wrote:
> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
>> The problem should be somewhere else.
 
> Well, perhaps it's the sin() function or the M_PI from
> math.h which causes this behaviour.
 
> All I know is that if you take my initial
> program (the one on nopaste) and simply substitute every
> "double" with "float", it magically starts to work
> as expected.

 input = malloc(N*sizeof(int));
 DST  = malloc(N*sizeof(int)); 

It might work better with sizeof(double) since input and DST
are (double)

-- glen

 

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> It might work better with sizeof(double) since input and DST
> are (double)

Ha! What a stupid mistake. The very first run I actually used integers,
but realized directly: woah, integers don't make any sense here. So I
*manually* updated the "int"s to "double"s. Then later I did change
types automatically, so I didn't bother to check everything and just
skimmed over the malloc()-line every time..

Note: on my system sizeof(float)==sizeof(int) evaluates to true (both 4),
while sizeof(double)==8; so that explains why it accidently worked with floats..

Thanks for the input, guys.

Rune Allnor wrote:
> On 13 Mar, 08:21, Newbie <no...@nospam.invalid> wrote:
> 
>> All I know is that if you take my initial
>> program (the one on nopaste) and simply substitute every
>> "double" with "float", it magically starts to work
>> as expected.
>>
>> So, I don't know, somewhere there's a problem with
>> the doubles.
> 
> Check the type declarations in the calling routine.
> Do you call your routine with float arrays instead
> of double arrays?
> 
> It's been a long time since I programmed C, so I don't
> remember if the compiler is supposed to complain if
> you call the routine with the wrong types of arguments.
> 
> Rune

If it is a GNU derivative compiler, turning on -Wall will
enable all warnings. This should highlight all such problems.
Emphasis "should".

--
Les Cargill