help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] normalization factor in Discrete Sine Transformation


From: Peter Melchior
Subject: [Help-gsl] normalization factor in Discrete Sine Transformation
Date: Wed, 16 Nov 2005 17:48:57 +0100

Hello everybody,

I implemented a discrete sine transformation using GSL's FFT for real
functions and the standard algorithm for storing the real numbers in the
complex array for doing the complex FFT.
For the discrete cosine transformation it works perfectly, but for the
sine transformation I get an additional factor 4, when doing the
transformation and its inverse after another.
This factor 4 doesn't depend on the length of the data set, so its not a
normal normalization issue.

Any suggestions what can cause this problem?

Best regards,

Peter Melchior



Here's the code:

void fft_sin_forward(NumVector<double>& f, NumVector<double>& F) {
  int J = f.size() - 1;
  int N = 2*J;
  NumVector<Complex> Fcomplex(N), fcomplex(N);
  for (int i =1; i < J; i++) {
    fcomplex(i) = f(i);
    fcomplex(N-i) = -f(i);
  }
  fft_forward(fcomplex,Fcomplex);
  F(0) = F(J) = 0;
  for (int i =1; i < J; i++)
    F(i) = -2 * imag(Fcomplex(i));
  // don't know why but I need factor 4 here
  F /= 4;
}

void fft_sin_inverse(NumVector<double>& F, NumVector<double>& f) {
  int J = F.size() - 1;
  int N = 2*J;
  NumVector<Complex> Fcomplex(N), fcomplex(N);
  for (int i = 1; i < J; i++) {
    Fcomplex(i) = -2.*I*F(i);
    Fcomplex(N-i) = 2.*I*F(i);
  }
  fft_inverse(Fcomplex,fcomplex);
  f(0) = f(J) = 0;
  for (int i =1; i < J; i++) 
    f(i) = real(fcomplex(i));
}




reply via email to

[Prev in Thread] Current Thread [Next in Thread]