WIP: filtered source

python/meep.i

@@ -1555,6 +1555,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
     from .source import (
         ContinuousSource,
         CustomSource,
+        FilteredSource,
         EigenModeSource,
         GaussianSource,
         Source,

python/source.py

@@ -2,6 +2,9 @@
 
 import meep as mp
 from meep.geom import Vector3, check_nonnegative
+from scipy import signal
+import numpy as np
+from scipy.special import erf
 
 
 def check_positive(prop, val):
@@ -65,6 +68,7 @@ def __init__(self, frequency=None, width=0, fwidth=float('inf'), start_time=0, c
         self.swigobj = mp.gaussian_src_time(self.frequency, self.width, self.start_time,
                                             self.start_time + 2 * self.width * self.cutoff)
         self.swigobj.is_integrated = self.is_integrated
+        self.peak_time = 0.5 * (self.start_time + self.start_time + 2 * self.width * self.cutoff)
 
     def fourier_transform(self, freq):
         return self.swigobj.fourier_transform(freq)
@@ -80,6 +84,90 @@ def __init__(self, src_func, start_time=-1.0e20, end_time=1.0e20, center_frequen
         self.swigobj = mp.custom_src_time(src_func, start_time, end_time, center_frequency)
         self.swigobj.is_integrated = self.is_integrated
 
+class FilteredSource(CustomSource):
+    def __init__(self,center_frequency,frequencies,frequency_response,dt,T,time_src):
+        dt = dt/2 # divide by two to compensate for staggered E,H time interval
+        self.dt = dt
+        self.frequencies=frequencies
+        self.N = np.round(T/dt)
+        f = self.func()
+
+        # calculate dtft of input signal
+        signal_t = np.array([time_src.swigobj.current(t,dt) for t in np.arange(0,T,dt)]) # time domain signal
+        signal_dtft = np.exp(1j*2*np.pi*frequencies[:,np.newaxis]*np.arange(0,signal_t.size)[np.newaxis,:]*dt)@signal_t # vectorize dtft for speed        
+
+        # multiply sampled dft of input signal with filter transfer function
+        H = signal_dtft * frequency_response
+
+        # estimate the impulse response using a sinc function RBN
+        self.nodes, self.err = self.estimate_impulse_response(H)
+
+        # initialize super
+        super(FilteredSource, self).__init__(end_time=T,src_func=f,center_frequency=center_frequency,is_integrated=False)
+
+    def cos_window_td(self,a,n,f0):
+        omega0 = 2*np.pi*f0
+        cos_sum = np.sum([(-1)**k * a[k] * np.cos(2*np.pi*n*k/self.N )for k in range(len(a))])
+        return np.where(n < 0 or n > self.N,0,
+            np.exp(-1j*omega0*n*self.dt) * (cos_sum)
+            )
+
+    def cos_window_fd(self,a,f,f0):
+        omega = 2*np.pi*f
+        omega0 = 2*np.pi*f0
+        df = 1/(self.N*self.dt)
+        cos_sum = np.zeros((f*f0).shape,dtype=np.complex128)#a[0]*self.sinc(f,f0)
+        for k in range(len(a)):
+            cos_sum += (-1)**k * a[k]/2 * self.sinc(f,f0-k*df) + (-1)**k * a[k]/2 * self.sinc(f,f0+k*df)
+        return cos_sum
+
+    def sinc(self,f,f0):
+        omega = 2*np.pi*f
+        omega0 = 2*np.pi*f0
+        num = np.where(omega == omega0, self.N, (1-np.exp(1j*(self.N+1)*(omega-omega0)*self.dt)))
+        den = np.where(omega == omega0, 1, (1-np.exp(1j*(omega-omega0)*self.dt)))
+        return num/den
+
+    def rect(self,t,f0):
+        n = int(np.round((t)/self.dt))
+        omega0 = 2*np.pi*f0
+        return np.where(n < 0 or n > self.N,0,np.exp(-1j*omega0*(n)*self.dt))
+
+    def hann(self,t,f0):
+        n = int(np.round((t)/self.dt))
+        a = [0.5, 0.5]
+        return self.cos_window_td(a,n,f0)
+
+    def hann_dtft(self,f,f0):
+        a = [0.5, 0.5]
+        return self.cos_window_fd(a,f,f0)
+
+    def nuttall(self,t,f0):
+        n = int(np.round((t)/self.dt))
+        a = [0.355768, 0.4873960, 0.144232, 0.012604]
+        return self.cos_window_td(a,n,f0)
+
+    def nutall_dtft(self,f,f0):
+        a = [0.355768, 0.4873960, 0.144232, 0.012604]
+        return self.cos_window_fd(a,f,f0)
+
+    def __call__(self,t):
+        vec = self.nuttall(t,self.frequencies)
+        return np.dot(vec,self.nodes)
+
+    def func(self):
+        def _f(t): 
+            return self(t)
+        return _f
+
+    def estimate_impulse_response(self,H):
+        # Use vandermonde matrix to calculate weights of each gaussian.
+        # Each sinc is centered at each frequency point
+        vandermonde = self.nutall_dtft(self.frequencies[:,np.newaxis],self.frequencies[np.newaxis,:])
+        nodes = np.matmul(np.linalg.pinv(vandermonde),H)
+        H_hat = np.matmul(vandermonde,nodes)
+        l2_err = np.sum(np.abs(H-H_hat)**2/np.abs(H)**2)
+        return nodes, l2_err
 
 class EigenModeSource(Source):
 

tests/temp.py

@@ -0,0 +1,106 @@
+import meep as mp
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import signal
+import meep_adjoint as mpa
+
+# --------------------------- #
+# Design filter
+# --------------------------- #
+
+freq_min = 1/1.7
+freq_max = 1/1.4
+nfreq = 100
+freqs = np.linspace(freq_min,freq_max,nfreq)
+resolution = 10
+run_time = 500
+
+# --------------------------- #
+# Run normal simulation
+# --------------------------- #
+
+cell = [16,8,0]
+geometry = [mp.Block(mp.Vector3(mp.inf,1,mp.inf),
+                     center=mp.Vector3(),
+                     material=mp.Medium(epsilon=12))]
+fcen = 1/1.55
+gauss_src = mp.GaussianSource(frequency=fcen,fwidth=0.1/1.55,is_integrated=False)
+sources = [mp.EigenModeSource(gauss_src,
+                     eig_band=1,
+                     size=[0,8],
+                     center=[-6,0])]
+pml_layers = [mp.PML(1.0)]
+
+sim = mp.Simulation(cell_size=cell,
+                    boundary_layers=pml_layers,
+                    geometry=geometry,
+                    sources=sources,
+                    resolution=resolution)
+mon = sim.add_dft_fields([mp.Ez],freq_min,freq_max,nfreq,center=[0,0],size=[0,1])
+
+sim.run(until_after_sources=run_time)
+
+fields = np.zeros((nfreq,),dtype=np.complex128)
+# just store one spatial point for each freq
+for f in range(nfreq):
+    fields[f] = sim.get_dft_array(mon,mp.Ez,f)[10]
+
+# build simple bandpass filter
+dt = sim.fields.dt
+fs = 1/dt
+freqs_scipy = freqs * np.pi / (fs/2)
+num_taps = 321
+taps = signal.firwin(num_taps,[freq_min/(fs/2), freq_max/(fs/2)], pass_zero='bandstop',window='boxcar')
+
+w,h = signal.freqz(taps,worN=freqs_scipy)
+
+# frequency domain calculation
+desired_fields = h * fields
+
+# --------------------------- #
+# Run filtered simulation
+# --------------------------- #
+T = 1200
+filtered_src = mp.FilteredSource(fcen,freqs,h,dt,T,gauss_src)
+
+sources = [mp.EigenModeSource(filtered_src,
+                     eig_band=1,
+                     size=[0,8],
+                     center=[-6,0])]
+
+sim.reset_meep()
+sim.change_sources(sources)
+
+mon = sim.add_dft_fields([mp.Ez],freq_min,freq_max,nfreq,center=[0,0],size=[0,1])
+sim.run(until_after_sources=run_time)
+
+fields_filtered = np.zeros((nfreq,),dtype=np.complex128)
+# just store one spatial point for each freq
+for f in range(nfreq):
+    fields_filtered[f] = sim.get_dft_array(mon,mp.Ez,f)[10]
+
+# --------------------------- #
+# Compare results
+# --------------------------- #
+
+plt.figure()
+plt.subplot(2,1,1)
+plt.semilogy(freqs,np.abs(desired_fields),label='Frequency Domain')
+plt.semilogy(freqs,np.abs(fields_filtered),'-.',label='Time Domain')
+plt.grid()
+plt.xlabel('Frequency')
+plt.ylabel('Magnitude')
+plt.legend()
+
+plt.subplot(2,1,2)
+plt.plot(freqs,np.unwrap(np.angle(desired_fields)),label='Frequency Domain')
+plt.plot(freqs,np.unwrap(np.angle(fields_filtered)),'-.',label='Time Domain')
+plt.grid()
+plt.xlabel('Frequency')
+plt.ylabel('Angle')
+plt.legend()
+
+plt.tight_layout()
+plt.savefig('filtered_{}.png'.format(resolution))
+
+plt.show()