From fc71ae905f997d5650ece1be804f8f3ba13f718b Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 18 Feb 2021 21:38:22 -0500 Subject: Code for all-pass panning --- panap.bqn | 253 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 253 insertions(+) create mode 100644 panap.bqn diff --git a/panap.bqn b/panap.bqn new file mode 100644 index 0000000..0817387 --- /dev/null +++ b/panap.bqn @@ -0,0 +1,253 @@ +# Stereo panning + +o ← ≠◶⟨•Import∘"options.bqn", ⊑⟩ •args +mix‿fil ← •Import¨ "mix.bqn"‿"filter.bqn" +SP ← (4⊸×÷5⊸-) 4׬⊸× +Sin ← (⊢-∘⊢⍟(>⟜1)·SP 1⊸|) 2|÷⟜π +Cos ← Sin (π÷2)⊸- + +# The goal of stereo panning is to make sounds appear to come from +# different directions. A well-panned sound for music will have the +# following three features: +# +# 1. The high frequencies (and preferably only the high frequencies) are +# louder on the near side. +# 2. The low frequencies begin slightly earlier on the near side. +# 3. It still sounds good when the channels are combined. +# +# Requirements 1 and 2 come from physics: these are properties acquired by +# sounds as they travel around a human head, and they are the features the +# brain uses to identify their source. Requirement 3 comes from the fact +# that music is sometimes played in mono and should still have acceptable +# quality when mixed this way. +# +# Requirement 1 causes no problems and is simple to implement by moving part +# of the far channel to the near side (that is, subtracting a fraction of +# that channel, possibly with a high-pass filter, from the far side and +# adding it to the near side). Requirements 2 and 3, however, are in direct +# conflict. The obvious solution to 2 is to delay the far channel; however +# doing so will cause comb filtering when channels are re-combined. +# +# Typically the answer is to drop one of the requirements. Acceptable stereo +# localization can still be obtained without any delays, thus dropping +# requirement 2. And panning with a delay, or more generally applying a +# head-related transfer function (HRTF), yields the best possible +# localization at the expense of a weak mono mix. +# +# However, it is possible to compromise, obtaining excellent (but not +# perfect) stereo and good mono. The key is to realize that only the lower +# frequencies must be delayed, a feat which can be performed with all-pass +# filters. For delays under 1/2000 of a second we can combine a forwards +# filter with a backwards filter to delay a wide range of frequencies and then +# undelay those over 2000Hz, since those frequencies are not used to detect +# delay. This leaves a signal largely without artifacts--the only issue is +# a small amount of attenuation at frequencies near the cutoff when the two +# channels are combined. + +# Pan with all-pass filters. Arguments are the same as pan +# Delays the input by slightly more than 60 samples. +PanIR ← { + AP ← fil.AllPass + near ← 0.5>𝕨 + off ← | 1-˜2×𝕨 + 𝕩 ↩ (-⊸≍ off × 5000 fil.Lp1 ⊏)⊸+ 𝕩 + a‿b ← Get_panap_coeff Shift_coeff off×25 + near ⌽ (b⊸AP⌾⌽ a⊸AP)⌾⊏ 𝕩 +} +Window ← (2÷˜1+Cos π×1↓↕⊸÷6)⊸(⊣ ⌽⊸mix.Fadefront mix.Fadeback) # Tukey window +PanAP ⇐ Window∘PanIR⍟(≠⟜0.5)⟜(≍˜60=↕180)⊸mix.Reverb + +# Empirically determined good coefficients for Panap +# y is a single parameter giving the delay for low frequencies. +Shift_coeff ← { + f ← 732 + 79000 ÷ 1.6 ⋆˜ 10 + 𝕩 + rs ← 1.097‿¯2.054‿2.16 {+⟜(𝕩⊸×)´𝕨} 1e4 ÷˜ f + 𝕩‿rs‿f +} + +# Data for Shift_coeff +# Values of rs which make the third derivative of the phase transfer at +# zero equal to zero. (rs × 1000) is given. +# frequency of max phase difference +# 2000 1750 1500 1250 1100 1000 +# s 5 782 804 827 852 868 878 +# h 10 816 829 845 864 876 885 +# i 15 849 857 866 878 887 894 +# f 20 875 879 885 893 899 904 +# t 25 895 897 901 906 910 913 +# +# "Good" shift, frequency, and rs values. +# Frequency is determined from the formula in Shift_coeff. +# In Shift_coeff, rs is modelled as quadratic w.r.t frequency. +# ˜0 2716 0.698 +# 0.1 2685 0.701 +# 1 2436 0.725 +# 2.5 2121 0.759 +# 5 1769 0.802 +# 10 1387 0.853 +# 15 1190 0.882 +# 20 1074 0.900 +# 25 1000 0.913 +# 30 948 0.924 + +# ========================================================= +# Argument is ⟨number of samples to shift, magnitude, stop frequency⟩. +# Magnitude should be a real number strictly between 0 and 1, +# and controls how pointy the end result is in phase space. +# Returns two all-pass filter inputs matching the specifications. +# +# The polynomials here should probably be treated as black magic, but +# a derivation from the constraints is included below anyway. +Get_panap_coeff ← { 𝕊 n‿rr‿f: + PM ← +´¨+⌜○↕○≠⊔○⥊×⌜ # Multiply polynomials + + r ← טrr + c1‿c2 ← 2×Cos¨ 1‿2 × 2×π×f÷o.freq + + q ← (r-1)×2÷n + poly ← +´ ⟨ + 4×⟨1+r+q,¯2⟩ PM (⟨r,1+r,0⟩×⟨0,6,0⟩+c2-2×c1) + ⟨1+טr,0,4⟩×1-c1 + ⟨1+r,¯2⟩ PM˜⊸PM (2+c1)×⟨1+r,-c1⟩ + ⟩ + + P ← {+⟜(𝕩⊸×)´poly} + v ← P¨ bnd←0‿rr + ! 0≥×´v + s ← -×⊑v + re ← ⊑ (s×v) {u←s×P m←2÷˜+´𝕩⋄𝕊⍟(1e¯15<-˜´∘⊢)´u‿m((0≤u)⊑⌈‿⌊)¨𝕨‿𝕩} bnd + + sv ← (ט ÷ (4×q) + ×⟜((3-r)+2×re)) 1+r-2×re + + GetZ ← ⊢ ≍ √∘-⟜(ט) # 𝕩 is real part and 𝕨 is square magnitude + (GetZ○⊑ ≍○< GetZ○(⊑+sv×+´))´ ⟨r‿¯1,re‿1⟩ +} + +# --------------------------------------------------------- +" # Derivation +# Transfer function for an all-pass filter: argument is z, +# and rs and re are (ט∘|) and (9&o.) of the complex parameter +H = (1 + (rs×ט) - 2×re&×) ÷ (ט + rs - 2×re&×) +# Derivative with respect to z +Hp = (2 × (re-˜rs&×) - H×(1-re)) ÷ (ט+rs-2×re&×) + = (2 × (× rs-H) - re×1-˜H) ÷ (ט+rs-2×re&×) + +# Definition of the phase transfer function T +(⋆⍳T ω) = (H ⋆⍳ω) +# Derivative in terms of H and its derivative Hp +⍳(⋆⍳T(ω))×Tp(ω) = Hp(⋆⍳ω)×⍳⋆⍳ω +Tp(ω) = Hp(⋆⍳ω) × ⋆⍳ω-T(ω) + = Hp(z) × z ÷ H(z) + = (2×z × (z×rs-H) - re×H-1) ÷ (1 + (rs×טs) - 2×re×z) + +# Simplified at 1 and ¯1 (0 and maximum frequency) +H(1) = 1 +T(0) = 0 +Tp(0) = (2×rs-1) ÷ (1+rs-2×re) + +H(¯1) = 1 +T(π) = 0 +Tp(π) = (2×rs-1) ÷ (1+rs+2×re) + + +# For a delay +H(z) = z⋆-n +T(ω) = -n×ω +Tp(ω) = -n + +# Constraints +# We refer to the forward and backwards filters with +# postfixes of 1 and 2. +Tp1(0) = Tp2(0) - n +Tp1(π) ≃ Tp2(π) +Tp1(f) = Tp2(f) + +# In a symmetric notation: [a,b] = (a1÷b1) - (a2÷b2) +# where an and bn are the values for the nth filter. +# z2 = טz +(-n÷2)= [rs-1 , 1+rs-2×re] +0 = [rs-1 , 1+rs+2×re] +0 = [((z×rs-H) - re×H-1) , (1 + (z2×rs) - (2×z×re))] + +# Useful identity +# Define (reC = re2 - re1) and (rsC = rs2 - rs1) +rs1×re2 - rs2×re1 = (rs1×re1 + rs1×reC) - (rs1×re1 + rsC×re1) + = rs1×reC - rsC×re1 + +# Second constraint, where [12] indicates the value with 1 and 2 swapped +((rs1-1) × 1+rs2+2×re2) = [12] +((rs1 × 1+2×re2) - (rs2+2×re2)) = [12] +(rs1 + re1 + rs1×re2) = [12] +0 = rsC + reC + (re1×rsC - reC×rs1) + = (rsC × re1+1) - (reC × rs1-1) +(reC ÷ re1+1) = (rsC ÷ rs1-1) + +# Thus, we can define s so that +reC = s×re1+1 +rsC = s×rs1-1 + +rs1×re2 - rs2×re1 = rs1×reC - rsC×re1 + = s × (rs1×re1+1) - (re1×rs1-1) + = s × rs1 + re1 + +# Cross-multiply first constraint and subtract second +(rs2-rs1) = (n÷8)×(1+rs1-2×re1)×(1+rs2-2×re2) +(s×rs1-1) = (n÷8)×(1+rs1-2×re1)×((1+rs1-2×re1) + s×(rs1-1)-2×(re1+1)) +(s × (rs1-1)×(8÷n)) = (ט 1+rs1-2×re1) - s×(1+rs1-2×re1)×(rs1-˜3+2×re1) +s = (ט 1+rs1-2×re1) ÷ (((rs1-1)×(8÷n)) + (1+rs1-2×re1)×(rs1-˜3+2×re1)) + + +# Third constraint +# Value of Tp, dropping factor of 2×z +((z×rs1-h1) + re1×h1-1) × (1 + (z2×rs2) - (2×z×re2)) = [12] +(re1 -˜ z×rs1 + h1×re1-z) × (1 + (z2×rs2) - (2×z×re2)) = [12] +((z×rs1) - (re1 + z2×rs1×re2)) + (h1×re1-z)×(1 + (z2×rs2) - (2×z×re2)) = [12] +((z×rs1) - (re1 + z2×rs1×re2)) + (h1×h2×re1-z)×(rs2 + z2 - 2×z×re2) = [12] +((z×rs1) + (z2×rs2×re1) - re1) + h1×h2×((rs2×re1) + (z×rs1) - z2×re1) = [12] +((z×rsC) + (z2×(rs1×re2 - rs2×re1)) - reC) + h1×h2×((rs1×re2 - rs2×re1) + (z×rsC) - z2×reC) = 0 +# Divide by s +((z×rs1-1) + (z2×(rs1+re1)) - re1+1) + h1×h2×((rs1+re1) + (z×rs1-1) - z2×re1+1) = 0 +((rs1×z2+z) + (re1×z2-1) - z+1) + h1×h2×((rs1×z+1) + (re1×1-z2) - z2+z) = 0 +# Divide by z+1 +((rs1×z) + (re1×z-1) - 1) + h1×h2×(rs1 + (re1×1-z) - z) = 0 +h1×h2 = ((rs1×z) + (re1×z-1) - 1) ÷ ((-rs1) + (re1×z-1) + z) + = 1 + (z+1)×(1-rs1) ÷ (rs1 + (re1×1-z) - z) + +# Value of h2 in terms of re1 and rs1 +# With these and the previous definition, we obtain a 4th-order +# polynomial in re1 and rs1. +a = 2×z×re1, r = rs1 +h1 = ((1 + (z2×r) - a) ) ÷ ((z2 + r - a) ) +h2 = ((1 + (z2×r) - a) + s×(((z2×r) - a) - (z×z+2))) ÷ ((z2 + r - a) + s×((r - a) - (1+2×z))) + +# Algebra +h2 = 1 + ((r-1)×(z2-1)×(1+s)) ÷ ((-a×1+s) + (z2 + r) + s×(r - (1+2×z))) + = 1 + (r-1)×(z2-1) ÷ (((ט1+z)÷1+s) + (r - (1+2×z)) - a) + = 1 + (r-1)×(z2-1) ÷ ((z2+r-a) - (ט1+z)×s÷1+s) +h1 = 1 + (r-1)×(z2-1) ÷ (z2+r-a) +s÷1+s = (ט1+r-2×e) ÷ (4 × ((r-1)×(2÷n)) + 1+r-2×e) + = (ט1+r-2×e) ÷ (4 × q + 1+r-2×e) + +# Put each of our values in a regular form, defining k, a=b+d, b, c +h2 = (1+k÷a) = 1 + (1-r)×(1-z2) ÷ ((z2+r-2×z×e) - (ט1+z)×s÷1+s) +h1 = (1+k÷b) = 1 + (1-r)×(1-z2) ÷ (z2+r-2×z×e) +h1×h2 = (1+k÷c) = 1 + (1-r)×(1-z2) ÷ (1-z)×((r-z) + (1-z)×e) + +# Simplify the overall form +(1+k÷a)×(1+k÷b) = (1+k÷c) +c×(a+k)×(b+k) = ab×(c+k) # Cross-multiply +c×(a+b+k) = ab # Subtract abc +c×k+b = a×b-c + +c×k+b = (1-z)×((r-z)+(1-z)×e) × ((1-r)×(1-z2))+z2+r-2×z×e + = (1-z)×((r-z)+(1-z)×e) × (1+z2×r)-2×z×e +a×b-c = ((z2+r-2×z×e) - (ט1+z)×s÷1+s)×((z2+r-2×z×e)-(1-z)×((r-z) + (1-z)×e)) + = ((z2+r-2×z×e) - (ט1+z)×s÷1+s)×((z×1+r) - e×1+z2) + +q←(r-1)×(2÷n) +((1-z)×⟨r-z,1-z⟩PM⟨1+z2×r,-2×z⟩PM⟨4×q+1+r,¯8⟩) + ⟨z×1+r,-1+z2⟩ PM (PM˜(1+z)×⟨1+r,¯2⟩) - ⟨z2+r,-2×z⟩PM⟨4×q+1+r,¯8⟩ +# Group q +{(4×⟨q+1+r,¯2⟩PM(𝕩PM⟨z2+r,-2×z⟩)-˜(1-z)×⟨r-z,1-z⟩PM⟨1+z2×r,-2×z⟩) + 𝕩 PM PM˜(1+z)×⟨1+r,¯2⟩}⟨z×1+r,-1+z2⟩ + +# Expand; divide by z2; simplify +zd←z+÷z ⋄ ze←z2+÷z2 +(4×⟨q+1+r,¯2⟩PM(⟨r,1+r,0⟩×⟨0,6,0⟩+ze-2×zd)+⟨1+טr,0,4⟩×1-zd) + (PM˜⟨1+r,¯2⟩)PM(2+zd)×⟨1+r,-zd⟩" -- cgit v1.2.3