Update wfs.py

- add WFS prefilter note
- bandlimited Dirac impulse for nicer plots
- finer grid for nicer plots
- use cmap='Blues' for nicer plots
- make point source normalization consistent with frequency domain

Author: fs446

sfs/td/wfs.py

@@ -8,7 +8,7 @@
     import matplotlib.pyplot as plt
     import numpy as np
     import sfs
-    from scipy.signal import unit_impulse
+    from scipy.signal import unit_impulse, firwin
 
     # Plane wave
     npw = sfs.util.direction_vector(np.radians(-45))
@@ -26,19 +26,27 @@
 
     # Impulsive excitation
     fs = 44100
-    signal = unit_impulse(512), fs
+    # either with: signal = unit_impulse(512), fs
+    # or with: bandlimited Dirac, i.e. a linear-phase lowpass FIR
+    Nlp = 2**10 + 1
+    tlp = (Nlp/2) / fs
+    signal = firwin(numtaps=Nlp,
+                    cutoff=16000,
+                    window=('kaiser', 4),
+                    pass_zero='lowpass',
+                    scale=True, fs=fs), fs
 
     # Circular loudspeaker array
     N = 32  # number of loudspeakers
     R = 1.5  # radius
     array = sfs.array.circular(N, R)
 
-    grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02)
+    grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.01)
 
     def plot(d, selection, secondary_source, t=0):
         p = sfs.td.synthesize(d, selection, array, secondary_source, grid=grid,
                               observation_time=t)
-        sfs.plot2d.level(p, grid)
+        sfs.plot2d.level(p, grid, cmap='Blues')
         sfs.plot2d.loudspeakers(array.x, array.n,
                                 selection * array.a, size=0.15)
 
@@ -108,7 +116,8 @@ def plane_25d(x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None):
         delays, weights, selection, secondary_source = \
             sfs.td.wfs.plane_25d(array.x, array.n, npw)
         d = sfs.td.wfs.driving_signals(delays, weights, signal)
-        plot(d, selection, secondary_source)
+        # note that WFS prefilter is not included
+        plot(d, selection, secondary_source, t=tlp)
 
     """
     if c is None:
@@ -187,8 +196,10 @@ def point_25d(x0, n0, xs, xref=[0, 0, 0], c=None):
 
         delays, weights, selection, secondary_source = \
             sfs.td.wfs.point_25d(array.x, array.n, xs)
+        weights *= 4*np.pi*rs  # normalize
         d = sfs.td.wfs.driving_signals(delays, weights, signal)
-        plot(d, selection, secondary_source, t=ts)
+        # note that WFS prefilter is not included
+        plot(d, selection, secondary_source, t=ts+tlp)
 
     """
     if c is None:
@@ -281,8 +292,10 @@ def point_25d_legacy(x0, n0, xs, xref=[0, 0, 0], c=None):
 
         delays, weights, selection, secondary_source = \
             sfs.td.wfs.point_25d(array.x, array.n, xs)
+        weights *= 4*np.pi*rs  # normalize
         d = sfs.td.wfs.driving_signals(delays, weights, signal)
-        plot(d, selection, secondary_source, t=ts)
+        # note that WFS prefilter is not included
+        plot(d, selection, secondary_source, t=ts+tlp)
 
     """
     if c is None:
@@ -363,8 +376,10 @@ def focused_25d(x0, n0, xs, ns, xref=[0, 0, 0], c=None):
 
         delays, weights, selection, secondary_source = \
             sfs.td.wfs.focused_25d(array.x, array.n, xf, nf)
+        weights *= 4*np.pi*rs  # normalize
         d = sfs.td.wfs.driving_signals(delays, weights, signal)
-        plot(d, selection, secondary_source, t=tf)
+        # note that WFS prefilter is not included
+        plot(d, selection, secondary_source, t=tf+tlp)
 
     """
     if c is None: