# How to verify the convolution theorem in Julia?

I’ve trying to prove convolution theorem in Julia but without success. Say I have two function f1 = sin(30πt) & f2 = sin(10πt). Then, I have the following code

``````fconv_fda = ifft(fft(f1).*fft(f2))
fconv_tda = DSP.conv(f1,f2)[1:length(f1)]
plot(fconv_tda)
plot!(real.(fconv_fda))
``````

But both the plots turn out to be different. What am I doing wrong here?

you can’t prove theorem with numerical calculation

1 Like

Okay, wrong choice of word. I just want to say that both values come out to be same, which is what I expected. Is there something I’m missing here?

1 Like

Welcome to the julia discourse @Minimum-Pollution-96!

I don’t know exactly what’s going on with your example above but I would guess it’s some kind of mismatch between the conventions being used by the DSP and FFTW packages (I’m assuming you’re using FFTW for your Fourier transforms. One thing I noticed is that if you do

``````t = 0:0.0001:0.5;
fconv_fda = ifft(fft(f1).*fft(f2)) ./ 2
``````

then you get

``````julia> norm(fconv_fda - fconv_tda)/norm(fconv_tda)
0.00021074444869044446
``````

so things are starting to match up well. One handy resource is to start digging into the conventions is: API · AbstractFFTs.jl

Another good read, from the docs of the underlying FFTW C code:

4 Likes

A good explanation is provided in the following lecture, in particular see chapter 4.2.6 Convolution of two finite-duration signals using the DFT

It basically boils down to pad the input signals with enough zeros:

``````using FFTW, DSP, Plots; gr()

Δt = 0.004; # sampling period (s)
t = collect(0:Δt:1.0)[1:end-1]
N = length(t)
f1, f2 = sin.(2π*15*t), sin.(2π*5*t)
f1p, f2p = [f1;zeros(N)], [f2;zeros(N)]
fconv_tda = DSP.conv(f1p, f2p)[1:N]
fconv_fda = ifft(fft(f1p).*fft(f2p))[1:N]
plot(t, fconv_tda,lw=3,lc=:red, xlabel="time /s",ylabel="Amplitude")
plot!(t, real.(fconv_fda),ls=:dash,lc=:black,lw=2)
``````

8 Likes

Sorry, totally off topic here, but @jling I think the people that proved the four color theorem would disagree with you

4 Likes

I know you’re joking, but I believe that four-color proof was more about combinatorics than numerics. I doubt if there was any floating point arithmetic there

2 Likes

that was not a “numerical calculation”.

I didn’t mean “computation done by computer can’t prove theorem”, we even have languages dedicated for writing formal proof.

The whole point of interval arithmetic methods is to do rigorous computations using floating-point arithmetic. Various important theorems have been proved using these methods, including the Kepler conjecture on sphere packings.

12 Likes

Thank you so much! If its not too much, can I ask you regarding the corollary as well?

``````fmult_tda = f2.*f3
# fmult_fda = real.(bfft(conv(fft(f2p),fft(f3p)))/(2N)^2)[1:N]
fmult_fda = real.(ifft(conv(fft(f2p),fft(f3p)))/(N))[1:N]
plot(fmult_tda, size=(1000, 400))
plot!(fmult_fda)
``````

I read in FFTW Website that fft is unnormalised and hence I divided by N. But still, the plots generated is weird and seem to have lower frequency (with normalised amplitude?).

On a side note, AbstractFFT.jl seems to have separate function `bfft` which is unnormalised but it needs me to divide by `(2N)^2` instead. Is AbstractFFT not FFTW’s julia implementation?

Sorry for posting too many questions in single post. Also, thank you for the link. I’m trying to read through it. I tried to find the videos for lecture but those seem unavailable. This is my first time doing any kind of Fourier transform beyond basic math. So, if there’s any further suggestion, kindly let me know

Can you share a link or any available reference for this proof? I am very interesting.

Thank you @Pbellive !

Now, if I consider only [1:N] portion of each after fft, I get the frequencies matched up but the amplitude is different.

``````Δt = 0.001; # sampling period (s)
t = collect(-1:Δt:1.0)[1:end-1]
N = length(t)
f1, f2 = sin.(2π*15*t), sin.(2π*5*t)
f1p, f2p = [f1;zeros(N)], [f2;zeros(N)]
fmult_tda = f1.*f2
fmult_fda = real.(ifft(conv(fft(f1p)[1:N],fft(f2p)[1:N])/0.5N))[1:N]
plot(fmult_tda, size=(1000, 400))
plot!(fmult_fda)
``````

The code and reference links are here:

Here are a couple of references. Unfortunately they’re not exactly light reading:

5 Likes

@Minimum-Pollution-96, for the corollary one needs to take care of the negative frequencies with `fftshift`:

``````# "Corollary":
using FFTW, DSP, Plots; gr()

Δt = 0.004; # sampling period (s)
t = collect(-1:Δt:1.0)[1:end-1]
N = length(t)
f1, f2 = sin.(2π*15*t), sin.(2π*5*t)
fmult_tda = f1 .* f2

F1, F2 = fftshift(fft(f1)), fftshift(fft(f2))
F1p, F2p = [F1;zeros(N)], [F2;zeros(N)]
fmult_fda = real.(ifft(DSP.conv(F1p,F2p)[1:N]))
plot(t,fmult_tda, lw=3,lc=:red, xlabel="time /s",ylabel="Amplitude")
plot!(t,fmult_fda/(0.5N), ls=:dash,lc=:black,lw=2)
``````

1 Like

Thank you!