开发者

Fast way to sum Fourier series?

开发者 https://www.devze.com 2023-04-09 09:38 出处:网络
I have generated the coefficients using FFTW, now I want to reconstruct the original data, but using only the first numCoefs coefficients rather than all of them. At the moment I\'m using the below co

I have generated the coefficients using FFTW, now I want to reconstruct the original data, but using only the first numCoefs coefficients rather than all of them. At the moment I'm using the below code, which is very slow:

for ( unsigned int i = 0; i < length; ++i )
{
    double sum = 0;
    for ( unsigned int j = 0; j < numCoefs; ++j )
    {
        sum开发者_如何学Python += ( coefs[j][0] * cos( j * omega * i ) ) + ( coefs[j][1] * sin( j * omega * i ) );
    }
    data[i] = sum;
}

Is there a faster way?


A much simpler solution would be to zero the unwanted coefficients and then do an IFFT with FFTW. This will be a lot more efficient than doing an IDFT as above.

Note that you may get some artefacts in the time domain when you do this kind of thing - you're effectively multiplying by a step function in the frequency domain, which is equivalent to convolution with a sinc function in the time domain. To reduce the resulting "ringing" in the time domain you should use a window function to smooth out the transition between non-zero and zero coeffs.


If your numCoefs is anywhere near or greater than log(length), then an IFFT, which is O(n*log(n)) in computational complexity, will most likely be faster, as well as pre-optimized for you. Just zero all the bins except for the coefficients you want to keep, and make sure to also keep their negative frequency complex conjugates as well if you want a real result.

If your numCoefs is small relative to log(length), then other optimizations you could try include using sinf() and cosf() if you don't really need more than 6 digits of precision, and pre-calculating omega*i outside the inner loop (although your compiler should do be doing that for you unless you have the optimization setting low or off).

0

精彩评论

暂无评论...
验证码 换一张
取 消

关注公众号