I am trying to understand scaling of the irfft function of the FFTW.jl package. I have a function that defines a continuous frequency domain signal (in radians). The complex portion of the signal has odd symmetry and the real part of the signal has even symmetry. I sample the positive half of the frequency domain signal and end up with an array of N elements. I then apply the irfft() function to the signal. I get the right shape out, but it is not scaled properly.
Can anyone advise me on the scaling factors that are applicable to this function?
Would it also be necessary to convert frequency domain signal to Hertz rather than radians before applying the irfft()?
It performs the “backward” transform defined here multiplied by a scale factor of 1/n where n is the length of the real output array, so that it is the inverse of rfft.
No. The frequencies are not an input to a DFT, and it is a linear operation so that it doesn’t matter what scaling/units you use.
Ok, so does the sampling interval in the frequency domain affect the scaling? In the example below, I have a fixed number of samples, N. When I change ωmax (the maximum frequency in the range), the amplitude of the ifft result changes even though N does not change. The plots below are for ωmax=4 and ωmax=20 respectively. The time domain signal amplitude seems inversely related to the frequency sampling interval, despite N remaining constant. To recover the correct time domain signal amplitude it seems I need to include some other scaling factor related to frequency.
using FFTW
using SpecialFunctions
using Gaston #plotting package
#function to generate a frequency domain signal
function K0diff(ω,r,r′,z)
if ω == 0
result = log(r′/r) + 0im
else
result = exp(1im*ω*z)*(besselk(0,abs(ω)*r) - besselk(0,abs(ω)*r′))
end
return result
end
#create a frequency domain function and apply irfft
N = 2^10
Ndiv2 = Int(floor(N/2))
ωmax = 4
ω = range(-ωmax,ωmax,length = N + 1)
frequencyDomainSignal = K0diff.(ω,5,20,1)
p = plan_irfft(frequencyDomainSignal[Ndiv2+1:end],N,1)
timeDomainSignal = fftshift(p*frequencyDomainSignal[Ndiv2+1:end],1)
#output plots using Gaston
p1 = plot(real.(frequencyDomainSignal),Axes(xlabel = "'Frequency Index'",ylabel = "'Amplitude'"), handle = 1)
plot!(imag.(frequencyDomainSignal))
p2 = plot(timeDomainSignal,Axes(xlabel = "'Time Index'",ylabel = "'Amplitude'"), handle = 2)
plot([p1 ; p2])
If you are thinking of the DFT as an approximation for a continuous Fourier transform, then you should multiply it by Δω in order for the sum to be an approximation for \int d\omega. (The finite window is also an approximation.)
This has nothing to do with Julia or the FFTW.jl package. Any FFT function computes the discrete Fourier transform (DFT), not the continuous Fourier transform, and you need to understand how the two relate for your application.
Ok, that makes sense and fixed my problem. DSP related functions are outside my comfort zone and there is so much information available. I was barking up the wrong tree in my research. Thanks for the redirect.
