Add user sources from external wave files

Author: skhelladi

CMakeLists.txt

@@ -88,10 +88,16 @@ endif()
 FIND_PACKAGE(VTK COMPONENTS IOXML REQUIRED)
 SET(VTK_DIR "/usr/lib/x86_64-linux-gnu/cmake/vtk-9.1" CACHE PATH "VTK directory override" FORCE)
 
+SET(WAVE_INCLUDE src/wave)
+INCLUDE_DIRECTORIES(${WAVE_INCLUDE})
 
+SET(WAVE_LD ${PROJECT_SOURCE_DIR}/build/bin)
+LINK_DIRECTORIES(${WAVE_LD})
 
 # Libraries
 INCLUDE_DIRECTORIES(include)
 
 # Sources
+ADD_SUBDIRECTORY(src/wave)
 ADD_SUBDIRECTORY(src)
+

README.md

@@ -113,9 +113,10 @@ c0 = 100
 
 # Source:
 # name = fct,x,y,z,size,intensity,frequency,phase,duration
-# - fct supported = [monopole, dipole, quadrupole, formula, file (coming soon)]
+# - fct supported = [monopole, dipole, quadrupole, formula, file (csv, wav)]
 # - if fct = formula => name = fct,"formula expr",x,y,z,size,duration (ex: formulat = 0.1 * sin(2 * pi * 50 * t))
 # - if fct = file => name = fct,"filename",x,y,z,size
+#	suported file formats are : csv, wav
 # - (x,y,z) = source position
 # - intensity = source intensity
 # - frequency = source frequency
@@ -124,6 +125,7 @@ c0 = 100
 source1 = formula, "0.1 * sin(2 * pi * 50 * t)", 0.0,0.0,0.0, 0.1, 0.1
 source2 = monopole, 0.0,0.0,0.0, 0.1, 0.1,50,0,0.1
 source3 = file,"data/data.csv", 0.0,0.0,0.0, 0.1
+source4 = file,"data/data.wav", 0.0,0.0,0.0, 0.1
 # source4 = udf, "-0.1 * sin(2 * pi * 50 * t)", -0.5,0.0,0.0, 0.1, 0.1
 
 # Initial condition:

data/data1.wav

@@ -7,8 +7,8 @@
 
 # Initial, final time and time step(t>0)
 timeStart=0
-timeEnd=0.075
-timeStep=0.00001
+timeEnd=1.0
+timeStep=0.00005
 
 # Saving rate:
 timeRate=0.001
@@ -39,25 +39,24 @@ c0 = 100
 
 # Source:
 # name = fct,x,y,z,size,intensity,frequency,phase,duration
-# - fct supported = [monopole, dipole, quadrupole, formula, file (coming soon)]
-# - if fct = formula => name = fct,"formula expr",x,y,z,size,duration (ex: formula = 0.1 * sin(2 * pi * 50 * t))
+# - fct supported = [monopole, dipole, quadrupole, formula, file (csv, wav)]
+# - if fct = formula => name = fct,"formula expr",x,y,z,size,duration (ex: formulat = 0.1 * sin(2 * pi * 50 * t))
 # - if fct = file => name = fct,"filename",x,y,z,size
+#	suported file formats are : csv, wav
 # - (x,y,z) = source position
 # - intensity = source intensity
 # - frequency = source frequency
 # NB: Extended source or Multiple sources are supported.
 #     (source1 = ..., source2 = ...) indice must change.
-# source1 = monopole, 0.0,0.0,0.0, 0.1, 0.1,50,0,0.1
-# source1 = formula, "0.1 * sin(2 * pi * 50 * t)", 0.0,0.0,0.0, 0.1, 0.2
 source1 = file,"data/data.csv", 0.0,0.0,0.0, 0.1
-# source2 = udf, "-0.1 * sin(2 * pi * 50 * t)", -0.5,0.0,0.0, 0.1, 0.1
+# source2 = file,"data/data2.wav", 0.0,0.0,0.0, 0.1
 
 # Observers position:
 # name = x,y,z, size
 # - (x,y,z) = position
 # NB: Multiple observers are supported and recursively added.
 #     (observer1 = ..., observer2 = ...)
