# Repeated Convolutions With Large 1D Arrays

**URL:** <https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194>\
**Category:** Performance\
**Tags:** question, dsp, fft\
**Created:** [March 20, 2025, 7:33pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194 "2025-03-20T19:33:57Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [March 20, 2025, 7:33pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/1 "2025-03-20T19:33:57Z")

</div>

What’s the optimal way to compute repeated convolutions with large 1D arrays? the filters, (g\_x,g\_y, & g\_z) stay the same each iteration, but the “signals” change. To be clear, the computation of the f’s is fast, the convolutions are what is slowing me down.

```julia
for i in 1:iters

    #input function for x
    f_x .= x_fun(s,ẑ, a_xx,b_zx) # 2500001-element Vector{Float64}
    #input function for y
    f_y .= y_fun(x̂,new_z, a_xy, b_zy) # 2500001-element Vector{Float64} 
    #input function for z 
    f_z .= z_fun(ŷ, a,b, z_p) # 2500001-element Vector{Float64}

    #Get estimates
    x̂ .= fftfilt(g_x,δ*f_x)[1:end] # g_x: 1001-element Vector{Float64} 
    ŷ .= fftfilt(g_y,δ*f_y)[1:end] # g_y: 1001-element Vector{Float64}
    ẑ .= fftfilt(g_z,δ*f_z)[1:end] # g_z: 100001-element Vector{Float64}
    
end

```

Thanks,

DS

---

<div class="post-metadata">

**Author:** ![Jeff\_Emanuel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff_emanuel/32/15440_2.png) [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)\
**Post date:** [March 20, 2025, 8:21pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/2 "2025-03-20T20:21:23Z")

</div>

This isn’t going to be the improvement you need, but you can omit `[1:end]` because it is completely unnecessary. It’s making a copy of the filter result before the element-wise assignment to the left-hand side. It costs time and increases memory pressure and GC overhead.

You could use `fftfilt!` to put the filter results directly into your left-hand variables and avoid `fftfilt`’s internal result allocation.

There are also small bits of redundant work in `fftfilt`. Each call redundantly computes the optimal FFT length, plans FFTs, and applies the FFT to the `g` filters. These could be extracted, but they are still small compared to the multitude of FFT/inverse FFT applications on the overlapped windows of the signals. Likely the biggest possible gain would be to parallelize those.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [March 20, 2025, 8:52pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/3 "2025-03-20T20:52:38Z")

</div>

Ahh. Thx for the tips, [Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)! Will implement. Didn’t realize the [1:end] was an issue. Using threads seems like an obvious choice - will give it a shot.

What about precomputing the FFT’s for the filters outside of the loop with FFTW.jl? Or using Metal.jl (I have a mac)?

Thanks again,

DS

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 20, 2025, 10:07pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/4 "2025-03-20T22:07:17Z")

</div>

