Speech detection and separation, noise cancellation AIMs are submitted.

AnalysisTransform/analysistransform.cpp

@@ -25,7 +25,7 @@ analysistransform::analysistransform(int windowsize, int signallength, int hopsi
     this->nofmics = nofmics;
 
     // Create a hamming window of appropriate length
-    window = (double*)malloc(sizeof(double) * this->windowsize);
+    window = (double*)malloc(sizeof(double) * (this->windowsize));
 
     hamming(this->windowsize, window);
 }
@@ -33,9 +33,9 @@ analysistransform::analysistransform(int windowsize, int signallength, int hopsi
 // Create a hamming window of windowLength samples in buffer
 void analysistransform::hamming(int windowLength, double* buffer) {
 
-    for (int i = 0; i < windowLength; i++) {
+    for (int i = 0; i <= windowLength; i++) {
 
-        buffer[i] = 0.54 - (0.46 * cos(2 * PI * (i / (((double)windowLength - 1.0) * 1.0))));
+        buffer[i] = (0.53836 - (0.46164 * cos(2 * PI * ((double)i / (((double)windowLength) * 1.0)))));
     }
 }
 
@@ -84,7 +84,7 @@ void* analysistransform::FFT(double* signal, double* realpart, double* imagpart)
                 if (readIndex < this->signallength * this->nofmics) {
 
                     // Note the windowing!
-                    data[i] = signal[readIndex] * window[i];
+                    data[i] = signal[readIndex] * window[i+1];
 
                 }
                 else  {
@@ -108,7 +108,7 @@ void* analysistransform::FFT(double* signal, double* realpart, double* imagpart)
             // Matlab's spectrogram routine works.
 
 
-            for (i = 0; i < this->windowsize / 2; i++) {
+            for (i = 0; i <= this->windowsize / 2; i++) {
 
                 *(realpart + outputcounter) = fft_result[i][0];
                 *(imagpart + outputcounter) = fft_result[i][1];
@@ -131,3 +131,98 @@ void* analysistransform::FFT(double* signal, double* realpart, double* imagpart)
     return 0;
 
 }
+
+/*
+Signal in : x[0] to x[signalLength - 1]
+nofmics : the number of the microphones
+realpart : fft_real[0] to fft_real[N/2 -1]
+imagpart : imagpart[0] to imagpart[N/2 -1]
+*/
+void* analysistransform::FFT(double* signal, float* realpart, float* imagpart) {
+
+
+    double* data;
+    fftw_complex* fft_result;
+    fftw_plan       plan_forward;
+    int             i;
+
+    data = (double*)fftw_malloc(sizeof(double) * this->windowsize);
+
+    fft_result = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * this->windowsize );
+
+    plan_forward = fftw_plan_dft_r2c_1d(this->windowsize, data, fft_result, FFTW_ESTIMATE);
+
+    int readIndex;
+
+    // Should we stop reading in chunks?
+    int bStop = 0;
+
+    int channelindex = 0;
+
+    // Process each chunk of the signal
+
+    for (channelindex = 0; channelindex < this->nofmics; channelindex++)
+    {
+        int outputcounter = channelindex;
+        int chunkPosition = 0;
+
+        while ((chunkPosition + this->windowsize) <= this->signallength && !bStop) {
+
+            // Copy the chunk into our buffer
+            for (i = 0; i < this->windowsize; i++) {
+
+
+                readIndex = (chunkPosition + i) * this->nofmics + channelindex;
+
+                if (readIndex < this->signallength * this->nofmics) {
+
+                    // Note the windowing!
+                    data[i] = signal[readIndex]; //   *window[i];
+
+                }
+                else {
+
+                    // we have read beyond the signal, so zero-pad it!
+
+                    data[i] = 0.0;
+
+                    bStop = 1;
+
+                }
+            }
+
+            // Perform the FFT on our chunk
+            fftw_execute(plan_forward);
+
+
+            // Copy the first (windowSize/2 + 1) data points into your spectrogram.
+            // We do this because the FFT output is mirrored about the nyquist
+            // frequency, so the second half of the data is redundant. This is how
+            // Matlab's spectrogram routine works.
+
+
+            for (i = 0; i < this->windowsize / 2; i++) {
+//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
+               *(realpart + outputcounter) = (float)fft_result[i][0] / (float)((this->windowsize));
+              *(imagpart + outputcounter) = (float)fft_result[i][1] / (float)((this->windowsize));
+//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
+              //  *(realpart + outputcounter) = ((float)fft_result[i][0] * window[i]) / (float)((this->windowsize) / 2.0);
+              //  *(imagpart + outputcounter) = ((float)fft_result[i][1] * window[i]) / (float)((this->windowsize) / 2.0);
+                outputcounter += this->nofmics;
+            }
+
+
+            chunkPosition += this->overlap;
+
+        } // Excuse the formatting, the while ends here.
+        bStop = 0;
+
+    }
+    fftw_destroy_plan(plan_forward);
+
+    fftw_free(data);
+    fftw_free(fft_result);
+
+    return 0;
+
+}

AnalysisTransform/analysistransform.h

@@ -15,6 +15,7 @@ public:
     analysistransform();
     analysistransform(int windowsize, int signallength, int hopsize,int nofmics);
     void* FFT(double* signal, double* realpart, double* imagpart);
+    void* FFT(double* signal, float* realpart, float* imagpart);
 private:
     void hamming(int windowLength, double* buffer);
     int signallength;

NoiseCancellation/include/healpixkernel_nc.h

+#pragma once
+
+#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string> 
+#include <vector>
+#include <complex>
+
+#include "pixelinfo_nc.h"
+//#include "MatrixSnf.h"
+
+#include "boost/geometry.hpp"
+#include "boost/math/special_functions/bessel.hpp"
+#include "boost/math/special_functions/bessel_prime.hpp"
+#include "boost/math/special_functions/legendre.hpp"
+#include "boost/math/special_functions/hankel.hpp"
+#include "boost/math/special_functions/spherical_harmonic.hpp"
+#include "boost/variant/get.hpp"
+
+#include <armadillo>
+
+#define Order0_HP 12
+#define Order1_HP 48
+#define Order2_HP 192
+#define Order3_HP 762
+#define PI 3.14159
+#define Kappa 4.0
+
+#define HPMaxvalue 19.2
+
+
+using namespace std;
+using namespace boost::math;
+using namespace boost::geometry;
+using namespace arma;
+
+class healpixkernel_nc
+{
+public:
+	healpixkernel_nc(int order);
+	void gethealpixcoordinates(int id, float& theta, float& phi);
+	int* gethealpixneighbours(int id);
+	int gethealpixorder();
+	int gethealpixnumberofpixels();
+	float* calculatelegendrekernels(int max_sphorder);
+	float* getpinvlegendrekernels(int hpnumber);
+	float* getpinvkernels(int hpnumber);
+	void* getmultiplicationresult();
+	double* create_vonmises(int hpnumber_center);
+	double* create_vonmisesfull();
+	void* calculateorthogonalLK();
+	void* calculateSRFMultiplier(int max_sphorder, int nofmicrophones,int freqrange, vector<float> micanglestheta, vector<float> micanglesphi,float* realpart, float* imagpart);
+	void ConvertSphericalToCartesian(float& distance, float& theta, float& phi, float& x, float& y, float& z);
+	void ConvertCartesionToSpherical(float& distance, float& theta, float& phi, float& x, float& y, float& z);
+	float FindAngleBtw(float& x1, float& y1, float& z1, float& mag1, float& x2, float& y2, float& z2, float& mag2);
+	float def_vonmises(float kappa, float costheta);
+	float* calculatepw_magnitude(int Micnumber, float* micanglestheta, float* micanglesphi);
+	void gethealpixindice(int& id, float theta, float phi);
+	void* calc_pinv_matrices(int nofsources, vector<float> Az, vector<float> El);
+private:
+	int healpixorder;
+	vector<pixelinfo_nc*> pixels_invector;
+	double* legendrekernels;
+	float* legendrekernels_infloat;
+	complex<float>* SRFMultiplier;
+	vector<fmat> pinv_legendrekernels;
+	vector<fmat> res_legendrekernels;
+	float* piv_legendre;
+	double* vonmises_values;
+	float* multiplicationresult;
+	float* pinv_multiplier;
+	double* pinv_multiplier_indouble;
+};
+

NoiseCancellation/include/kernel_noisecancel.cuh

+#include "cuda_runtime.h"
+#include "device_launch_parameters.h"
+
+#include <stdio.h>
+#include <math.h>
+#include <chrono>
+
+
+// dev_a : SRF (AudioBlock x NofFrequencybins x Legendre )
+// dev_c : AudioBlock x NofFrequencybins x Legendre x Legendre
+cudaError_t Device_CovarianceSRF(float* dev_c, float* dev_a, float* dev_a_im, unsigned int freqanalysiswindow, unsigned int BlockSize, unsigned int tfheight, unsigned int freqstart, unsigned int legendresamples);
+// OUT : dev_a : AudioBlock x NofFrequencybins x Nofsource x Legendre
+// IN : dev_c : AudioBlock x NofFrequencybins x Legendre x Legendre
+// IN : dev_wpfcoeff : Nofsource x Legendre
+cudaError_t Device_LeftMultiplier(float* dev_c, float* dev_a, float* dev_wpfcoeff, unsigned int freqanalysiswindow, unsigned int BlockSize, unsigned int tfheight, unsigned int freqstart, unsigned int legendresamples, unsigned int nofsources);
+// OUT : dev_a : AudioBlock x NofFrequencybins x Nofsource x Nofsource
+// IN : dev_c : AudioBlock x NofFrequencybins x Nofsource x Legendre
+// IN : dev_wpfcoeff : Nofsource x Legendre
+cudaError_t Device_RightMultiplier(float* dev_c, float* dev_a, float* dev_wpfcoeff, unsigned int freqanalysiswindow, unsigned int BlockSize, unsigned int tfheight, unsigned int freqstart, unsigned int legendresamples, unsigned int nofsources);
+// IN : dev_tfmultiplier : Nofsource x Nofsource
+// IN : dev_a : AudioBlock x NofFrequencybins x Nofsources
+// OUT : dev_a : AudioBlock x NofFrequencybins x Nofsources
+cudaError_t Device_TFBinMasking(float* dev_a, float* dev_a_im, float* dev_tfmultiplier, float* dev_sumtfmultiplier, unsigned int freqanalysiswindow, unsigned int BlockSize, unsigned int tfheight, unsigned int freqstart, unsigned int legendresamples, unsigned int nofsources);
+// IN : dev_tfmultiplier : AudioBlock x NofFrequencybins x Nofsource x Nofsource
+// IN : dev_tfmultiplier_im : AudioBlock x NofFrequencybins x Nofsource x Nofsource
+// OUT : dev_a : AudioBlock x NofFrequencybins 
+cudaError_t Device_TFBinSummation(float* dev_sum, float* dev_tfmultiplier, unsigned int freqanalysiswindow, unsigned int BlockSize, unsigned int tfheight, unsigned int freqstart, unsigned int legendresamples, unsigned int nofsources);

NoiseCancellation/include/noisecancellation.h

+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <nlohmann/json.hpp>
+#include "kernel_noisecancel.cuh"
+#include "healpixkernel_nc.h"
+
+#define HEALPixOrder 2
+#define MaxSphericalHarmonicsOrder 4
+
+using json = nlohmann::json;
+
+#pragma once
+
+	
+
+struct complex_num_out {
+	const float* real = 0;
+	const float* imag = 0;
+};
+
+class noisecancellation
+{
+
+public:
+	noisecancellation(int windowsize, int blocksize);
+	void getwienermatrices(json audioscenegeometry);
+	complex_num_out Operate_NoiseCancellation(void* srf_multiplier, void* srf_multiplier_im,void* separated_sources, void* separated_sources_im);
+private:
+
+	healpixkernel_nc* hp_kernel;
+	float* magkernels;
+	float* pinv_multiplier = 0;
+
+	unsigned int AudioBlockSize;
+	unsigned int NofFrequency;
+	unsigned int nofsources;
+	unsigned int NofSamples;
+
+	// GPU device malloc parameters
+	const float* dev_pinv_multiplier = 0;
+	const float* dev_covariance = 0;
+	const float* dev_multip = 0;
+	const float* dev_multip_im = 0;
+
+	const float* SRF_multiplier = 0;
+	const float* dev_sum_tfmultiplier = 0;
+};
+

NoiseCancellation/include/pixelinfo_nc.h

+#pragma once
+
+class pixelinfo_nc
+{
+public:
+	int pixelid;
+	float azimuth_in_radian;
+	float elevation_in_radian;
+	float* legendre_kernel_pressure;
+	int neighbour[8];
+	pixelinfo_nc(int pixid,float az,float el,int* neighbours)
+	{
+		this->pixelid = pixid;
+		this->azimuth_in_radian = az;
+		this->elevation_in_radian = el;
+		for (int kjk = 0; kjk < 8; kjk++)
+		{
+			this->neighbour[kjk] = neighbours[kjk];
+		}
+	}
+};

SoundFieldDescription/micarray.h

@@ -6,6 +6,21 @@
 
 #define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
 
+enum class ArrayType
+{
+	Spherical,
+	Circular,
+	Planar,
+	Linear,
+	Other
+};
+
+enum class ArrayScattering
+{
+	Rigid,
+	Open,
+	OtherArray
+};
 
 using namespace std;
 class micarray

SoundFieldDescription/soundfielddescription.h

@@ -26,22 +26,6 @@ using namespace boost::geometry;
 using namespace boost::math;
 using json = nlohmann::json;
 
-enum class ArrayType
-{
-	Spherical,
-	Circular,
-	Planar,
-	Linear,
-	Other
-};
-
-enum class ArrayScattering
-{
-	Rigid,
-	Open,
-	OtherArray
-};
-
 class soundfielddescription
 {
 public:

SpeechDetectionandSeparation/include/healpixkernel_sds.h

+#pragma once
+
+#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string> 
+#include <vector>
+#include <complex>
+
+#include "pixelinfo_sds.h"
+
+#include "boost/geometry.hpp"
+#include "boost/math/special_functions/bessel.hpp"
+#include "boost/math/special_functions/bessel_prime.hpp"
+#include "boost/math/special_functions/legendre.hpp"
+#include "boost/math/special_functions/hankel.hpp"
+#include "boost/math/special_functions/spherical_harmonic.hpp"
+#include "boost/variant/get.hpp"
+
+#include <armadillo>
+
+#define Order0_HP 12
+#define Order1_HP 48
+#define Order2_HP 192
+#define Order3_HP 762
+#define PI 3.14159
+#define Kappa 4.0
+
+#define HPMaxvalue 19.2
+
+
+using namespace std;
+using namespace boost::math;
+using namespace boost::geometry;
+using namespace arma;
+
+class healpixkernel_sds
+{
+public:
+	healpixkernel_sds(int order);
+	void gethealpixcoordinates(int id, float& theta, float& phi);
+	int* gethealpixneighbours(int id);
+	int gethealpixorder();
+	int gethealpixnumberofpixels();
+	float* calculatelegendrekernels(int max_sphorder);
+	float* calculatelegendrekernels_bnkra(int max_sphorder);
+	float* getpinvlegendrekernels(int hpnumber);
+	float* getpinvkernels(int hpnumber);
+	void* getmultiplicationresult(int max_sphorder);
+	void* return_multiplicationresult_real();
+	void* return_multiplicationresult_imag();
+
+	double* create_vonmises(int hpnumber_center);
+	double* create_vonmisesfull();
+	void* calculateorthogonalLK();
+	void* calculateSRFMultiplier(int max_sphorder, int nofmicrophones,int freqrange, vector<float> micanglestheta, vector<float> micanglesphi,float* realpart, float* imagpart);
+	void ConvertSphericalToCartesian(float& distance, float& theta, float& phi, float& x, float& y, float& z);
+	void ConvertCartesionToSpherical(float& distance, float& theta, float& phi, float& x, float& y, float& z);
+	float FindAngleBtw(float& x1, float& y1, float& z1, float& mag1, float& x2, float& y2, float& z2, float& mag2);
+	float def_vonmises(float kappa, float costheta);
+	float* calculatepw_magnitude(int Micnumber, float* micanglestheta, float* micanglesphi);
+	void gethealpixindice(int& id, float theta, float phi);
+	void* calc_pinv_matrices(int nofsources, vector<float> Az, vector<float> El);
+private:
+	int healpixorder;
+	vector<pixelinfo_sds*> pixels_invector;
+	double* legendrekernels;
+	double* legendrekernels_bnkra;
+	double* legendrekernels_bnkra_imag;
+	float* legendrekernels_infloat;
+	float* legendrekernels_infloat_bnkra;
+	float* legendrekernels_infloat_bnkra_imag;
+	complex<float>* SRFMultiplier;
+	vector<fmat> pinv_legendrekernels;
+	vector<fmat> res_legendrekernels;
+	float* piv_legendre;
+	double* vonmises_values;
+	float* multiplicationresult;
+	float* multiplicationresult_imag;
+	float* pinv_multiplier;
+	double* pinv_multiplier_indouble;
+};
+

SpeechDetectionandSeparation/include/kernel.cuh

+
+#include "cuda_runtime.h"
+#include "device_launch_parameters.h"
+
+#include <stdio.h>
+#include <math.h>
+#include <chrono>
+
+
+
+cudaError_t addWithCuda(int* c, const int* a, const int* b, unsigned int size);
+cudaError_t multiplyWithCuda(float* c, float* a, float* b, unsigned int size);
+cudaError_t ComplexmultiplyWithCuda(float* c, float* c_im, float* a, float* b, float* a_im, float* b_im, unsigned int size);
+cudaError_t Device_DComplexmultiplyWithCuda(float* dev_c, float* dev_c_im, float* dev_a, float* dev_b, float* dev_a_im, float* dev_b_im, unsigned int micsize, unsigned int tfheight, unsigned int size, unsigned int BlockSize, unsigned int freqstart, unsigned int legendresamples, float* dev_kr_matrix);
+cudaError_t Device_DmultiplyWithCuda(float* dev_c, float* dev_c_im, float* dev_a, float* dev_b, float* dev_a_im, unsigned int freqanalysiswindow, unsigned int BlockSize, unsigned int tfheight, unsigned int freqstart, unsigned int legendresamples,float* dev_kr_matrix);
+cudaError_t AddWithCuda(float* c, float* c_im, float* a, float* b, float* a_im, float* b_im, unsigned int size);
+cudaError_t SubtractWithCuda(float* c, float* c_im, float* a, float* b, float* a_im, float* b_im, unsigned int size);
+cudaError_t Device_SubtractWithCuda(float* dev_c, float* dev_c_im, float* dev_a, float* dev_b, float* dev_a_im, float* dev_b_im, unsigned int size);
+cudaError_t Device_AddWithCuda(float* dev_c, float* dev_c_im, float* dev_a, float* dev_b, float* dev_a_im, float* dev_b_im, unsigned int size);
+cudaError_t Device_GetHistogram(float* dev_next, float* dev_next_im, float* dev_a, unsigned int freqanalysis, unsigned int legendresamples, unsigned int BlockSize);
+cudaError_t RescomplexmatrixMultiplicationKernel(float* c, float* c_im, float* a, float* b, float* a_im, float* hist, unsigned int freqnumber, unsigned int legendresamples, unsigned int BlockSize, float* dev_kr_matrix);
+cudaError_t MagcomplexmatrixMultiplicationKernel(float* c, float* c_im, float* a, float* b, float* a_im, float* hist, unsigned int noffrequency, unsigned int legendresamples, unsigned int BlockSize, float* dev_kr_matrix);
+cudaError_t Swapmatrices(float* c, float* c_im, float* a, float* a_im, unsigned int noffrequency, unsigned int BlockSize, unsigned int legendresamples);
+cudaError_t Swapmatrices2D(float* c, float* c_im, float* a, float* a_im, unsigned int size, unsigned int BlockSize);
+cudaError_t RENTHistogramKernel(float* c, float* c_im, float* a, float* a_im, float* hist, float* RENThist, unsigned int legendresamples, unsigned int BlockSize, unsigned int noffreqbins);
+cudaError_t CalculateMultiplicationMatrices(float* dev_multiplication_kernel, float* dev_general_histogram, float* dev_kernels_vonmises, float* dev_kernels, unsigned int BlockSize, unsigned int size, unsigned int iterationcount);
+cudaError_t OMPSeparatedSources(float* dev_result, float* dev_result_im, float* dev_multiplication_kernel, float* dev_multiplication_kernel_imag, float* dev_magnitudes, float* dev_magnitudes_im, float* dev_general_histogram, unsigned int BlockSize, unsigned int legendresamples, unsigned int NofFrequencyBins, unsigned int WindowLength, unsigned int nofiterations, float* sourcedoa, unsigned int nofsources,float* dev_pressurelevel,float* dev_krmatrix);

SpeechDetectionandSeparation/include/micarray_sds.h

+#pragma once
+
+#include <string> 
+#include <vector>
+#include <complex>
+
+#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
+
+enum class ArrayType_sds
+{
+	Spherical,
+	Circular,
+	Planar,
+	Linear,
+	Other
+};
+
+enum class ArrayScattering_sds
+{
+	Rigid,
+	Open,
+	OtherArray
+};
+
+using namespace std;
+class micarray_sds
+{
+public:
+	micarray_sds();
+	vector<float> getmicthetas();
+	vector<float> getmicphis();
+	void addnewmicrophone(float mic_x,float mic_y,float mic_z, float mic_ld_x, float mic_ld_y, float mic_ld_z,float ma_x,float ma_y,float ma_z);
+private:
+	void ConvertCartesionToSpherical(float& distance, float& theta, float& phi, float& x, float& y, float& z);
+	void ConvertSphericalToCartesian(float& distance, float& theta, float& phi, float& x, float& y, float& z);
+	int numberofmicrophones;
+	vector<float> microphonetheta;
+	vector<float> microphonephi;
+	vector<float> microphonedistance;
+	vector<float> microphone_x;
+	vector<float> microphone_y;
+	vector<float> microphone_z;
+};
+

SpeechDetectionandSeparation/include/pixelinfo_sds.h

+#pragma once
+
+class pixelinfo_sds
+{
+public:
+	int pixelid;
+	float azimuth_in_radian;
+	float elevation_in_radian;
+	float* legendre_kernel_pressure;
+	int neighbour[8];
+	pixelinfo_sds(int pixid,float az,float el,int* neighbours)
+	{
+		this->pixelid = pixid;
+		this->azimuth_in_radian = az;
+		this->elevation_in_radian = el;
+		for (int kjk = 0; kjk < 8; kjk++)
+		{
+			this->neighbour[kjk] = neighbours[kjk];
+		}
+	}
+};

SpeechDetectionandSeparation/include/resource.h

+//{{NO_DEPENDENCIES}}
+// Microsoft Visual C++ generated include file.
+// Used by speechdetectionandseparation.rc
+
+// Yeni nesneler iin sonraki varsaylan deerler
+// 
+#ifdef APSTUDIO_INVOKED
+#ifndef APSTUDIO_READONLY_SYMBOLS
+#define _APS_NEXT_RESOURCE_VALUE        101
+#define _APS_NEXT_COMMAND_VALUE         40001
+#define _APS_NEXT_CONTROL_VALUE         1001
+#define _APS_NEXT_SYMED_VALUE           101
+#endif
+#endif

SpeechDetectionandSeparation/include/speechdetectionandseparation.h

+
+#include "kernel.cuh"
+
+#include "healpixkernel_sds.h"
+#include "micarray_sds.h"
+
+#define NOMINMAX
+#define _USE_MATH_DEFINES
+#include <cmath>
+
+#include <nlohmann/json.hpp>
+#include <boost/uuid/uuid.hpp>
+#include <boost/uuid/uuid_generators.hpp>
+#include <boost/uuid/uuid_io.hpp>
+
+#define HEALPixOrder 2
+#define MicNumber 32
+#define MaxSphericalHarmonicsOrder 4
+#define AngleDiff 0.79
+#define RENTThreshold 0.5
+#define PI 3.14159
+#define RENTIIteration 3;
+
+using json = nlohmann::json;
+struct RENTData {
+	float* Histogram;
+	float* MaxValuedHistogram;
+	int NofPixels;
+	int PlaneDim;
+	int NofFrequencybinAnalysis;
+};
+
+struct DOA_Data {
+	int hpnumber;
+	float theta;
+	float phi;
+	float distance;
+	float x;
+	float y;
+	float z;
+	bool firstin;
+	int activeframe;
+	string uuid;
+};
+
+struct complex_num {
+	float* real = 0;
+	float* imag = 0;
+	float* pressurelevel = 0;
+};
+
+
+class  speechdetectionandseparation
+{
+public:
+	speechdetectionandseparation(char microphone_array_geometry[],int blocksize);
+	void* initialization(int windowsize);
+	int RENTComputation(float* stft_ch_real, float* stft_ch_imag);
+	complex_num SpeechSeparation();
+	complex_num Return_BlockSRF();
+	json audioscenegeometry_out();
+	void hammingwindow(int windowLength);
+private:
+
+	int* SortInDescendingOrder(int NofPixels, int* inp_mt, float average);
+	int healpix_samplingpoints;
+	float* magkernels;
+	float* magkernels_img;
+	float* multiplicationres;
+	float* multiplicationres_imag;
+	float* kr_1_real;
+	float* kr_1_imag;
+	float* srf_magkernels;
+	float* srf_magkernels_img;
+	vector<float> micanglestheta;
+	vector<float> micanglesphi;
+	float* hammingwin;
+
+
+	healpixkernel_sds* hp_kernel;
+	micarray_sds micarr;
+	vector<DOA_Data> source_doas;
+
+	//json file parameters
+	int nofmicrophones;
+	int nofsamples;
+	ArrayType_sds array_type;
+	ArrayScattering_sds array_scattering_type;
+	int sampling_rate;
+	int sample_type;
+	int block_size;
+	int samplingrates[8] = { 16000,24000,32000,44100,48000,64000,96000,192000 };
+	int sampletypes[3] = { 16,24,32 };
+	int blocksizes[8] = { 64,128,256,512,1024,2048,4096,8192 };
+	int AudioBlockSize;
+	int MaxNumberofSeparatedObjects;
+	int FrequencybinAnalysisStart = 0;
+	int NumberofFrequencybins; // = 1024;
+	int RENT_iteration; // = RENTIIteration;
+	float* max_hist;
+// Block information parameters
+	int BlockIndex;
+	int BlockStart;
+	int BlockEnd;
+
+
+
+	// GPU device malloc parameters
+	const float* dev_a = 0;
+	const float* dev_b = 0;
+	const float* dev_c = 0;
+	const float* dev_kernel = 0;
+	const float* dev_kernels_vonmises = 0;
+	const float* dev_input = 0;
+
+	const float* dev_histogram = 0;
+	const float* dev_RENT_histogram = 0;
+	const float* dev_magnitudes = 0;
+	const float* dev_magnitudes_im = 0;
+	const float* dev_input_im = 0;
+
+	const float* dev_magnitudes_final = 0;
+	const float* dev_magnitudes_im_final = 0;
+
+	const float* dev_general_histogram = 0;
+	const float* dev_multiplication_kernel = 0;
+	const float* dev_multiplication_kernel_imag = 0;
+
+	const float* dev_a_im = 0;
+	const float* dev_b_im = 0;
+	const float* dev_c_im = 0;
+	const float* dev_kernel_im = 0;
+
+	const float* dev_result;
+	const float* dev_result_im;
+
+	const float* dev_res_kernel = 0;
+	const float* dev_pinv_kernel = 0;
+	const float* dev_kr_matrix = 0;
+
+	const float* dev_pressure_level = 0;
+
+	const float* source_doa_hp = 0;
+	const float* dev_hammingwin = 0;
+
+	int DD_InputDataFormat;
+	int PartitionedSamples;
+	int NofFrequencyBins;
+
+	float* RENThist;
+
+	float* c_hist;
+	float* c_hist_im;
+
+	int DD_ResultDataFormat;
+
+	float* actual_result;
+	float* actual_result_im;
+
+	float* host_source_doa_hp;
+
+	int WindowSize;
+	float* kr_matrix;
+
+	float* SRF_real;
+	float* SRF_imag;
+};
+

SynthesisTransform/ReadMe

-###############################Introduction##################################################
-Synthesis tranform is an AIM that is used for transferring the time-frequency domain interleaved microphone array signals into the time domain.
-IN:
-realpart : fft_real[0] to fft_real[N/2 -1]
-imagpart : imagpart[0] to imagpart[N/2 -1]
-OUT:
-Signal : x[0] to x[signalLength - 1]
-
-
-###############################Windows (x64) Installation##################################################
-#include \include directory
--------Debug--------
-libs \libs\Windows_x64\Debug\fftw3.lib 
--------Release--------
-libs \libs\Windows_x64\Release\fftw3.lib   
-
-###############################Windows (x86) Installation##################################################
-#include \include directory
--------Debug--------
-libs 
--------Release--------
-libs  
-
-###############################Linux Installation##################################################
-#include \include directory
--------Debug--------
-libs 
--------Release--------
-libs  
-
-###############################Standalone usage of this AIM (code script)##################################################
-// For analysis of interleaved Nofchannelsx10240 samples with fftsize=2048 and overlap=1024
-// The output will be Interleaved Nofchannelsx20480 in real numbers.
-// Input and output data formats can be followed from MPAI-CAE v1.3 
- 
-synthesistransform* m_fft = new synthesistransform();
-
-float* signal; // Interleaved data output
-float* realpart; // Real part to be malloc
-float* imagpart; // Imaginary part to be malloc
-
-m_fft->IFFT(signal, realpart, imagpart);
-
-

SynthesisTransform/synthesistransform.cpp

@@ -24,7 +24,7 @@ synthesistransform::synthesistransform(int windowsize, int signallength, int hop
     this->overlap = hopsize;
     this->nofmics = nofmics;
 
-    window = (double*)malloc(sizeof(double) * this->windowsize);
+    window = (double*)malloc(sizeof(double) * (this->windowsize+1));
 
     hamming(this->windowsize, window);
 }
@@ -32,26 +32,33 @@ synthesistransform::synthesistransform(int windowsize, int signallength, int hop
 // Create a hamming window of windowLength samples in buffer
 void synthesistransform::hamming(int windowLength, double* buffer) {
 
-    for (int i = 0; i < windowLength; i++) {
+    for (int i = 0; i <= windowLength; i++) {
 
-        buffer[i] = 0.54 - (0.46 * cos(2 * PI * (i / (((double)windowLength - 1.0) * 1.0))));
+        buffer[i] = (0.53836 - (0.46164 * cos(2 * PI * ((double)i / (((double)windowLength) * 1.0))))) ;
     }
 }
 
+void synthesistransform::SetSourceSize(int nofsource)
+{
+
+    this->nofmics = nofsource;
+
+}
+
 
 void* synthesistransform::IFFT(double* istft_out, double* realpart, double* imagpart) {
 
 
-    fftw_complex* data, * ifft_output;
+    complex<double>* data, * ifft_output;
     fftw_plan       plan_backward;
     int             i;
 
-    data = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * this->windowsize);
-    ifft_output = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * this->windowsize);
+    data = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);
+    ifft_output = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);
 
     int outputcounter = 0;
 
-    plan_backward = fftw_plan_dft_1d(this->windowsize, data, ifft_output, FFTW_BACKWARD, FFTW_ESTIMATE);
+    plan_backward = fftw_plan_dft_1d(this->windowsize, reinterpret_cast<fftw_complex*>(data), reinterpret_cast<fftw_complex*>(ifft_output), FFTW_BACKWARD, FFTW_ESTIMATE);
 
     int readIndex;
 
@@ -61,6 +68,7 @@ void* synthesistransform::IFFT(double* istft_out, double* realpart, double* imag
     // Process each chunk of the signal
     int channelindex = 0;
 
+
     for (channelindex = 0; channelindex < this->nofmics; channelindex++)
     {
         int outputcounter = channelindex;
@@ -72,17 +80,17 @@ void* synthesistransform::IFFT(double* istft_out, double* realpart, double* imag
             {
                 readIndex = (chunkPosition + i) * this->nofmics + channelindex;
 
-                data[i][0] = *(realpart + readIndex);
-                data[i][1] = *(imagpart + readIndex);
+                reinterpret_cast<fftw_complex*>(data)[i][0] = *(realpart + readIndex);
+                reinterpret_cast<fftw_complex*>(data)[i][1] = *(imagpart + readIndex); 
 
             }
 
             for (i = this->windowsize / 2; i < this->windowsize; i++)
             {
-                readIndex = (chunkPosition + (i - this->windowsize / 2)) * this->nofmics + channelindex;
+                //readIndex = (chunkPosition + (i - this->windowsize / 2)) * this->nofmics + channelindex;
 
-                data[i][0] = 0.0; // *(realpart + readIndex);
-                data[i][1] = 0.0; // *(imagpart + readIndex);
+                reinterpret_cast<fftw_complex*>(data)[i][0] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][0]; // *(realpart + readIndex);
+                reinterpret_cast<fftw_complex*>(data)[i][1] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][1]; // *(imagpart + readIndex);
 
             }
             // Perform the FFT on our chunk
@@ -97,12 +105,118 @@ void* synthesistransform::IFFT(double* istft_out, double* realpart, double* imag
 
             for (i = 0; i < this->windowsize/2; i++)
             {
-                istft_out[outputcounter] = (ifft_output[i][0] / window[i]) / (this->windowsize / 2); // / windowSize; // The i term should be non-existent.
+                 
+                istft_out[outputcounter] = (double)(ifft_output[i].real());/* / (double)window[i]*/; // / (double)((this->windowsize / 2.0)); // / windowSize; // The i term should be non-existent.
+
+                outputcounter += this->nofmics;
+            }
+
+            chunkPosition += this->overlap;
+
+        } // Excuse the formatting, the while ends here.
+
+        bStop = 0;
+    }
+
+    fftw_destroy_plan(plan_backward);
+
+    fftw_free(data);
+    fftw_free(ifft_output);
+    return 0;
+
+}
+
+void* synthesistransform::IFFT(double* istft_out, float* realpart, float* imagpart) {
+
+
+    complex<double>* data, * ifft_output;
+    fftw_plan       plan_backward;
+    int             i;
+
+    data = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);
+    ifft_output = (complex<double>*)fftw_malloc(sizeof(complex<double>) * this->windowsize);
+
+    int outputcounter = 0;
+
+    plan_backward = fftw_plan_dft_1d(this->windowsize, reinterpret_cast<fftw_complex*>(data), reinterpret_cast<fftw_complex*>(ifft_output), FFTW_BACKWARD, FFTW_ESTIMATE);
+
+    int readIndex;
+
+    // Should we stop reading in chunks?
+    int bStop = 0;
+
+    // Process each chunk of the signal
+    int channelindex = 0;
+    int zeropadding = 0;
+
+    for (channelindex = 0; channelindex < this->nofmics; channelindex++)
+    {
+        int outputcounter = channelindex;
+        int chunkPosition = 0;
+
+       /* for (int ki = 0; ki < (this->signallength * this->nofmics); ki++)
+        {
+            istft_out[ki] = 0.0;
+        }*/
+
+        //float pressurelev = 0.00001 / pressurelevel[channelindex];
+
+        while (chunkPosition < this->signallength && !bStop) {
+
+            for (i = 0; i < zeropadding; i++)
+            {
+                //readIndex = (chunkPosition + (i - this->windowsize / 2)) * this->nofmics + channelindex;
+
+                reinterpret_cast<fftw_complex*>(data)[i][0] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][0]; // *(realpart + readIndex);
+                reinterpret_cast<fftw_complex*>(data)[i][1] = 0.0; // reinterpret_cast<fftw_complex*>(data)[this->windowsize - i - 1][1]; // *(imagpart + readIndex);
+            }
 
+            for (i = zeropadding; i < this->windowsize / 2; i++)
+            {
+                readIndex = chunkPosition  * this->nofmics + i * this->nofmics + channelindex;
+                // If there is not an analysis use directly
+                reinterpret_cast<fftw_complex*>(data)[i][0] = *(realpart + readIndex); // *window[i]; // *(float)((this->windowsize));
+                reinterpret_cast<fftw_complex*>(data)[i][1] = *(imagpart + readIndex); // *window[i]; // *(float)((this->windowsize));
+                
+//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
+              //reinterpret_cast<fftw_complex*>(data)[i][0] = *(realpart + readIndex);
+              //reinterpret_cast<fftw_complex*>(data)[i][1] = *(imagpart + readIndex);
+//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
+
+            }
+            for (i = (this->windowsize / 2); i < this->windowsize; i++)
+            {
+                readIndex = (chunkPosition + (this->windowsize -i - 1)) * this->nofmics + channelindex;
+
+                reinterpret_cast<fftw_complex*>(data)[i][0] = 0.0; // *(realpart + readIndex)* window[(this->windowsize / 2) - i];
+                reinterpret_cast<fftw_complex*>(data)[i][1] = 0.0; //  *(imagpart + readIndex)* (-1) * window[this->windowsize - i];
+
+            }
+            // Perform the FFT on our chunk
+            fftw_execute(plan_backward);
+
+            // Copy the first (windowSize/2 + 1) data points into your spectrogram.
+            // We do this because the FFT output is mirrored about the nyquist
+            // frequency, so the second half of the data is redundant. This is how
+            // Matlab's spectrogram routine works.
+
+            //outputcounter = chunkPosition * this->nofmics + channelindex;
+
+
+            
+            for (i = 0; i < this->windowsize  ; i++)
+            {
+                // If there is not an analysis use directly
+                istft_out[outputcounter] = (double)(ifft_output[i].real()); // / windowSize; // The i term should be non-existent.
+//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
+                // Without analysis use this one:
+              //  istft_out[outputcounter] = ((double)(ifft_output[i].real()) / window[i]); // / (double)((this->windowsize));
+//USE THIS ONE !/////////////////////////////////////////////////////////////////////////////////////////
                 outputcounter += this->nofmics;
             }
 
             chunkPosition += this->overlap;
+            outputcounter -= this->overlap * this->nofmics;
 
         } // Excuse the formatting, the while ends here.
 

SynthesisTransform/synthesistransform.h

@@ -15,6 +15,8 @@ public:
     synthesistransform();
     synthesistransform(int windowsize, int signallength, int hopsize, int nofmics);
     void* IFFT(double* signal, double* realpart, double* imagpart);
+    void* IFFT(double* istft_out, float* realpart, float* imagpart);
+    void SetSourceSize(int nofsource);
 private:
     void hamming(int windowLength, double* buffer);
     int signallength;

config/0_healpix_res.txt

+0,0.8410686705679303,0.7853981633974483,4,-1,3,2,1,-1,5,8
+1,0.8410686705679303,2.356194490192345,5,-1,0,3,2,-1,6,9
+2,0.8410686705679303,3.926990816987241,6,-1,1,0,3,-1,7,10
+3,0.8410686705679303,5.497787143782138,7,-1,2,1,0,-1,4,11
+4,1.5707963267948966,0.0,11,7,3,-1,0,5,8,-1
+5,1.5707963267948966,1.5707963267948966,8,4,0,-1,1,6,9,-1
+6,1.5707963267948966,3.141592653589793,9,5,1,-1,2,7,10,-1
+7,1.5707963267948966,4.71238898038469,10,6,2,-1,3,4,11,-1
+8,2.300523983021863,0.7853981633974483,11,-1,4,0,5,-1,9,10
+9,2.300523983021863,2.356194490192345,8,-1,5,1,6,-1,10,11
+10,2.300523983021863,3.926990816987241,9,-1,6,2,7,-1,11,8
+11,2.300523983021863,5.497787143782138,10,-1,7,3,4,-1,8,9