Fix density calculation

giuliovv · Oct 31, 2021 · b555663 · b555663
1 parent 57b555c
commit b555663
Show file tree

Hide file tree

Showing 5 changed files with 65 additions and 55 deletions.
diff --git a/examples/test.ipynb b/examples/test.ipynb
diff --git a/examples/tfest_example.ipynb b/examples/tfest_example.ipynb
diff --git a/setup.py b/setup.py
@@ -4,15 +4,15 @@
 setuptools.setup(
   name = 'tfest',
   packages = ['tfest'],
-  version = '0.1.11',
+  version = '0.1.12',
   license="""MIT""",
   description = """Transfer function estimation based on frequency response.""",
   long_description_content_type="text/markdown",
   long_description=README,
   author = 'Giulio Vaccari',
   author_email = 'io@giuliovaccari.it',
   url = 'https://github.com/giuliovv/tfest',
-  download_url = 'https://github.com/giuliovv/tfest/archive/refs/tags/v0.1.11-alpha.tar.gz',
+  download_url = 'https://github.com/giuliovv/tfest/archive/refs/tags/v0.1.12-alpha.tar.gz',
   keywords = ['tfest', 'frequency', 'matlab'],
   install_requires=[
           'matplotlib',

diff --git a/tfest/__pycache__/tfest.cpython-36.pyc b/tfest/__pycache__/tfest.cpython-36.pyc
diff --git a/tfest/tfest.py b/tfest/tfest.py
@@ -44,29 +44,30 @@ def frequency_response(self, method="h1", time=None):
         if time == None:
             warnings.warn("Setting default time=1")
             time = 1
+        dt = time/len(self.u)
         if method == "fft":
             u_f = fft(self.u)
             y_f = fft(self.y)
             u_no_zero = u_f[np.nonzero(u_f)]
             y_no_zero = y_f[np.nonzero(u_f)]
-            frequency = fftfreq(u_f.size, d=time/len(self.u))[np.nonzero(u_f)]
+            frequency = fftfreq(u_f.size, d=dt)[np.nonzero(u_f)]
             H = y_no_zero/u_no_zero
         elif method == "h1":
             # https://dsp.stackexchange.com/questions/71811/understanding-the-h1-and-h2-estimators
-            cross_sd, frequency = csd(self.y, self.u, Fs=time/len(self.u))
-            power_sd, _ = psd(self.u, Fs=time/len(self.u))
+            cross_sd, frequency = csd(self.u, self.y, Fs=1/dt, NFFT=len(self.u))
+            power_sd, _ = psd(self.u, Fs=1/dt, NFFT=len(self.u))
             H = cross_sd/power_sd
         elif method == "h2":
-            cross_sd, frequency = csd(self.u, self.y, Fs=time/len(self.u))
-            power_sd, _ = psd(self.y, Fs=time/len(self.u))
+            cross_sd, frequency = csd(self.y, self.u, Fs=1/dt, NFFT=len(self.u))
+            power_sd, _ = psd(self.y, Fs=1/dt, NFFT=len(self.u))
             H = power_sd/cross_sd
         else:
             raise Exception("Unknown method")
         self.frequency = frequency
         self.H = H
         return self.H, frequency
 
-    def estimate(self, nzeros, npoles, init_value=1, options={'xatol': 1e-2, 'disp': True}, method="density", time=1):
+    def estimate(self, nzeros, npoles, init_value=1, options={'xatol': 1e-3, 'disp': True}, method="h1", time=None):
         """
         npoles: number of poles
         nzeros: number of zeros