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
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?
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:
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)
Sorry, totally off topic here, but @jling I think the people that proved the four color theorem would disagree with you 
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 
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.
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:
https://annals.math.princeton.edu/wp-content/uploads/annals-v162-n3-p01.pdf
@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)
Thank you!