Hi there cool people!
I am making some Julia experiments around fft, conv and normalization. Long
story short, by using signals for which it is simple to calculate
analytically transforms and convolutions, I just play with normalization
until I find match of the Julia calculated quantities with the analytical
ones. Here what I found so far:
- .fft output needs to be normalized by (1/fs), with fs the sample rate
of the input signal (ifft will be then normalized by fs)
- conv output needs to be normalized by (1/fs) as well
The fft normalization is not surprising (I was used to the same in Matlab).
What do you think? It seems to be working ok but I also feel like I am
missing something... Here the playground scripts:
# convplay.jl
#=
Make some experiment with known convolutions:
=#
using PyPlot
# Sampling Variables:
Fₛ = 96000 # Sample Rate
# Time windows of the signals:
T₁ = -20 # seconds
T₂ = 20 # seconds
# Define two rectangular pulses:
α = 2 # Amplitude
τ₁ = 4 # seconds
τ₂ = 6 # seconds
a = 1 # Amplitude
t₁ = 0 # seconds
t₂ = 4 # seconds
S₁ = ceil(T₁ * Fₛ) # Everything before the first sample is lost.
S₂ = floor(T₂ * Fₛ) # Everything after the last sample is lost.
t = (1 / Fₛ) * (S₁:S₂)
# Initialize and populate the signals:
x₁ = zeros(t)
x₂ = zeros(t)
for s = 1:length(t)
if t[s] >= τ₁ && t[s] <= τ₂
x₁[s] = α
end
if t[s] >= t₁ && t[s] <= t₂
x₂[s] = a
end
end
# Convolve:
c = (1 / Fₛ) * conv(x₁, x₂) # 4, the correct peak of c, is found with this
normalization.
# fftplay.jl
#=
Have some PHUN with fft.
=#
using PyPlot
# Close existing plot windows:
close("all")
# Define Signal stuff:
Fₛ = 96000 # Sample Rate. Hz.
P = 1 # Amplitude
t₁ = -0.5 # Initial instant of signal time window (with pads). s.
t₂ = 0.5 # Final instant of signal time window (with pads). s.
T = 3e-3 # Width of window containing the actual signal (N pulse). s.
Centered on 0 s.
# Define Processing Stuff:
Mprime = 10 # https://groups.google.com/forum/#!topic/julia-dev/MdTBxgVoab0
# ^^ fft is more efficient if the signal is operating on has a length that
can
# be expressed as a product of small prime numbers:
# p1^e1 * p2^e2 * ... ei integer exponents, pi prime numbes, i = 1,2,...
# Small: smaller than ten.
# Define the N pulse signal:
Δt = t₂ - t₁ # Global Signal Window. s.
Δtₛ = floor(Δt * Fₛ) # Same in samples.
t = (1 / Fₛ) * ((1 : Δtₛ) - floor(Δtₛ / 2)) # Time. Keep as range and don't
waste memory. s.
# Signal ranges in samples:
_ , T₁ = findmin(abs(t - (-T / 2)))
_ , T₂ = findmin(abs(t - (+T / 2)))
# Define Signal:
x = zeros(t)
x[T₁ : T₂] = (-2 * P / T) * t[T₁ : T₂]
# Plot the signal:
figure(1)
plot(t, x)
# Prepare for transform:
nfft = nextprod(primes(Mprime), Δtₛ) # this is pretty much how to optimize
fourier transform length.
x_pad = [x; zeros(nfft - length(x))]
# We are now ready for fft stuff:
# Set up multithread:
FFTW.set_num_threads(2) # I got two cores.
X = (1 / Fₛ) * fft(x_pad)
# Define Frequency Stuff:
f = linspace(0, Fₛ, nfft + 1) # A snake biting its tail...
f = f[1 : nfft]
# Nyquist Sample:
if mod(nfft, 2) == 0
nyq = convert(Int, floor(nfft/ 2))
else
nyq = convert(Int, floor(nfft/ 2) + 1)
end
# Positive frequencies and transforms:
f₊ = f[1 : nyq]
# Multiply for sqrt(2) so that the integral of the square magnitude over
the positive frequency
# gives the correct total power:
X₊ = sqrt(2) * X[1 : nyq]
figure(2)
semilogx(f₊, abs(X₊))
xlim(xmax = 20e3)
ylim(ymax = 2e-3)
# For real signals the positive frequency trasform can be calculated
directly:
X₊ᵈ = (sqrt(2) / Fₛ) * rfft(x_pad)
# This will contain the very last frequency sample (Nyquist frequency) as
well:
f₊ᵈ = [f₊; Fₛ/2]
# Check if any differences:
isdif_transf = isequal(X₊, X₊ᵈ[1 : nyq])
# Try to reconstruct the original signal:
x₁ = Fₛ * ifft(X) # This has an unexpected very small imaginary part.
# A good way to get rid of the imaginary part appears to be irfft, that
works assuming
# exact Hermitian symmetry. So, use only the positive frequency transform,
However, you must
# include all the info (so, up to Nyquist frequency included):
x₂ = Fₛ * irfft(X[1 : (nyq + 1)], length(x))
# Similarly:
x₃ = (Fₛ / sqrt(2)) * irfft(X₊ᵈ, length(x))
# Check if any differences:
isdif_sig_1 = isequal(x, x₂)
isdif_sig_2 = isequal(x, x₃)
# It is to be expected that all the isdif variables will be false.
# Update the figures with the new plots:
figure(1)
plot(t, x₂)
plot(t, x₃)
xlabel(L"$t$ s", fontsize = 20)
ylabel(L"$x\left(t\right)$", fontsize = 20)
xlim(xmin = -T - T/3, xmax = T + T/3)
ylim(ymin = -P - P/4, ymax = P + P/4)
grid(b = true)
figure(2)
plot(f₊, abs(X₊ᵈ[1 : nyq]))
xlabel(L"$f$ Hz", fontsize = 20)
ylabel(L"\left|X\left(f\right)\right|", fontsize = 20)
grid(b = true)
