diff --git a/src/rvoice/fluid_chorus.c b/src/rvoice/fluid_chorus.c
index c80dc9e98..a1168a7c6 100644
--- a/src/rvoice/fluid_chorus.c
+++ b/src/rvoice/fluid_chorus.c
@@ -169,10 +169,11 @@
 /* modulator */
 typedef struct
 {
-    fluid_real_t   a1;          /* Coefficient: a1 = 2 * cos(w) */
-    fluid_real_t   buffer1;     /* buffer1 */
-    fluid_real_t   buffer2;     /* buffer2 */
-    fluid_real_t   reset_buffer2;/* reset value of buffer2 */
+    // for sufficient precision members MUST be double! See https://github.com/FluidSynth/fluidsynth/issues/1331
+    double   a1;           /* Coefficient: a1 = 2 * cos(w) */
+    double   buffer1;      /* buffer1 */
+    double   buffer2;      /* buffer2 */
+    double   reset_buffer2;/* reset value of buffer2 */
 } sinus_modulator;
 
 /*-----------------------------------------------------------------------------
@@ -236,6 +237,10 @@ struct _fluid_chorus_t
 /*-----------------------------------------------------------------------------
  Sets the frequency of sinus oscillator.
 
+ For sufficient precision use double precision in set_sinus_frequency() computation !.
+ Never use: fluid_real_t , cosf(), sinf(), FLUID_COS(), FLUID_SIN(), FLUID_M_PI.
+ See https://github.com/FluidSynth/fluidsynth/issues/1331
+
  @param mod pointer on modulator structure.
  @param freq frequency of the oscillator in Hz.
  @param sample_rate sample rate on audio output in Hz.
@@ -244,16 +249,17 @@ struct _fluid_chorus_t
 static void set_sinus_frequency(sinus_modulator *mod,
                                 float freq, float sample_rate, float phase)
 {
-    fluid_real_t w = 2 * FLUID_M_PI * freq / sample_rate; /* initial angle */
-    fluid_real_t a;
+    double w = (2.0 * M_PI) * freq / sample_rate;  /* step phase between each sinus wave sample (in radian) */
+    double a; /* initial phase at which the sinus wave must begin (in radian) */
 
-    mod->a1 = 2 * FLUID_COS(w);
+    // DO NOT use potentially single precision cosf or FLUID_COS here! See https://github.com/FluidSynth/fluidsynth/issues/1331
+    mod->a1 = 2 * cos(w);
 
-    a = (2 * FLUID_M_PI / 360) * phase;
+    a = (2.0 * M_PI / 360.0) * phase;
 
-    mod->buffer2 = FLUID_SIN(a - w); /* y(n-1) = sin(-initial angle) */
-    mod->buffer1 = FLUID_SIN(a); /* y(n) = sin(initial phase) */
-    mod->reset_buffer2 = FLUID_SIN(FLUID_M_PI / 2 - w); /* reset value for PI/2 */
+    mod->buffer2 = sin(a - w); /* y(n-1) = sin(-initial angle) */
+    mod->buffer1 = sin(a); /* y(n) = sin(initial phase) */
+    mod->reset_buffer2 = sin((M_PI / 2.0) - w); /* reset value for PI/2 */
 }
 
 /*-----------------------------------------------------------------------------
@@ -264,21 +270,21 @@ static void set_sinus_frequency(sinus_modulator *mod,
  @param mod pointer on modulator structure.
  @return current value of the modulator sine wave.
 -----------------------------------------------------------------------------*/
-static FLUID_INLINE fluid_real_t get_mod_sinus(sinus_modulator *mod)
+static FLUID_INLINE double get_mod_sinus(sinus_modulator *mod)
 {
-    fluid_real_t out;
+    double out;
     out = mod->a1 * mod->buffer1 - mod->buffer2;
     mod->buffer2 = mod->buffer1;
 
-    if(out >= 1.0f) /* reset in case of instability near PI/2 */
+    if(out >= 1.0) /* reset in case of instability near PI/2 */
     {
-        out = 1.0f; /* forces output to the right value */
+        out = 1.0; /* forces output to the right value */
         mod->buffer2 = mod->reset_buffer2;
     }
 
-    if(out <= -1.0f) /* reset in case of instability near -PI/2 */
+    if(out <= -1.0) /* reset in case of instability near -PI/2 */
     {
-        out = -1.0f; /* forces output to the right value */
+        out = -1.0; /* forces output to the right value */
         mod->buffer2 = - mod->reset_buffer2;
     }