-# observer1 = 2.11792,0.00340081,0.0,0.1
+observer1 = 2.11792,0.00340081,0.0,0.1
 # observer2 = -2.11792,0.00340081,0.0,0.1
 
 # Initial condition:

data/data2.wav

@@ -1,10 +1,11 @@
-#include <configParser.h>
 #include <Mesh.h>
+#include <configParser.h>
 
 #ifndef DGALERKIN_SOLVER_H
 #define DGALERKIN_SOLVER_H
 
-namespace solver {
+namespace solver
+{
     /**
      * Solve using forward explicit scheme. O(h)
      *
@@ -22,6 +23,9 @@ namespace solver {
      * @param config
      */
     void rungeKutta(std::vector<std::vector<double>> &u, Mesh &mesh, Config config);
+
+    // std::vector<std::vector<float>> data4wave;
+
 }
 
 #endif

data/data3.wav

@@ -5,6 +5,11 @@
 #include <iomanip>
 #include <iostream>
 #include <vector>
+
+#include "file.h"
+#include "err.h"
+
+
 // #include<execution>
 // #include <assert.h>
 
@@ -37,6 +42,8 @@
 //////////////////////////////////////////////////////////////
 // using namespace std;
 
+std::string fileExtension(std::string file);
+
 namespace screen_display
 {
     void write_string(std::string text, std::string color = YELLOW);
@@ -88,7 +95,9 @@ namespace io
      * @return data as vector of vector of doubles.
      */
     std::vector<std::vector<double>> parseCSVFile(std::string inputFileName, char separator);
-
+    void writeWave(std::vector<float> V, std::string filename, uint32_t sample_rate, uint16_t bits_per_sample=16, uint16_t channel_number=1, size_t nb_sequence=1);
+    void readWave(std::string filename, std::vector<float> &V, uint32_t &sample_rate);
+    std::vector<std::vector<double>> parseWAVEFile(std::string inputFileName);
 }
 /////////////////////////
 

data/data4.wav

@@ -1,10 +1,8 @@
 # programs
-
 SET(SRCS dgalerkin.cpp configParser.cpp ../include/configParser.h Mesh.cpp ../include/Mesh.h utils.cpp ../include/utils.h solver.cpp ../include/solver.h ../include/eqEdit.h eqEdit.cpp)
 
 ADD_EXECUTABLE(dgalerkin ${SRCS})
-
-TARGET_LINK_LIBRARIES(dgalerkin ${GMSH_LIBRARIES} ${LAPACKBLAS_LIBRARIES} ${VTK_LIBRARIES} -fopenmp -ltbb -lpthread)
+TARGET_LINK_LIBRARIES(dgalerkin ${GMSH_LIBRARIES} ${LAPACKBLAS_LIBRARIES} ${VTK_LIBRARIES} -fopenmp -ltbb -lpthread -lwave)
 
 vtk_module_autoinit(
 	TARGETS dgalerkin

doc/config/Square.conf

@@ -73,7 +73,7 @@ namespace config
                 if (key.find("source") == 0)
                 {
                     std::vector<std::string> sep = split(iter->second, ',');
-                    
+
                     if (sep[0] != "formula" && sep[0] != "file")
                     {
                         double x = std::stod(sep[1]);
@@ -116,7 +116,7 @@ namespace config
                         }
                     }
                     else
-                    { 
+                    {
 
                         if (sep[0] == "formula")
                         {
@@ -148,9 +148,23 @@ namespace config
                             double z = std::stod(sep[4]);
                             double size = std::stod(sep[5]);
                             // double duration = std::stod(sep[6]);
-                            // double pole = -1;    
-                            Sources S("", {-1, x, y, z, size},io::parseCSVFile(filename,';'));
-                            config.sources.push_back(S);
+                            // double pole = -1;
+                            if (fileExtension(filename) == "csv")
+                            {
+                                Sources S("", {-1, x, y, z, size}, io::parseCSVFile(filename, ';'));
+                                config.sources.push_back(S);
+                            }
+                            else if (fileExtension(filename) == "wav")
+                            {
+                                Sources S("", {-1, x, y, z, size}, io::parseWAVEFile(filename));
+                                config.sources.push_back(S);
+                            }
+                            else
+                            {
+                                Fatal_Error("Not supported data")
+                            }
+
+                            // getchar();
                         }
                     }
                 }

include/solver.h

@@ -23,6 +23,8 @@ namespace solver
     std::vector<double> elStiffvector;
     std::vector<std::vector<std::vector<double>>> Flux;
 
+    std::vector<std::vector<float>> data4wave;
+
     /**
      * Perform a numerical step: u[t+1] = dt*M^-1*(S[u[t]]-F[u[t]]) + beta*u[t]
      * for all elements in mesh object.
@@ -139,6 +141,8 @@ namespace solver
         outfile << "time;res_p;res_rho;res_vx;res_vy;res_vz;elapsed_time" << std::endl;
 
         std::vector<std::ofstream> obs_outfile(config.observers.size());
+        data4wave.clear();
+        data4wave.resize(config.observers.size());
         for (int obs = 0; obs < config.observers.size(); ++obs)
         {
             std::string filename = "results/observers" + std::to_string(obs + 1) + ".txt";
@@ -192,7 +196,7 @@ namespace solver
 
             for (int src = 0; src < config.sources.size(); ++src)
             {
-                if (config.sources[src].formula == "")
+                if (config.sources[src].formula == "" && config.sources[src].data.empty())
                 {
                     double amp = config.sources[src].source[5];
                     double freq = config.sources[src].source[6];
@@ -204,10 +208,18 @@ namespace solver
                 }
                 else
                 {
-                    double duration = config.sources[src].source[5];
-                    if (t < duration)
+                    if (config.sources[src].data.empty())
+                    {
+                        double duration = config.sources[src].source[5];
+                        if (t < duration)
+                            for (int n = 0; n < srcIndices[src].size(); ++n)
+                                u[0][srcIndices[src][n]] = config.sources[src].value(t);
+                    }
+                    else
+                    {
                         for (int n = 0; n < srcIndices[src].size(); ++n)
-                            u[0][srcIndices[src][n]] = config.sources[src].value(t);
+                            u[0][srcIndices[src][n]] = config.sources[src].interpolate_value(t);
+                    }
                 }
             }
 
@@ -274,9 +286,13 @@ namespace solver
                 v[0] /= w_sum;
                 v[1] /= w_sum;
                 v[2] /= w_sum;
+                data4wave[obs].push_back(p);
                 obs_outfile[obs] << t << ";" << rho << ";" << p << ";" << v[0] << ";" << v[1] << ";" << v[2] << std::endl;
             }
         }
+        for (int obs = 0; obs < config.observers.size(); ++obs)
+            io::writeWave(data4wave[obs], "results/observer_" + std::to_string(obs + 1) + ".wav", 1.0 / config.timeStep, 16, 1, 1);
+
         outfile.close();
         for (int obs = 0; obs < config.observers.size(); ++obs)
             obs_outfile[obs].close();
@@ -381,7 +397,8 @@ namespace solver
         outfile << "time;log(res_p);log(res_rho);log(res_vx);log(res_vy);log(res_vz);elapsed_time(s)" << std::endl;
 
         std::vector<std::ofstream> obs_outfile(config.observers.size());
-
+        data4wave.clear();
+        data4wave.resize(config.observers.size());
         for (int obs = 0; obs < config.observers.size(); ++obs)
         {
             std::string filename = "results/observers" + std::to_string(obs + 1) + ".txt";
@@ -453,10 +470,11 @@ namespace solver
                         if (t < duration)
                             for (int n = 0; n < srcIndices[src].size(); ++n)
                                 u[0][srcIndices[src][n]] = config.sources[src].value(t);
-                    }else 
+                    }
+                    else
                     {
-                       for (int n = 0; n < srcIndices[src].size(); ++n)
-                                u[0][srcIndices[src][n]] = config.sources[src].interpolate_value(t); 
+                        for (int n = 0; n < srcIndices[src].size(); ++n)
+                            u[0][srcIndices[src][n]] = config.sources[src].interpolate_value(t);
                     }
                 }
             }
@@ -547,9 +565,13 @@ namespace solver
                 v[0] /= w_sum;
                 v[1] /= w_sum;
                 v[2] /= w_sum;
+                data4wave[obs].push_back(p);
                 obs_outfile[obs] << t << ";" << rho << ";" << p << ";" << v[0] << ";" << v[1] << ";" << v[2] << std::endl;
             }
         }
+        for (int obs = 0; obs < config.observers.size(); ++obs)
+            io::writeWave(data4wave[obs], "results/observer_" + std::to_string(obs + 1) + ".wav", 1.0 / config.timeStep, 16, 1, 1);
+
         outfile.close();
         for (int obs = 0; obs < config.observers.size(); ++obs)
             obs_outfile[obs].close();

include/utils.h

@@ -28,6 +28,14 @@ extern "C"
 
 //! added by Sofiane KHELLADI in 11/03/2022 /////////////////////
 
+std::string fileExtension(std::string file)
+{
+
+    std::size_t found = file.find_last_of(".");
+    return file.substr(found + 1);
+}
+
+
 namespace screen_display
 {
     void write_string(std::string text, std::string color)
@@ -168,6 +176,83 @@ namespace io
 
         return data;
     }
+
+    void writeWave(std::vector<float> V, std::string filename, uint32_t sample_rate, uint16_t bits_per_sample, uint16_t channel_number, size_t nb_sequence)
+    {
+        wave::File write_file;
+        // uint32_t sample_rate = 10000;//V.size();
+        // uint16_t bits_per_sample = 16;
+        // uint16_t channel_number = 1;
+        // filename += ".wav";
+
+        std::vector<float> content;
+
+        for (size_t i = 0; i < nb_sequence; i++)
+            for (auto v : V)
+            {
+                content.push_back(v);
+            }
+
+        wave::Error err = write_file.Open(filename, wave::kOut);
+        if (err)
+        {
+            Fatal_Error("Something went wrong in out open")
+        }
+        write_file.set_sample_rate(sample_rate);
+        write_file.set_bits_per_sample(bits_per_sample);
+        write_file.set_channel_number(channel_number);
+
+        err = write_file.Write(content);
+        if (err)
+        {
+            Fatal_Error("Something went wrong in write")
+        }
+    }
+
+    void readWave(std::string filename, std::vector<float> &V, uint32_t &sample_rate)
+    {
+        wave::File read_file;
+        wave::Error err = read_file.Open(filename, wave::kIn);
+        if (err)
+        {
+            Fatal_Error("Something went wrong in out open")
+        }
+
+        // ASSERT_EQ(read_file.sample_rate(), 44100);
+        // ASSERT_EQ(read_file.bits_per_sample(), 16);
+        // ASSERT_EQ(read_file.channel_number(), 2);
+        sample_rate = read_file.sample_rate();
+        err = read_file.Read(&V);
+        if (err)
+        {
+            Fatal_Error("Something went wrong in read")
+        }
+
+        std::cout << "source file="<< filename << std::endl;
+        std::cout << "sample_rate=" << read_file.sample_rate() << std::endl;
+        std::cout << "bits_per_sample=" << read_file.bits_per_sample() << std::endl;
+        std::cout << "channel_number=" << read_file.channel_number() << std::endl;
+    }
+
+    std::vector<std::vector<double>> parseWAVEFile(std::string inputFileName)
+    {
+        std::vector<std::vector<double>> data;
+        uint32_t sample_rate;
+        std::vector<float> V;
+
+        readWave(inputFileName,V,sample_rate);
+        double timeStep = 1.0/sample_rate;
+        double time(0);
+        for(auto v:V)
+        {
+          data.push_back({time,v});
+          time += timeStep;
+        }
+
+        // writeWave(V,"data/data_write.wav",sample_rate,16,1,1);
+
+        return data;
+    }
 }
 
 // Lapack with direct calling to fortran interface.