What do you feel is missing?
On Monday, March 14, 2016 at 6:20:39 PM UTC-4, CrocoDuck O'Ducks wrote:
>
> 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))
> ...
