diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index 6108b8c..c8acd75 100644
--- a/src/cosmopower.cpp
+++ b/src/cosmopower.cpp
@@ -269,6 +269,19 @@ double CosmoPower::matterTransferFunctionHu(double k)
   return T_k;
 }
 
+double CosmoPower::noWiggleTk(double k)
+{
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+
+  return T0;
+}
+
 double CosmoPower::powerHuBaryons(double k)
 {
   double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
@@ -473,6 +486,9 @@ void CosmoPower::setFunction(CosmoFunction f)
     case POWER_TEST:
       eval = &CosmoPower::powerTest;
       break;
+    case NOWIGGLE_TK:
+      eval = &CosmoPower::noWiggleTk;
+      break;
     default:
       abort();
     }
diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp
index 7b02f5f..5ed9b61 100644
--- a/src/cosmopower.hpp
+++ b/src/cosmopower.hpp
@@ -90,7 +90,8 @@ namespace CosmoTool {
       POWER_BDM,
       POWER_TEST,
       HU_WIGGLES_ORIGINAL,
-      SAMPLE_WIGGLES
+      SAMPLE_WIGGLES,
+      NOWIGGLE_TK
     };
 
     CosmoPower();
@@ -124,6 +125,7 @@ namespace CosmoTool {
     double powerTest(double k);
     double powerHuWigglesOriginal(double k);
     double sample_BAO(double k);
+    double noWiggleTk(double k);
   };
 
 };
diff --git a/src/powerSpectrum.cpp b/src/powerSpectrum.cpp
index 777654c..7786080 100644
--- a/src/powerSpectrum.cpp
+++ b/src/powerSpectrum.cpp
@@ -58,6 +58,7 @@ using namespace std;
 #define POWER_TEST 8
 
 #define POWER_SPECTRUM POWER_EFSTATHIOU
+#define SAMPLE_WIGGLES 9
 
 namespace Cosmology {
 
@@ -206,6 +207,31 @@ double powG(double y)
  */
 double powerSpectrum(double k, double normPower)
 {
+#if POWER_SPECTRUM == SAMPLE_WIGGLES
+  // BAO wiggle parameterization for reconstruction
+  // Babic et al. 2022, https://arxiv.org/abs/2203.06177
+
+  // No-wiggle transfer function
+
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+
+  // Wiggle parameterization
+
+  double A = 0;
+  double r_s = 10;
+  double k_D = 2 * M_PI / 100;
+
+  double param = 1 + A * sin(k * r_s) * exp(- k / k_D);
+
+  return normPower * pow(k, n) * T0 * T0 * param;
+#endif
+
 #if POWER_SPECTRUM == POWER_EFSTATHIOU
   double a = 6.4/Gamma0;
   double b = 3/Gamma0;