> [@donkeysaddle](#):
>
> What about precomputing the FFT’s for the filters outside of the loop with [FFTW.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/FFTW)? Or using [Metal.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Metal) (I have a mac)?

Yes, I would definitely recommend this — that will save you one of the 3 FFTs (i.e. ≈ 33%). Also precomputing an FFTW “plan” (and adding the `FFTW.MEASURE` or `FFTW.PATIENT` option, which takes longer to create but runs faster). Use real-to-complex and complex-to-real FFTs if your data are real.

(And if you have some freedom to choose the transform length, choose a highly composite size, i.e. factors of 2,3,5. Using a length `2500001`, which has a large prime factor of `18797`, will be many times slower for FFTs. If you are doing a linear convolution you should be zero-padding anyway, which lets you choose the length to pad to.)

Note that `δ*f_x` allocates a new array, so I would include `δ .* …` in the definition of your arrays. And `x̂ .= fftfilt(...)` still allocates a new array for the rhs. But if you use a precomputed-FFT plan, you can act it in-place on precomputed outputs with `mul!`.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [March 28, 2025, 7:41pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/5 "2025-03-28T19:41:21Z")

</div>

Thank you, @stevengj and @Jeff_Emanuel, for the thoughtful responses. Incredibly helpful. Really appreciate it. I was able to speedup the code considerably by following your advice and using threads.

This code is considerably faster:

```julia
FFTW.set_num_threads(Threads.nthreads() - 3) # 3 threads will be used for inner loop below

    for i in 1:iters        
        #Get estimates using threads to multiprocess
        Threads.@threads for j in 1:3 
            if j == 1 
                #input function for x
                f_x .= x_fun(s,ẑ,a_xx,b_zx) # f_x: 2_500_000-element Vector{Float64}                
                fftfilt!(x̂,g_x,f_x) # g_x: 1000-element Vector{Float64}
            elseif j == 2
                #input function for y
                f_y .= y_fun(x̂,ẑ, a_xy,b_zy) # 2_500_000-element Vector{Float64} 
                fftfilt!(ŷ,g_y,f_y) # g_y: 1000-element Vector{Float64}                
            else 
                #input function for z
                f_z .= z_fun(ŷ,a,b,z_p) # 2_500_000-element Vector{Float64}
                fftfilt!(ẑ,g_z,f_z) # g_z: 100000-element Vector{Float64}
            end
        end
    end

```

Interestingly, I didn’t notice any difference when I changed the vector lengths. I also have been unable to see much of an improvement by computing the FFT’s for the filters outside of the loops. I am sure I am doing something (many things…) wrong, but I have no idea what. Any advice would be greatly appreciated.

I have included a comparison example where despite using `FFTW.PATIENT ` and precomputing the FFT for the filter outside the loop I am unable to beat `fftfilt`:

```julia
using BenchmarkTools,DSP,FFTW,LinearAlgebra

end_time = 25000 - .01 # end time for signal
filt_end = 1000 - .01 # end time for filter
δ = .01 # granularity
g = 100.0*exp.(-100.0*(0.0:δ:filt_end)); #filter, length 100_000, 

# Method 1 Precompute Patient Plan:
x = sin.(0.0:δ:end_time).^2; #signal, length 2_500_000

lx =length(x) #length of signal: 2_500_000

# Compute the full convolution length
N = lx + length(g) - 1; #length of arrays minus 1

# Create zero-padded arrays:
xp = zeros(Float64, N); 
gp = zeros(Float64, N);

#populate padded arrays
xp[1:lx] .= x;
gp[1:length(g)] .= g;

# Create FFT plans with PATIENT flag:
plan_x = plan_rfft(xp, flags=FFTW.PATIENT);
plan_g = plan_rfft(gp, flags=FFTW.PATIENT);

# Precompute FFT's:
G_fft = plan_g * gp;
X_fft = plan_x * xp;

conv_full_x = irfft(X_fft .* G_fft, N)[1:lx]; # compute convolution

@benchmark begin
    for i in 1:10
        #update x
        xp[1:lx] .+= sin.(0.0:δ:end_time).^2

        #calculate x fft
        mul!(X_fft,plan_x, xp)
        
        #inversre fft
        conv_full_x .= irfft(X_fft .* G_fft, N)[1:lx];
    end
end

#Method 2 fftfilt!:
x_filt = sin.(0.0:δ:end_time).^2 #signal array length 2_500_000

conv_full_x_filt = fftfilt(g,x); # compute convolution 

@benchmark begin
    for i in 1:10   
        #update x
        x_filt .+= sin.(0.0:δ:end_time).^2
        
        #convolve 
        fftfilt!(conv_full_x_filt,g,x_filt)
    end
end

```

Thanks,

DS

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 28, 2025, 9:01pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/6 "2025-03-28T21:01:46Z")

</div>

