diff --git a/psrdada_cpp/effelsberg/edd/CMakeLists.txt b/psrdada_cpp/effelsberg/edd/CMakeLists.txt
index a7f75c58fd786f0c5d81ecde28039a7d054ac099..96ac85460dd3321e50da2b2c55258d4046a68f44 100644
--- a/psrdada_cpp/effelsberg/edd/CMakeLists.txt
+++ b/psrdada_cpp/effelsberg/edd/CMakeLists.txt
@@ -16,9 +16,7 @@ set(psrdada_cpp_effelsberg_edd_src
     src/EDDRoach.cpp
     src/EDDRoach_merge.cpp
     src/ScaledTransposeTFtoTFT.cu
-    src/SKRfiReplacement.cpp
     src/SKRfiReplacementCuda.cu
-    src/SpectralKurtosis.cpp
     src/SpectralKurtosisCuda.cu
     src/Tools.cu
     src/Unpacker.cu
diff --git a/psrdada_cpp/effelsberg/edd/SKRfiReplacementCuda.cuh b/psrdada_cpp/effelsberg/edd/SKRfiReplacementCuda.cuh
index 2ab940dddefb81a5a047204b2f2ce3ef44bf9466..9683181ca421e663b7bda2f6a6bde27d22332c0b 100644
--- a/psrdada_cpp/effelsberg/edd/SKRfiReplacementCuda.cuh
+++ b/psrdada_cpp/effelsberg/edd/SKRfiReplacementCuda.cuh
@@ -1,3 +1,6 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SKRFIREPLACEMENTCUDA_CUH
+#define PSRDADA_CPP_EFFELSBERG_EDD_SKRFIREPLACEMENTCUDA_CUH
+
 #include "psrdada_cpp/common.hpp"
 #include <thrust/host_vector.h>
 #include <thrust/device_vector.h>
@@ -15,8 +18,6 @@ namespace psrdada_cpp {
 namespace effelsberg {
 namespace edd {
 
-#define DEFAULT_NUM_CLEAN_WINDOWS 1 //number of clean windows used for computing DataStatistics 
-
 class SKRfiReplacementCuda{
 public:
     /**
@@ -32,11 +33,13 @@ public:
     /**
      * @brief    Replaces data in rfi_windows with replacement data (generated using statistics of data from clean_windows).
      *
-     * @param(in)          rfi_status    rfi_status of input data
-     * @param(in & out)    data          Data on which RFI has to be replaced. Returns the same but with RFI replaced.   
+     * @param(in)          rfi_status      rfi_status of input data
+     * @param(in & out)    data            Data on which RFI has to be replaced. Returns the same but with RFI replaced.   
+     * @param(in)          clean_windows   number of clean windows used for computing data statistics.
      */
-    void replace_rfi_data(const thrust::device_vector<int> &rfi_status,
-                          thrust::device_vector<thrust::complex<float>> &data);
+    void replace_rfi_data(const thrust::device_vector<int> &rfi_status, 
+                          thrust::device_vector<thrust::complex<float>> &data,
+                          std::size_t clean_windows = 5);
 
 private:
     /**
@@ -70,6 +73,7 @@ private:
     thrust::device_vector<int> _rfi_status;
     std::size_t _window_size;
     std::size_t _nwindows, _nrfi_windows, _nclean_windows;
+    std::size_t _nclean_windows_stat; //number of clean windows used for computing DataStatistics 
     thrust::device_vector<int> _rfi_window_indices;
     thrust::device_vector<int> _clean_window_indices;
     thrust::device_vector<thrust::complex <float>> _clean_data;
@@ -78,3 +82,5 @@ private:
 } //edd
 } //effelsberg
 } //psrdada_cpp
+
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SKRFIREPLACEMENTCUDA_CUH
diff --git a/psrdada_cpp/effelsberg/edd/SpectralKurtosisCuda.cuh b/psrdada_cpp/effelsberg/edd/SpectralKurtosisCuda.cuh
index df5ae9624cddfbe9d278dd1bf8db8cfbfcc71743..a3dc3ca016af056422db0036b0ddd9f13f9c55ac 100644
--- a/psrdada_cpp/effelsberg/edd/SpectralKurtosisCuda.cuh
+++ b/psrdada_cpp/effelsberg/edd/SpectralKurtosisCuda.cuh
@@ -1,3 +1,6 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISCUDA_CUH
+#define PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISCUDA_CUH
+
 #include "psrdada_cpp/common.hpp"
 #include <thrust/device_vector.h>
 #include <thrust/host_vector.h>
@@ -14,8 +17,6 @@ namespace psrdada_cpp {
 namespace effelsberg {
 namespace edd {
 
-#define DEFAULT_SK_MIN 0.9
-#define DEFAULT_SK_MAX 1.1
 #define BLOCK_DIM 1024
 
 struct RFIStatistics{
@@ -33,27 +34,28 @@ public:
      * @param(in)  sk_min        minimum value of spectral kurtosis.
      * @param(in)  sk_max        maximum value of spectral kurtosis.
      */
-    SpectralKurtosisCuda(std::size_t nchannels, std::size_t window_size, float sk_min = DEFAULT_SK_MIN,
-                         float sk_max = DEFAULT_SK_MAX);
+    SpectralKurtosisCuda(std::size_t nchannels, std::size_t window_size, float sk_min = 0.8,
+                         float sk_max = 1.2);
     ~SpectralKurtosisCuda();
 
     /**
-     * @brief      computes spectral kurtosis for the given data and returns its rfi statistics.
+     * @brief      computes spectral kurtosis for the given data using thrust calls and returns its rfi statistics.
+     *             This is used to compare and verify  the results generated by compute_sk - an optimized kernel implementation.
      *
      * @param(in)  data          input data
      * @param(out) stats         RFI statistics
      *
      */
-    void compute_sk(thrust::device_vector<thrust::complex<float>> const& data, RFIStatistics &stats);
+    void compute_sk_thrust(thrust::device_vector<thrust::complex<float>> const& data, RFIStatistics &stats);
 
     /**
-     * @brief      computes spectral kurtosis (using optimized kernel fucntion) for the given data and returns its rfi statistics.
+     * @brief      computes spectral kurtosis (using optimized kernel function) for the given data and returns its rfi statistics.
      *
      * @param(in)  data          input data
      * @param(out) stats         RFI statistics
      *
      */
-    void compute_sk_k(thrust::device_vector<thrust::complex<float>> &data, RFIStatistics &stats);
+    void compute_sk(thrust::device_vector<thrust::complex<float>> &data, RFIStatistics &stats);
 
 private:
     /**
@@ -67,7 +69,10 @@ private:
     std::size_t _sample_size; //size of input data
     float _sk_min, _sk_max;
     thrust::device_vector<float> _d_s1, _d_s2;
+    cudaStream_t _stream;
 };
 } //edd
 } //effelsberg
 } //psrdada_cpp
+
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISCUDA_CUH
diff --git a/psrdada_cpp/effelsberg/edd/src/SKRfiReplacementCuda.cu b/psrdada_cpp/effelsberg/edd/src/SKRfiReplacementCuda.cu
index a3edc67bdd95e61ffd4e39e930b3f5300a002d4d..9d9d914f9656a3ddde9b3bd68cb36e52d9f2b88d 100644
--- a/psrdada_cpp/effelsberg/edd/src/SKRfiReplacementCuda.cu
+++ b/psrdada_cpp/effelsberg/edd/src/SKRfiReplacementCuda.cu
@@ -66,7 +66,13 @@ void SKRfiReplacementCuda::get_clean_window_indices()
 {
     nvtxRangePushA("get_clean_window_indices");
     _nclean_windows = thrust::count(_rfi_status.begin(), _rfi_status.end(), 0);
-    _clean_window_indices.resize(DEFAULT_NUM_CLEAN_WINDOWS);
+    if (_nclean_windows > _nclean_windows_stat) _clean_window_indices.resize(_nclean_windows_stat);
+    else{
+        BOOST_LOG_TRIVIAL(info) << "Total number of clean windows is less than number of clean windows (_nclean_windows_stat)"
+                                << "chosen to compute data statistics.\n"
+                                <<" Therefore choosing total clean windows for statistics computation.\n";
+        _clean_window_indices.resize(_nclean_windows);
+    }
     thrust::copy_if(thrust::make_counting_iterator<int>(0),
                     thrust::make_counting_iterator<int>(_nwindows),
                     _rfi_status.begin(),
@@ -79,8 +85,8 @@ void SKRfiReplacementCuda::get_clean_data_statistics(const thrust::device_vector
 {
     nvtxRangePushA("get_clean_data_statistics");
     _window_size = data.size() / _nwindows;
-    _clean_data.resize(DEFAULT_NUM_CLEAN_WINDOWS * _window_size);
-    for(std::size_t ii = 0; ii < DEFAULT_NUM_CLEAN_WINDOWS; ii++){
+    _clean_data.resize(_nclean_windows_stat * _window_size);
+    for(std::size_t ii = 0; ii < _nclean_windows_stat; ii++){
         std::size_t window_index = _clean_window_indices[ii];
         std::size_t ibegin = window_index * _window_size;
         std::size_t iend = ibegin + _window_size - 1;
@@ -97,6 +103,8 @@ void SKRfiReplacementCuda::compute_clean_data_statistics()
 {
     nvtxRangePushA("compute_clean_data_statistics");
     std::size_t length = _clean_data.size();
+    //The distribution of both real and imag have same mean and standard deviation. 
+    //Therefore computing _ref_mean, _ref_sd for real distribution only.  
     _ref_mean = (thrust::reduce(_clean_data.begin(), _clean_data.end(), thrust::complex<float> (0.0f, 0.0f))). real() / length;
     _ref_sd = std::sqrt(thrust::transform_reduce(_clean_data.begin(), _clean_data.end(), mean_subtraction_square(_ref_mean),
                         0.0f, thrust::plus<float> ()) / length);
@@ -106,15 +114,17 @@ void SKRfiReplacementCuda::compute_clean_data_statistics()
 }
 
 void SKRfiReplacementCuda::replace_rfi_data(const thrust::device_vector<int> &rfi_status,
-                                            thrust::device_vector<thrust::complex<float>> &data)
+                                            thrust::device_vector<thrust::complex<float>> &data, 
+                                            std::size_t clean_windows)
 {
     nvtxRangePushA("replace_rfi_data");
     _rfi_status = rfi_status;
     thrust::device_vector<thrust::complex<float>> replacement_data;
     //initialize data members of the class
+    _nclean_windows_stat = clean_windows; //no. of clean windows used for computing statistics
     init();
     //RFI present and not in all windows
-    if((_nrfi_windows > 0) && (_nrfi_windows < _nwindows)){
+    if(_nclean_windows < _nwindows){
         get_clean_data_statistics(data);
         //Replacing RFI
         thrust::counting_iterator<unsigned int> sequence_index_begin(0);
diff --git a/psrdada_cpp/effelsberg/edd/src/SpectralKurtosisCuda.cu b/psrdada_cpp/effelsberg/edd/src/SpectralKurtosisCuda.cu
index a8718e861e7d7a80590a4c1d96cc11c5490fa771..ee07975b74788ba4f95666fa78d7893f91e22958 100755
--- a/psrdada_cpp/effelsberg/edd/src/SpectralKurtosisCuda.cu
+++ b/psrdada_cpp/effelsberg/edd/src/SpectralKurtosisCuda.cu
@@ -97,11 +97,13 @@ SpectralKurtosisCuda::SpectralKurtosisCuda(std::size_t nchannels, std::size_t wi
       _sk_max(sk_max)
 {
     BOOST_LOG_TRIVIAL(debug) << "Creating new SpectralKurtosisCuda instance... \n";
+    cudaStreamCreate(&_stream);
 }
 
 SpectralKurtosisCuda::~SpectralKurtosisCuda()
 {
     BOOST_LOG_TRIVIAL(debug) << "Destroying SpectralKurtosisCuda instance... \n";
+    cudaStreamDestroy(_stream);
 }
 
 void SpectralKurtosisCuda::init()
@@ -112,17 +114,17 @@ void SpectralKurtosisCuda::init()
         throw std::runtime_error("Data(sample) size is not a multiple of window_size. Give different window size. \n");
     }
     _nwindows = _sample_size /_window_size;
-    _d_s1.resize(_nwindows);
-    _d_s2.resize(_nwindows);
 }
 
-void SpectralKurtosisCuda::compute_sk(const thrust::device_vector<thrust::complex<float>> &data, RFIStatistics &stats){
+void SpectralKurtosisCuda::compute_sk_thrust(const thrust::device_vector<thrust::complex<float>> &data, RFIStatistics &stats){
     nvtxRangePushA("compute_sk");
     _sample_size = data.size();
     BOOST_LOG_TRIVIAL(debug) << "Computing SK (thrust version) for sample_size " << _sample_size
                              << " and window_size " << _window_size <<".\n";
     //initializing class variables
     init();
+    _d_s1.resize(_nwindows);
+    _d_s2.resize(_nwindows);
     //computing _d_s1 for all windows
     nvtxRangePushA("compute_sk_reduce_by_key_call");
     thrust::reduce_by_key(thrust::device, 
@@ -149,20 +151,18 @@ void SpectralKurtosisCuda::compute_sk(const thrust::device_vector<thrust::comple
     nvtxRangePop();
 }
 
-void SpectralKurtosisCuda::compute_sk_k(thrust::device_vector<thrust::complex<float>> &data, RFIStatistics &stats){
+void SpectralKurtosisCuda::compute_sk(thrust::device_vector<thrust::complex<float>> &data, RFIStatistics &stats){
     nvtxRangePushA("compute_sk_kernel");
     _sample_size = data.size();
     BOOST_LOG_TRIVIAL(debug) << "Computing SK (kernel version) for sample_size " << _sample_size
                              << " and window_size " << _window_size <<".\n";
-
-    _nwindows = _sample_size / _window_size;
+    //initializing variables
+    init();
     stats.rfi_status.resize(_nwindows);
-    thrust::complex<float> *k_data = thrust::raw_pointer_cast(data.data());
     int *k_rfi_status = thrust::raw_pointer_cast(stats.rfi_status.data());
-    int blockSize = BLOCK_DIM;
-    int gridSize = _nwindows; 
+    thrust::complex<float> *k_data = thrust::raw_pointer_cast(data.data());
     nvtxRangePushA("compute_sk_kernel_call");
-    compute_sk_kernel<<<gridSize, blockSize>>> (k_data, _sample_size, _window_size, _sk_max, _sk_min, k_rfi_status);
+    compute_sk_kernel<<<_nwindows, BLOCK_DIM, 0, _stream>>> (k_data, _sample_size, _window_size, _sk_max, _sk_min, k_rfi_status);
     cudaDeviceSynchronize();
     nvtxRangePop();
     nvtxRangePushA("compute_sk_kernel_rfi_fraction");
diff --git a/psrdada_cpp/effelsberg/edd/test/CMakeLists.txt b/psrdada_cpp/effelsberg/edd/test/CMakeLists.txt
index 8502287d909c6143f87e5ee0bc22ef272425165d..d665eee7dd64fb798f86012dee0f2354c7e1040d 100644
--- a/psrdada_cpp/effelsberg/edd/test/CMakeLists.txt
+++ b/psrdada_cpp/effelsberg/edd/test/CMakeLists.txt
@@ -14,6 +14,8 @@ set(
     src/VLBITest.cu
     src/EDDPolnMergeTester.cpp
     src/SKTestVector.cpp
+    src/SpectralKurtosis.cpp
+    src/SKRfiReplacement.cpp
     src/SpectralKurtosisTester.cpp
     src/SpectralKurtosisCudaTester.cu
 )
diff --git a/psrdada_cpp/effelsberg/edd/test/ReadData.hpp b/psrdada_cpp/effelsberg/edd/test/ReadData.hpp
deleted file mode 100644
index f5df08f2575d31c24d85ef8998f021691494d9af..0000000000000000000000000000000000000000
--- a/psrdada_cpp/effelsberg/edd/test/ReadData.hpp
+++ /dev/null
@@ -1,27 +0,0 @@
-#include "psrdada_cpp/common.hpp"
-
-#include <fstream>
-#include <complex>
-#include <vector>
-
-namespace psrdada_cpp {
-namespace effelsberg {
-namespace edd {
-
-# define DADA_HDR_SIZE 4096
-
-class ReadData
-{
-public:
-    std::vector<std::complex<float>> _pol0, _pol1;
-    std::size_t _sample_size;
-    //int _sample_size;
-    ReadData(std::string filename);
-    void read_file();
-private:
-    std::string _filename;
-};
-
-} //edd
-} //effelsberg
-} //psrdada_cpp
diff --git a/psrdada_cpp/effelsberg/edd/SKRfiReplacement.hpp b/psrdada_cpp/effelsberg/edd/test/SKRfiReplacement.hpp
similarity index 94%
rename from psrdada_cpp/effelsberg/edd/SKRfiReplacement.hpp
rename to psrdada_cpp/effelsberg/edd/test/SKRfiReplacement.hpp
index fc85234c0181f2a2972f1ae458b68e8b708b573d..7abf6c289fc7161ee8bea7d327f172199c7283f2 100644
--- a/psrdada_cpp/effelsberg/edd/SKRfiReplacement.hpp
+++ b/psrdada_cpp/effelsberg/edd/test/SKRfiReplacement.hpp
@@ -1,3 +1,6 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SKRFIREPLACEMENT_HPP
+#define PSRDADA_CPP_EFFELSBERG_EDD_SKRFIREPLACEMENT_HPP
+
 #include "psrdada_cpp/common.hpp"
 #include <vector>
 #include <complex>
@@ -89,3 +92,5 @@ private:
 } //edd
 } //effelsberg
 } //psrdada_cpp
+
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SKRFIREPLACEMENT_HPP
diff --git a/psrdada_cpp/effelsberg/edd/test/SKTestVector.hpp b/psrdada_cpp/effelsberg/edd/test/SKTestVector.hpp
index a103a7017c3573f8ed5f5ff9b0652a44cd11191e..b115c7220db8c0c8672a758041f969d9693f0008 100644
--- a/psrdada_cpp/effelsberg/edd/test/SKTestVector.hpp
+++ b/psrdada_cpp/effelsberg/edd/test/SKTestVector.hpp
@@ -1,3 +1,6 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SKTESTVECTOR_HPP
+#define PSRDADA_CPP_EFFELSBERG_EDD_SKTESTVECTOR_HPP
+
 #include "psrdada_cpp/common.hpp"
 
 #include <complex>
@@ -67,3 +70,4 @@ private:
 } //effelsberg
 } //psrdada_cpp
 
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SKTESTVECTOR_HPP
diff --git a/psrdada_cpp/effelsberg/edd/SpectralKurtosis.hpp b/psrdada_cpp/effelsberg/edd/test/SpectralKurtosis.hpp
similarity index 86%
rename from psrdada_cpp/effelsberg/edd/SpectralKurtosis.hpp
rename to psrdada_cpp/effelsberg/edd/test/SpectralKurtosis.hpp
index 3f09ac810c58808216b3a8c0e4e9631d71a9b0ba..0f576940fad4537a92bde833800536844b2229db 100644
--- a/psrdada_cpp/effelsberg/edd/SpectralKurtosis.hpp
+++ b/psrdada_cpp/effelsberg/edd/test/SpectralKurtosis.hpp
@@ -1,3 +1,6 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSIS_HPP
+#define PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSIS_HPP
+
 #include "psrdada_cpp/common.hpp"
 #include <complex>
 #include <vector>
@@ -7,9 +10,6 @@ namespace psrdada_cpp {
 namespace effelsberg {
 namespace edd {
 
-#define DEFAULT_SK_MIN 0.9
-#define DEFAULT_SK_MAX 1.1
-
 struct RFIStatistics{
     std::vector<int>  rfi_status;
     float rfi_fraction;
@@ -25,8 +25,8 @@ public:
      *            sk_min        minimum value of spectral kurtosis.
      *            sk_max        maximum value of spectral kurtosis.
      */
-    SpectralKurtosis(std::size_t nchannels, std::size_t window_size, float sk_min = DEFAULT_SK_MIN,
-		     float sk_max = DEFAULT_SK_MAX);
+    SpectralKurtosis(std::size_t nchannels, std::size_t window_size, float sk_min = 0.9,
+		     float sk_max = 1.1);
     ~SpectralKurtosis();
 
     /**
@@ -56,3 +56,4 @@ private:
 } //effelsberg
 } //psrdada_cpp
 
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSIS_HPP
diff --git a/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisCudaTester.cuh b/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisCudaTester.cuh
index 529c500889ac50f602f7679483fc77ae500d90a8..9544a749c75a62ffe0e4475e07546bb27861b0f7 100644
--- a/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisCudaTester.cuh
+++ b/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisCudaTester.cuh
@@ -1,3 +1,6 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISCUDATESTER_CUH
+#define PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISCUDATESTER_CUH
+
 #include "psrdada_cpp/common.hpp"
 #include "psrdada_cpp/effelsberg/edd/test/SKTestVector.hpp"
 #include "psrdada_cpp/effelsberg/edd/SpectralKurtosisCuda.cuh"
@@ -40,6 +43,4 @@ protected:
 } //effelsberg
 } //psrdada_cpp
 
-
-
-
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISCUDATESTER_CUH
diff --git a/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisTester.hpp b/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisTester.hpp
index 39f303abbafabf036260bdd719b94f0b41718f5c..bac613a9cc386bae6ea99f38fffeb2ecc5c25de1 100644
--- a/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisTester.hpp
+++ b/psrdada_cpp/effelsberg/edd/test/SpectralKurtosisTester.hpp
@@ -1,6 +1,9 @@
+#ifndef PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISTESTER_HPP
+#define PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISTESTER_HPP
+
 #include "psrdada_cpp/common.hpp"
-#include "psrdada_cpp/effelsberg/edd/SKRfiReplacement.hpp"
-#include "psrdada_cpp/effelsberg/edd/SpectralKurtosis.hpp"
+#include "psrdada_cpp/effelsberg/edd/test/SKRfiReplacement.hpp"
+#include "psrdada_cpp/effelsberg/edd/test/SpectralKurtosis.hpp"
 #include "psrdada_cpp/effelsberg/edd/test/SKTestVector.hpp"
 #include <gtest/gtest.h>
 
@@ -38,6 +41,4 @@ protected:
 } //effelsberg
 } //psrdada_cpp
 
-
-
-
+#endif //PSRDADA_CPP_EFFELSBERG_EDD_SPECTRALKURTOSISTESTER_HPP
diff --git a/psrdada_cpp/effelsberg/edd/src/SKRfiReplacement.cpp b/psrdada_cpp/effelsberg/edd/test/src/SKRfiReplacement.cpp
similarity index 98%
rename from psrdada_cpp/effelsberg/edd/src/SKRfiReplacement.cpp
rename to psrdada_cpp/effelsberg/edd/test/src/SKRfiReplacement.cpp
index dd6d6ed00274961d4e5e1a4575ec3cabf86aed31..7623efed2bd99bf55b1f01e7d011befd49f1a7c9 100644
--- a/psrdada_cpp/effelsberg/edd/src/SKRfiReplacement.cpp
+++ b/psrdada_cpp/effelsberg/edd/test/src/SKRfiReplacement.cpp
@@ -1,4 +1,4 @@
-#include "psrdada_cpp/effelsberg/edd/SKRfiReplacement.hpp"
+#include "psrdada_cpp/effelsberg/edd/test/SKRfiReplacement.hpp"
 
 namespace psrdada_cpp {
 namespace effelsberg {
diff --git a/psrdada_cpp/effelsberg/edd/src/SpectralKurtosis.cpp b/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosis.cpp
similarity index 97%
rename from psrdada_cpp/effelsberg/edd/src/SpectralKurtosis.cpp
rename to psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosis.cpp
index 8fa9aba2b6c96c4de9ff1db605437e0c02fba8df..e49480cc6b04b40148171c70a384d41f56b4a568 100644
--- a/psrdada_cpp/effelsberg/edd/src/SpectralKurtosis.cpp
+++ b/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosis.cpp
@@ -1,4 +1,4 @@
-#include "psrdada_cpp/effelsberg/edd/SpectralKurtosis.hpp"
+#include "psrdada_cpp/effelsberg/edd/test/SpectralKurtosis.hpp"
 
 namespace psrdada_cpp {
 namespace effelsberg {
diff --git a/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisCudaTester.cu b/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisCudaTester.cu
index ca4e561f15764a7198b0e186643f7490d734af72..d6129e99e0a8339c1a0fcdcd946cb9c2555bcef9 100644
--- a/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisCudaTester.cu
+++ b/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisCudaTester.cu
@@ -49,20 +49,32 @@ TEST_F(SpectralKurtosisCudaTester, sk_window_size_check)
 {
     std::vector<int> rfi_window_indices{};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(4000, 150, 0, 0, 0, rfi_window_indices, samples);
+    std::size_t sample_size = 4000;
+    std::size_t window_size = 150;
+    bool with_rfi = 0;
+    float rfi_freq = 0;
+    float rfi_amp = 0;
+    test_vector_generation(sample_size, window_size, with_rfi, rfi_freq, rfi_amp, rfi_window_indices, samples);
 
     RFIStatistics stat;
-    EXPECT_THROW(sk_computation(1, 150, samples, stat), std::runtime_error);
+    std::size_t nch = 1;
+    EXPECT_THROW(sk_computation(nch, window_size, samples, stat), std::runtime_error);
 }
 
 TEST_F(SpectralKurtosisCudaTester, sk_withoutRFI)
 {
     std::vector<int> rfi_window_indices{};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(40000, 400, 0, 0, 0, rfi_window_indices, samples);
+    std::size_t sample_size = 40000;
+    std::size_t window_size = 400;
+    bool with_rfi = 0;
+    float rfi_freq = 0;
+    float rfi_amp = 0;
+    test_vector_generation(sample_size, window_size, with_rfi, rfi_freq, rfi_amp, rfi_window_indices, samples);
 
     RFIStatistics stat;
-    sk_computation(1, 400, samples, stat);
+    std::size_t nch = 1;
+    sk_computation(nch, window_size, samples, stat);
     float expected_rfi_fraction = 0;
     EXPECT_EQ(expected_rfi_fraction, stat.rfi_fraction);
 }
@@ -71,34 +83,41 @@ TEST_F(SpectralKurtosisCudaTester, sk_withRFI)
 {
     std::vector<int> rfi_window_indices{3, 4, 6, 7, 8, 20, 30, 40, 45, 75};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(40000, 400, 1, 30, 10, rfi_window_indices, samples);
+    std::size_t sample_size = 40000;
+    std::size_t window_size = 400;
+    bool with_rfi = 1;
+    float rfi_freq = 30;
+    float rfi_amp = 10;
+    test_vector_generation(sample_size, window_size, with_rfi, rfi_freq, rfi_amp, rfi_window_indices, samples);
 
     RFIStatistics stat;
-    sk_computation(1, 400, samples, stat);
-    float expected_rfi_fraction = (rfi_window_indices.size()/float(40000/400));
+    std::size_t nch = 1;
+    sk_computation(nch, window_size, samples, stat);
+    float expected_rfi_fraction = (rfi_window_indices.size()/float(sample_size/window_size));
     EXPECT_EQ(expected_rfi_fraction, stat.rfi_fraction);
 }
 
 TEST_F(SpectralKurtosisCudaTester, sk_RFIreplacement)
 {
     std::size_t sample_size = 128* 1024 * 1024;
-    //std::size_t sample_size = 160000000;
     std::size_t window_size = 1024 * 2;
     //Test vector generation
     std::vector<int> rfi_window_indices{3, 4, 6, 7, 8, 20, 30, 40, 45, 75};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(sample_size, window_size, 1, 30, 10, rfi_window_indices, samples);
+    bool with_rfi = 1;
+    float rfi_freq = 30;
+    float rfi_amp = 10;
+    test_vector_generation(sample_size, window_size, with_rfi, rfi_freq, rfi_amp, rfi_window_indices, samples);
 
     //SK computation
     thrust::host_vector<thrust::complex<float>> h_samples(samples);
     thrust::device_vector<thrust::complex<float>> d_samples(h_samples);
     float sk_min = 0.8;
     float sk_max = 1.2;    
-    SpectralKurtosisCuda sk(1, window_size, sk_min, sk_max);
+    std::size_t nch = 1;
+    SpectralKurtosisCuda sk(nch, window_size, sk_min, sk_max);
     RFIStatistics stat;
     sk.compute_sk(d_samples, stat);
-    //float expected_rfi_fraction = (rfi_window_indices.size()/float(sample_size/window_size));
-    //EXPECT_EQ(expected_rfi_fraction, stat.rfi_fraction);
 
     //RFI replacement
     BOOST_LOG_TRIVIAL(info) <<"RFI replacement..\n";
@@ -118,20 +137,22 @@ TEST_F(SpectralKurtosisCudaTester, sk_kernel)
     std::size_t window_size = 2000;
     std::size_t nwindows = sample_size / window_size;
     //Test vector generation
-    std::vector<int> rfi_window_indices{3, 4, 6, 7, 8, 20, 30, 40, 45, 75};
+    std::vector<int> rfi_window_indices{1, 4, 6, 7, 8, 20, 30, 40, 45, 75};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(sample_size, window_size, 1, 30, 10, rfi_window_indices, samples);
+    bool with_rfi = 1;
+    float rfi_freq = 30;
+    float rfi_amp = 10;
+    test_vector_generation(sample_size, window_size, with_rfi, rfi_freq, rfi_amp, rfi_window_indices, samples);
 
     //SK computation
     thrust::host_vector<thrust::complex<float>> h_samples(samples);
     thrust::device_vector<thrust::complex<float>> d_samples(h_samples);
 
-    float sk_min = 0.8;
-    float sk_max = 1.2;    
-    SpectralKurtosisCuda sk(1, window_size, sk_min, sk_max);
+    std::size_t nch = 1;    
+    SpectralKurtosisCuda sk(nch, window_size);
     RFIStatistics stat, stat_k;
-    sk.compute_sk(d_samples, stat);
-    sk.compute_sk_k(d_samples, stat_k);
+    sk.compute_sk_thrust(d_samples, stat);
+    sk.compute_sk(d_samples, stat_k);
     for (int ii = 0; ii < nwindows; ii++){
         EXPECT_EQ(stat.rfi_status[ii], stat_k.rfi_status[ii]);
     }
@@ -140,11 +161,12 @@ TEST_F(SpectralKurtosisCudaTester, sk_kernel)
     //RFI replacement
     BOOST_LOG_TRIVIAL(info) <<"RFI replacement..\n";
     SKRfiReplacementCuda rr;
-    rr.replace_rfi_data(stat_k.rfi_status, d_samples);
+    std::size_t clean_windows = 100; //no. of clean windows used for computing data statistics
+    rr.replace_rfi_data(stat_k.rfi_status, d_samples, clean_windows);
 
     //SK computation after RFI replacement
     BOOST_LOG_TRIVIAL(info) <<"computing SK after replacing the RFI data..\n";
-    sk.compute_sk_k(d_samples, stat_k);
+    sk.compute_sk(d_samples, stat_k);
     float expected_val_after_rfi_replacement = 0;
     EXPECT_EQ(expected_val_after_rfi_replacement, stat_k.rfi_fraction);
 }
diff --git a/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisTester.cpp b/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisTester.cpp
index f68a5c36fb582b126283fb31963d7233fdad787d..9e45a87b456dc7c915a23309bc605c57a412c0b8 100644
--- a/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisTester.cpp
+++ b/psrdada_cpp/effelsberg/edd/test/src/SpectralKurtosisTester.cpp
@@ -47,21 +47,27 @@ TEST_F(SpectralKurtosisTester, sk_window_size_check)
 {
     std::vector<int> rfi_window_indices{};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(4000, 150, 0, 0, 0, rfi_window_indices, samples);
+    std::size_t sample_size = 4000;
+    std::size_t window_size = 150;
+    test_vector_generation(sample_size, window_size, 0, 0, 0, rfi_window_indices, samples);
 
     RFIStatistics stat;
-    EXPECT_THROW(sk_computation(1, 150, samples, stat), std::runtime_error);
+    std::size_t nch = 1;
+    EXPECT_THROW(sk_computation(nch, window_size, samples, stat), std::runtime_error);
 }
 
 TEST_F(SpectralKurtosisTester, sk_withoutRFI)
 {
     std::vector<int> rfi_window_indices{};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(40000, 400, 0, 0, 0, rfi_window_indices, samples);
+    std::size_t sample_size = 4000;
+    std::size_t window_size = 400;
+    test_vector_generation(sample_size, window_size, 0, 0, 0, rfi_window_indices, samples);
 
     RFIStatistics stat;
-    sk_computation(1, 400, samples, stat);
-    float expected_rfi_fraction = 0.01;
+    std::size_t nch = 1;
+    sk_computation(nch, window_size, samples, stat);
+    float expected_rfi_fraction = 0;
     EXPECT_EQ(expected_rfi_fraction, stat.rfi_fraction);
 }
 
@@ -69,11 +75,14 @@ TEST_F(SpectralKurtosisTester, sk_withRFI)
 {
     std::vector<int> rfi_window_indices{3, 4, 6, 7, 8, 20, 30, 40};
     std::vector<std::complex<float>> samples;
-    test_vector_generation(40000, 400, 1, 30, 10, rfi_window_indices, samples);
+    std::size_t sample_size = 40000;
+    std::size_t window_size = 400;
+    test_vector_generation(sample_size, window_size, 1, 10, 30, rfi_window_indices, samples);
 
     RFIStatistics stat;
-    sk_computation(1, 400, samples, stat);
-    float expected_rfi_fraction = (rfi_window_indices.size()/float(40000/400)) + 0.01;
+    std::size_t nch = 1;
+    sk_computation(nch, window_size, samples, stat);
+    float expected_rfi_fraction = (rfi_window_indices.size()/float(sample_size/window_size)) + 0.01;
     EXPECT_EQ(expected_rfi_fraction, stat.rfi_fraction); //To check: fails inspite of actual and expected values being same.
 }
 
@@ -99,7 +108,7 @@ TEST_F(SpectralKurtosisTester, sk_replacement)
     RFIStatistics stat;
     SpectralKurtosis sk(nch, window_size, sk_min, sk_max);
     sk.compute_sk(samples, stat); //computing SK
-    float expected_rfi_fraction = (rfi_window_indices.size()/float(40000/400)) + 0.01;
+    float expected_rfi_fraction = (rfi_window_indices.size()/float(sample_size/window_size)) + 0.01;
     EXPECT_EQ(expected_rfi_fraction, stat.rfi_fraction);
 
     //RFI replacement