Hi,

I need a pure julia FFT implementation as I am planning to use it in correspondence with the ForwardDiff package. Does this exist?

# Does there exists 1D FFT/IFFT implementation in pure Julia (i.e. no FFTW)

**vgdev**#1

**dlfivefifty**#2

You could adapt the `BigFloat`

fft implementation in FastTransforms.jl: https://github.com/MikaelSlevinsky/FastTransforms.jl/blob/master/src/fftBigFloat.jl

**vgdev**#4

Hmmm, it seems I will be running into problems with forwarddiffs dual numbers and the resulting complex numbers from the FFT.

Complex numbers are indeed not supported by ForwardDiff. If they were, a better way to differentiate through FFT would be to teach ForwardDiff the derivatives of FFTs: https://github.com/JuliaDiff/ForwardDiff.jl/pull/165 and refs

**stevengj**#6

Yes, since it is a linear transformation, the derivative of an FFT with respect to its inputs is simply another FFT.

**RoyiAvital**#7

For a vector x ou can always describe its Fourier Transform as y = D x where D is the DFT Matrix.

You can even generalize it for cases which are â€śZero Paddedâ€ť.

Since the Derivative is just {D}^{T} all you need is a function which generates those matrices.

Is there such built in function for that in Julia?

**vgdev**#8

Okey, so what I am really after is a way to find the derivative of the Volterra process \int_0^tg_\alpha(t-s) dW_s wrt. \alpha \in (-0.5,0) where g_\alpha(t-s) = (t-s)^\alpha and W a brownian motion. I am approximating this integral by looking at at the Riemann sum approximation \tilde{Y}(\frac{i}{n})= \sum_{k =1}^i g(\frac{b_k}{n})W_{i-k} on an n sized grid, i = 0,1,2,\dots. This can be interpreted as a discrete convolution \tilde{Y}(\frac{i}{n}) = (g*W)_i, which can be computed efficiently by FFT (Ihave gotten this to work with ForwardDiff using a direct naive convolution, but its too slow). Furthermore, if I need to do several runs of this numerical scheme to find an estimate of \mathbb{E}[\tilde{Y}(\frac{i}{n})] \sim \frac{1}{N}\sum_j^N\tilde{Y}^{(j)}(\frac{i}{n}). I have to perform this convolution each time I simulate a path of this process. So to speed up things I would like to compute the FFT of g(b_k/n) once (since its deterministic) and multiply it with the FFT of a new set of random variates W, and then take the IFFT to arrive at the convolution.

The code to my current implementation (with some slight changes to the actual approximation scheme):

```
function gen_volterra(a,s,n,G,rand_1,rand_2)
Z_1 = zeros(eltype(a),s); #first normal
Z_2 = zeros(eltype(a),s); # 2nd normal
#G = zeros(eltype(a),s+1); #Kernel function
@fastmath @inbounds @simd for i = 1:s
#G[i+1] = b(a,i+1)/(n^a); #fill kernel vector <<----this is where the error occurs
Z_1[i] = get_sigma_Z1(n)*rand_1[i] #set correct variance of random variate
Z_2[i] = get_sigma_Z2(a,n)*(get_rho(a)*rand_1[i]+sqrt(1-get_rho(a)*get_rho(a))*rand_2[i]); #set correct variance of random variate
end
return [0;(direct_conv(Z_1,G)[1:s]+Z_2)]; #Convolute
#return G
end
function direct_conv{T}(a::Array{T}, b::Array{T})
m = length(a)
n = length(b)
c = zeros(T,m+n-1)
@fastmath @inbounds @simd for j=1:m
@fastmath @inbounds @simd for k=1:n
c[j+k-1] += a[j]*b[k]
end
end
return c
end
```

and addtional helper functions:

```
function b(a, k)
return (k^(a+1.0)-(k-1.0)^(a+1.0))/(a+1.0);
end
function get_sigma_Z1(n)
return 1.0/sqrt(n);
end
function get_sigma_Z2(a,n)
return 1.0/sqrt(((2.0*a+1.0)*n^(2.0*a+1.0)))
end
function get_rho(a)
return sqrt(2.0*a+1.0)/(a+1.0)
end
function gen_G(a,n,s)
i = range(1,s)
return [0;map(i->b(a,i+1)/(n^a),i)]
end
```

Example of generating 1 path:

```
julia> a= -0.43
julia> n= 500
julia> s= 500
julia> G = gen_G(a,n,s)
julia> rand_1 = randn(s)
julia> rand_2 = randn(s)
julia> path = gen_volterra(a,s,n,G,rand_1,rand_2)
```

But the operation below will then be quite slow due to the number of floating point operations in the direct_conv method.

```
m = ForwardDiff.derivative(a ->gen_volterra(a,s,n,G,rand_1,rand_2),-0.43)
```

**stevengj**#9

You can just FFT the columns of the identity matrix. But usually it is preferable to work with the DFT matrix implicitly, since multiplying by it (or its transpose â€” it is symmetric) is computed efficiently by the `fft`

function (and the inverse or conjugate by `ifft`

).

**antoine-levitt**#10

Why do you need AD for this? Just do the derivative wrt alpha, which is again a convolution that you can do with an FFT.

**vgdev**#11

Because this process is used as an input in a SDE for another process, which I didnâ€™t include here for clairity. It is actually wrt to this other process I want the gradient, so the stuff outlined above is a subproblem as this is where the performance bottleneck is. But a solution might be to use this