> [@donkeysaddle](#):
>
> I have included a comparison example where despite using `FFTW.PATIENT ` and precomputing the FFT for the filter outside the loop I am unable to beat `fftfilt`:

Don’t benchmark loops in global scope. Performance-sensitive code should be in a function. This is the number one rule when using Julia — it is literally the [first performance tip](https://docs.julialang.org/en/v1/manual/performance-tips/).

> [@donkeysaddle](#):
>
> ```julia
> # Compute the full convolution length
> N = lx + length(g) - 1; #length of arrays minus 1
> 
> ```

You are using `N = 2599999`, which is a prime number. You should zero-pad to the next highly composite size ≥ this.

> [@donkeysaddle](#):
>
> ```julia
> #inversre fft
> conv_full_x .= irfft(X_fft .* G_fft, N)[1:lx];
> 
> ```

You should also precompute the `irfft` plan. (Also, I would use a `brfft` instead, and include the 1/N normalization factor in the precomputation of `G_fft`.

The `.=` is not helping you here — there’s nothing for it to fuse with on the right-hand side. You’re still allocating lots of temporary arrays (for `X_fft .* G_fft`, and for the slice `[1:lx]` since you’re not using views. Read the Julia performance tips.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [April 3, 2025, 3:57pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/7 "2025-04-03T15:57:17Z")

</div>

@stevengj thank you for the lesson(s)! Code is significantly faster. Appreciate your time and patience. Incredibly helpful.

DS

Code for example discussed:

```julia
using AbstractFFTs, BenchmarkTools, DSP, FFTW, LinearAlgebra

function fftw_test()
    end_time = 25000 - .01 
    filt_end = 1000 - .01 
    δ = .01
    g = 100.0*exp.(-100.0*(0.0:δ:filt_end)); #Filter array, length 100_000

    x = sin.(0.0:δ:end_time).^2; #Input signal array, length 2_500_000
    lx = length(x) #2_500_000

    #Pad size, greater than min required size, but factor of small primes 
    N = lx + length(g); 

    #Create zero-padded vectors:
    xp = zeros(N); 
    gp = zeros(N);

    xp[1:lx] .= x;
    gp[1:length(g)] .= g;

    #Create FFT plans with PATIENT flag:
    plan_x = plan_rfft(xp, flags=FFTW.PATIENT);
    plan_g = plan_rfft(gp, flags=FFTW.PATIENT);

    fft_length = Int(N/2+1) #Length of FFT vectors
    
    #Vectors for FFT's 
    G_fft = Vector{ComplexF64}(undef,fft_length) 
    X_fft = Vector{ComplexF64}(undef,fft_length)

    #Compute FFT's:
    mul!(G_fft,plan_g,gp/N)
    mul!(X_fft,plan_x,xp)

    #Create inverse FFT plan with PATIENT flag:
    plan_xi = plan_brfft(X_fft,N, flags=FFTW.PATIENT)
   
    #Vector to store solutions
    conv_full_xb = Vector{Float64}(undef,N)

    #Convolution 
    X_fft .*= G_fft
    
    #inverse FFT
    mul!(conv_full_xb,plan_xi,X_fft)

    for i in 1:10
        #Update 
        xp[1:lx] .+= x.^2

        #Perform FFT
        mul!(X_fft,plan_x,xp)
        
        #Convolve
        X_fft .*= G_fft
        
        #Inverse FFT
        mul!(conv_full_xb,plan_xi,X_fft)
    end

    return conv_full_xb[1:lx]

end

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [April 3, 2025, 4:55pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/8 "2025-04-03T16:55:59Z")

</div>

> [@donkeysaddle](#):
>
> `fft_length = Int(N/2+1) #Length of FFT vectors`

This is not a safe way to calculate the (integer) length. It does a floating point division followed by conversion to `Int`, but will fail for odd N.

Julia supports native integer arithmetic, which is faster and better for this usecase. Here you should use `div` to perform integer division:

```julia
fft_length = div(N, 2) + 1

```

or

```julia
fft_length = (N >> 1) + 1

```

This gives a valid fft length for both even and odd lengths (though even-length rfft is anyway faster).

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [April 3, 2025, 4:58pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/9 "2025-04-03T16:58:57Z")

</div>

> [@donkeysaddle](#):
>
> `conv_full_xb[1:lx]`

Alternatively, you could do

```julia
resize!(conv_full_xb, lx)

```

to avoid allocating new memory.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [April 3, 2025, 5:08pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/10 "2025-04-03T17:08:12Z")

</div>

DNF, thank you! This makes sense.

Would `@view conv_full_xb[1:lx]` also be acceptable?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [April 3, 2025, 5:30pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/11 "2025-04-03T17:30:38Z")

</div>

I guess… But then you are leaving a ‘hidden’ chunk of memory around under the hood, which seems a bit less pleasing to me. With `resize!` you free up that piece of memory, and you have a plain Vector, which is perhaps ‘cleaner’. What if you pass the vector on to some code that only accepts Vector?

It’s possibly not a very big deal, though.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [April 3, 2025, 5:38pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/12 "2025-04-03T17:38:27Z")

</div>

Makes sense! Thanks!

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [April 4, 2025, 9:05pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/14 "2025-04-04T21:05:14Z")

</div>

Regarding `fft_length`: The [`nextprod`](https://docs.julialang.org/en/v1/base/math/#Base.nextprod) function is perfect for finding a good FFT length. E.g.

```julia
fft_length = nextprod((2, 3, 5), N ÷ 2 + 1)

```

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [April 4, 2025, 9:26pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/15 "2025-04-04T21:26:34Z")

</div>

@PeterSimon, Wouldn’t it be `N = nextprod((2, 3, 5), N ÷ 2 + 1)`? If so, it actually seems to increase the time slightly. I assume the added “efficiency” is outweighed by the increased number of computations (N is already divisible by 2 and 5). But that’s just a guess, as you can see from this thread, I really don’t know what the hell I am talking about…

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [April 4, 2025, 9:52pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/16 "2025-04-04T21:52:21Z")

</div>

For specific filters, like `X` and `G` in the example code, there are tailored algos which would be closer to O(N) instead of the FFT-based O(N\*log(N)) and possibly with better constants.  
So if the filters are general, FFT is definitely the way to go, otherwise this lemon might be squeezed more.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [April 4, 2025, 10:29pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/17 "2025-04-04T22:29:23Z")

</div>

> [@donkeysaddle](#):
>
> Wouldn’t it be `N = nextprod((2, 3, 5), N ÷ 2 + 1)`?

Sorry, you’re correct.

> [@donkeysaddle](#):
>
> If so, it actually seems to increase the time slightly.

Your original value of `N` is `2_600_000`, so assuming that your expression `N ÷ 2 + 1` is the minimum allowable length of the padded arrays this is `1_300_001` and `N = nextprod((2, 3, 5), N ÷ 2 + 1)` yields a new value of `N = 1_310_720`. I am surprised to hear that using this length instead of `2_600_000` is slower, since it’s about half as large. And checking the factors:

```julia
julia> using Primes

julia> factor(2_600_000)
2^6 * 5^5 * 13

julia> factor(1_310_720)
2^18 * 5

```

it still seems like the result of `nextprod` should be faster, since it only involves 2 prime factors, rather than 3. Maybe others can explain this effect.

---

<div class="post-metadata">

**Author:** ![donkeysaddle](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@donkeysaddle](https://discourse.julialang.org/u/donkeysaddle)\
**Post date:** [April 8, 2025, 4:13pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/18 "2025-04-08T16:13:26Z")

</div>

@Dan cool! Where might I find these tailored algos?

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [April 8, 2025, 5:12pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/19 "2025-04-08T17:12:17Z")

</div>

With a fixed kernel it’s trivial to make a direct implementation of the convolution in linear time, but at these sizes it won’t be competitive with an FFT implementation. To improve the speed you need look into IIR filters. Those are generally more difficult to design but if `g = 100.0*exp.(-100.0*(0.0:δ:filt_end))` is a representative filter you are in luck, as this is the simplest of all IIR filters.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [April 8, 2025, 7:32pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/20 "2025-04-08T19:32:38Z")

</div>

> [@donkeysaddle](#):
>
> ```julia
> filt_end = 1000 - .01 
> δ = .01
> g = 100.0*exp.(-100.0*(0.0:δ:filt_end)); #Filter array, length 100_000
> 
> ```

The IIR approach is likely still best but it seems excessive to use a kernel of length 100000 when 99254 of the coefficients are exactly zero and all but the 29 first are smaller than 1e-10. Maybe this isn’t a real world example.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [April 8, 2025, 7:46pm UTC](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194/21 "2025-04-08T19:46:27Z")

</div>

Just to make things concrete regarding faster filter, here is an example:

```julia
N = 1000 # length of vector
n = 100 # length of kernel

v = rand(N) # a bit of random data to process

# naive implementation of convolution (delta = 0.1) 
ema = [sum(v[i-j]*exp(-0.1*j) for j in 0:(n-1)) for i in n:N];

@show ema[1]
@show length(ema)

# alternative implementation without FFT, using O(N) time
ema2 = similar(ema);
ma = sum(v[100-j]*exp(-0.1*j) for j in 0:(n-1))

ema2[1] = ma

for i in 2:(N-n+1)
    global ma
    ma -= v[i-1]*exp(-0.1*(n-1))
    ma *= exp(-0.1)
    ma += v[i+n-1]
    ema2[i] = ma
end

@show ema2[1]
@show ema2[2]
@show ema[2]
@show ema ≈ ema2

```

BTW EMA stands for exponential moving average, which is a common name for this kernel.

The other kernel in the OP can also have a tailored algorithm.

[Next page](https://discourse.julialang.org/t/repeated-convolutions-with-large-1d-arrays/127194.md?page=2)