| 1958 | |
| 1959 | NOINLINE WARN_UNUSED_RESULT |
| 1960 | static int fftblue_fft(fftblue_plan plan, double c[], int isign, double fct) |
| 1961 | { |
| 1962 | size_t n=plan->n; |
| 1963 | size_t n2=plan->n2; |
| 1964 | double *bk = plan->bk; |
| 1965 | double *bkf = plan->bkf; |
| 1966 | double *akf = RALLOC(double, 2*n2); |
| 1967 | if (!akf) return -1; |
| 1968 | |
| 1969 | /* initialize a_k and FFT it */ |
| 1970 | if (isign>0) |
| 1971 | for (size_t m=0; m<2*n; m+=2) |
| 1972 | { |
| 1973 | akf[m] = c[m]*bk[m] - c[m+1]*bk[m+1]; |
| 1974 | akf[m+1] = c[m]*bk[m+1] + c[m+1]*bk[m]; |
| 1975 | } |
| 1976 | else |
| 1977 | for (size_t m=0; m<2*n; m+=2) |
| 1978 | { |
| 1979 | akf[m] = c[m]*bk[m] + c[m+1]*bk[m+1]; |
| 1980 | akf[m+1] =-c[m]*bk[m+1] + c[m+1]*bk[m]; |
| 1981 | } |
| 1982 | for (size_t m=2*n; m<2*n2; ++m) |
| 1983 | akf[m]=0; |
| 1984 | |
| 1985 | if (cfftp_forward (plan->plan,akf,fct)!=0) |
| 1986 | { DEALLOC(akf); return -1; } |
| 1987 | |
| 1988 | /* do the convolution */ |
| 1989 | if (isign>0) |
| 1990 | for (size_t m=0; m<2*n2; m+=2) |
| 1991 | { |
| 1992 | double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m]; |
| 1993 | akf[m ] = akf[m]*bkf[m] + akf[m+1]*bkf[m+1]; |
| 1994 | akf[m+1] = im; |
| 1995 | } |
| 1996 | else |
| 1997 | for (size_t m=0; m<2*n2; m+=2) |
| 1998 | { |
| 1999 | double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m]; |
| 2000 | akf[m ] = akf[m]*bkf[m] - akf[m+1]*bkf[m+1]; |
| 2001 | akf[m+1] = im; |
| 2002 | } |
| 2003 | |
| 2004 | /* inverse FFT */ |
| 2005 | if (cfftp_backward (plan->plan,akf,1.)!=0) |
| 2006 | { DEALLOC(akf); return -1; } |
| 2007 | |
| 2008 | /* multiply by b_k */ |
| 2009 | if (isign>0) |
| 2010 | for (size_t m=0; m<2*n; m+=2) |
| 2011 | { |
| 2012 | c[m] = bk[m] *akf[m] - bk[m+1]*akf[m+1]; |
| 2013 | c[m+1] = bk[m+1]*akf[m] + bk[m] *akf[m+1]; |
| 2014 | } |
| 2015 | else |
| 2016 | for (size_t m=0; m<2*n; m+=2) |
| 2017 | { |
no test coverage detected