Hi there, sorry for the late reply. I get slightly confused when using DSP.jl <https://github.com/JuliaDSP/DSP.jl>. I am trying to implement the inverse filter technique to measure linear systems impulse response. So, I define a filter with DSP.jl and the test signal x (a log swept sine), with sample rate fs. Then, I calculate the inverse filter z (which is an operation on x, pretty much a time reversal with an exponential amplitude envelope). If y is the output of the linear filter the impulse response is given by the convolution of y with z. So, I am running this algorithm and comparing with the results from impz in order to understand whether it works correctly. My results:
- Manual Convolution: Y = (1 / fs) * fft(y) Z = (1 / fs) * fft(z) H = Y .* Z h = fs * ifft(H) - Convolution with conv(): h = (1 / fs) * conv(y, z) - Impz: h = impz(filter, t) where t is the same time array resulting from the other operations. What happens is that the first two operations agree, while the results from impz need to be multiplied by fs to agree with them. However, calling freqz produces good agreement: H = freqz(filter, f, fs) I am tempted to conclude that freqz makes the normalization (knowing fs), while impz doesn't. But I was also wondering whether my normalizations are wrong. By the way, if I calculate h with conv() and then H as (1 / fs) * fft(h) the magnitudes are ok, but the phase of H is completely different from the expected one. In this case, the manual convolution works way better. Do you have an idea why this is happening? (maybe I will open a separate thread for this) By the way, I am following this <http://ant-novak.com/swept-sine.php> to define test signal and inverse filter. On Saturday, 19 March 2016 10:01:17 UTC, Jeffrey Sarnoff wrote: > > 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)) >> ... > >
