diff --git a/src/nmreval/models/bds.py b/src/nmreval/models/bds.py
index 7e770e7..9ccb922 100644
--- a/src/nmreval/models/bds.py
+++ b/src/nmreval/models/bds.py
@@ -10,14 +10,20 @@ class _AbstractBDS:
     equation = ''
     params = [r'\Delta\epsilon', r'\tau_{0}']
     bounds = [(0, None), (0, None)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
     susceptibility = None
     iscomplex = True
 
     @classmethod
-    def func(cls, x, *args, complex_mode: int = 0, **kwargs):
+    def func(cls, x, *args, complex_mode: int = 0, xaxis: str = 'freq', **kwargs):
         # args[0] : Delta epsilon
         # args[1:] : every other parameter
-        chi = args[0] * cls.susceptibility(2*np.pi*x, *args[1:], **kwargs)
+
+        _w = _convert_x_to_omega(x, xaxis=xaxis)
+
+        chi = args[0] * cls.susceptibility(_w, *args[1:], **kwargs)
         if complex_mode == 0:
             return chi
         elif complex_mode == 1:
@@ -64,11 +70,16 @@ class HavriliakNegamiAlphaGammaBDS:
     equation = r'\Delta\epsilon / [1-(i\omega\tau)^{\gamma}]^{\alpha}'
     params = [r'\Delta\epsilon', r'\tau_{0}', r'\alpha', r'\alpha\gamma']
     bounds = [(0, None), (0, None), (0, 1), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
     iscomplex = True
 
     @staticmethod
-    def func(x, deps, tau, alpha, alphagamma, complex_mode: int = 0, **kwargs):
-        chi = deps * HavriliakNegami.susceptibility(2*np.pi*x, tau, alpha, alphagamma/alpha, **kwargs)
+    def func(x, deps, tau, alpha, alphagamma, complex_mode: int = 0, xaxis: str = 'freq', **kwargs):
+        _w = _convert_x_to_omega(x, xaxis=xaxis)
+
+        chi = deps * HavriliakNegami.susceptibility(_w, tau, alpha, alphagamma/alpha, **kwargs)
         if complex_mode == 0:
             return chi
         elif complex_mode == 1:
@@ -115,13 +126,17 @@ class HNWithHF:
     equation = r'\Delta\epsilon HN(\omega, \tau, \alpha, \gamma) / CD(\omega, \tau_{c}, \alpha\gamma-\delta)'
     params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma', r'\tau_{c}', '\delta']
     bounds = [(0, None), (0, None), (0, 1), (0, 1), (0, None), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
     iscomplex = True
 
     @staticmethod
-    def func(x, deps, tau, alpha, gamma, tauc, delta, complex_mode: int = 0):
-        w = 2 * np.pi * x
-        numerator = (1 + 1j*w*tauc) ** (gamma-delta)
-        denominator = (1 + (1j*w*tau)**alpha) ** gamma
+    def func(x, deps, tau, alpha, gamma, tauc, delta, complex_mode: int = 0, xaxis: str = 'freq'):
+        _w = _convert_x_to_omega(x, xaxis=xaxis)
+
+        numerator = (1 + 1j*_w*tauc) ** (gamma-delta)
+        denominator = (1 + (1j*_w*tau)**alpha) ** gamma
 
         epsilon = deps * np.conjugate(numerator / denominator)
 
@@ -144,8 +159,8 @@ class _CCWithHF:
     iscomplex = True
 
     @staticmethod
-    def func(x, deps, tau, alpha, tauc, delta, complex_mode: int = 0):
-        return HNWithHF.func(x, deps, tau, alpha, 1, tauc, delta, complex_mode=complex_mode)
+    def func(x, deps, tau, alpha, tauc, delta, complex_mode: int = 0, xaxis: str = 'freq'):
+        return HNWithHF.func(x, deps, tau, alpha, 1, tauc, delta, complex_mode=complex_mode, xaxis=xaxis)
 
 
 class _CDWithHF:
@@ -157,8 +172,8 @@ class _CDWithHF:
     iscomplex = True
 
     @staticmethod
-    def func(x, deps, tau, gamma, tauc, delta, complex_mode: int = 0):
-        return HNWithHF.func(x, deps, tau, 1, gamma, tauc, delta, complex_mode=complex_mode)
+    def func(x, deps, tau, gamma, tauc, delta, complex_mode: int = 0, xaxis: str = 'freq'):
+        return HNWithHF.func(x, deps, tau, 1, gamma, tauc, delta, complex_mode=complex_mode, xaxis=xaxis)
 
 
 class PowerLawBDS:
@@ -167,16 +182,21 @@ class PowerLawBDS:
     equation = r'A / (i\omega)^{n}'
     params = ['A', 'n']
     bounds = [(None, None), (None, None)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
     iscomplex = True
 
     @staticmethod
-    def func(x, a, n, complex_mode: int = 0):
+    def func(x, a, n, complex_mode: int = 0, xaxis: str = 'freq'):
+        _w = _convert_x_to_omega(x, xaxis=xaxis)
+
         if complex_mode == 0:
-            ret_val = np.exp(1j*n*np.pi/2) * a / x**n
+            ret_val = np.exp(1j*n*np.pi/2) * a / _w**n
         elif complex_mode == 1:
-            ret_val = np.cos(n*np.pi/2) * a / x**n
+            ret_val = np.cos(n*np.pi/2) * a / _w**n
         elif complex_mode == 2:
-            ret_val = np.sin(n*np.pi/2) * a / x**n
+            ret_val = np.sin(n*np.pi/2) * a / _w**n
         else:
             raise ValueError(f'{complex_mode!r} is not 0, 1, 2')
 
@@ -189,17 +209,22 @@ class DCCondBDS:
     equation = r'i\sigma_{dc}/\epsilon_{0}\omega'
     params = [r'\sigma_{dc}']
     bounds = [(0, None)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
     iscomplex = True
 
     @staticmethod
-    def func(x, sigma, complex_mode: int = 0):
+    def func(x, sigma, complex_mode: int = 0, xaxis: str = 'freq'):
+        _w = _convert_x_to_omega(x, xaxis=xaxis)
+
         if complex_mode == 0:
             ret_val = np.zeros(x.shape, dtype=complex)
-            ret_val += 1j * sigma / x / epsilon0
+            ret_val += 1j * sigma / _w / epsilon0
         elif complex_mode == 1:
             ret_val = np.zeros(x.shape)
         elif complex_mode == 2:
-            ret_val = sigma / x / epsilon0
+            ret_val = sigma / _w / epsilon0
         else:
             raise ValueError(f'{complex_mode!r} is not 0, 1, 2')
 
@@ -210,11 +235,16 @@ class DerivativeHavriliakNegami:
     name = 'Derivative HN'
     type = 'Dielectric Spectroscopy'
     params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma']
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
     bounds = [(0, None), (0, None), (0, 1), (0, 1)]
 
     @staticmethod
-    def func(x, eps, tau, a, g):
-        omtau = 2*np.pi*x * tau
+    def func(x, eps, tau, a, g, xaxis: str = 'freq'):
+        w = _convert_x_to_omega(x, xaxis=xaxis)
+
+        omtau = w * tau
         theta = np.arctan2(np.sin(np.pi*a/2.), omtau**(-a) + np.cos(np.pi*a/2.))
 
         numer = a*g * omtau**a * np.cos(np.pi*a/2. - (1+g)*theta)
@@ -228,10 +258,15 @@ class DerivativeColeCole:
     type = 'Dielectric Spectroscopy'
     params = [r'\Delta\epsilon', r'\tau', r'\alpha']
     bounds = [(0, None), (0, None), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
 
     @staticmethod
-    def func(x, eps, tau, alpha):
-        omtau = 2*np.pi*x * tau
+    def func(x, eps, tau, alpha, xaxis: str = 'freq'):
+        w = _convert_x_to_omega(x, xaxis=xaxis)
+
+        omtau = w * tau
         theta = np.arctan2(np.sin(np.pi*alpha/2.), omtau**(-alpha) + np.cos(np.pi*alpha/2.))
 
         numer = alpha * omtau**alpha * np.cos(np.pi*alpha/2. - 2*theta)
@@ -245,10 +280,14 @@ class DerivativeColeDavidson:
     type = 'Dielectric Spectroscopy'
     params = [r'\Delta\epsilon', r'\tau', r'\gamma']
     bounds = [(0, None), (0, None), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
 
     @staticmethod
-    def func(x, eps, tau, g):
-        omtau = 2*np.pi*x * tau
+    def func(x, eps, tau, g, xaxis: str = 'freq'):
+        w = _convert_x_to_omega(x, xaxis=xaxis)
+        omtau = w * tau
         numer = g * omtau * np.sin((1+g)*np.arctan(omtau))
         denom = (1 + omtau**2)**((1+g)/2.)
 
@@ -260,10 +299,13 @@ class _DerivativeHNWithHF:
     type = 'Dielectric Spectroscopy'
     params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma', r'\tau_{c}', r'\delta']
     bounds = [(0, None), (0, None), (0, 1), (0, 1), (0, None), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
 
     @staticmethod
-    def func(x, deps, tau, alpha, gamma, tauc, delta):
-        w = 2*np.pi*x
+    def func(x, deps, tau, alpha, gamma, tauc, delta, xaxis: str = 'freq'):
+        w = _convert_x_to_omega(x, xaxis=xaxis)
 
         a_pi2 = alpha*np.pi/2
         w_lamb = (w*tau)**alpha
@@ -287,10 +329,13 @@ class DerivativeCCWithHF:
     type = 'Dielectric Spectroscopy'
     params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\tau_{c}', r'\delta']
     bounds = [(0, None), (0, None), (0, 1), (0, None), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
 
     @staticmethod
-    def func(x, deps, tau, alpha, tauc, delta):
-        return _DerivativeHNWithHF.func(x, deps, tau, alpha, 1, tauc, delta)
+    def func(x, deps, tau, alpha, tauc, delta, xaxis: str = 'freq'):
+        return _DerivativeHNWithHF.func(x, deps, tau, alpha, 1, tauc, delta, xaxis=xaxis)
 
 
 class DerivativeCDWithHF:
@@ -298,7 +343,22 @@ class DerivativeCDWithHF:
     type = 'Dielectric Spectroscopy'
     params = [r'\Delta\epsilon', r'\tau', r'\gamma', r'\tau_{c}', r'\delta']
     bounds = [(0, None), (0, None), (0, 1), (0, None), (0, 1)]
+    choices = [
+        ('x axis', 'xaxis', {'Frequency': 'freq', 'Omega': 'omega'})
+    ]
 
     @staticmethod
-    def func(x, deps, tau, gamma, tauc, delta):
-        return _DerivativeHNWithHF.func(x, deps, tau, 1, gamma, tauc, delta)
+    def func(x, deps, tau, gamma, tauc, delta, xaxis: str = 'freq'):
+        return _DerivativeHNWithHF.func(x, deps, tau, 1, gamma, tauc, delta, xaxis=xaxis)
+
+
+def _convert_x_to_omega(x, xaxis: str = 'freq'):
+    if xaxis not in ['freq', 'omega']:
+        raise ValueError(f'Argument `xaxis` is `freq` or `omega`, given is {xaxis!r}')
+
+    if xaxis == 'freq':
+        _w = 2 * np.pi * x
+    else:
+        _w = x
+
+    return _w