* Updating synth.c with a Gaussian white noise effect?
@ 2018-10-08 21:59 James Trammell
0 siblings, 0 replies; only message in thread
From: James Trammell @ 2018-10-08 21:59 UTC (permalink / raw)
To: sox-devel
[-- Attachment #1.1: Type: text/plain, Size: 1274 bytes --]
Hi,
I have tried to write my own SoX effect to generate Gaussian white noise as
an alternative to the uniform white noise that it generates currently. Even
though I have a copy of "Numerical Recipes in C", a fair amount of C
experience, and the internet at my disposal, I'm still stuck.
Would anyone mind looking at the code I have come up with so far, and
hopefully tell me where I am screwing up?
A patch file for synth.c is attached; patch your v14.4.2 synth.c file with
my changes using:
cd sox-14.4.2/src
patch -p1 <synth.diff
(If you would like me to send you my modified synth.c itself, let me know.)
I tried a number of ways to accomplish this; that is why there are so many
commented out blocks. The last block of code is where I left off.
The new effect is invoked like this:
sox -V -b 24 -r 96000 -n gaussianwhite.wav synth 60 gaussianwhitenoise
gaussianwhitenoise gain -12
Thanks,
James
case synth_gaussianwhitenoise: {
static double V1, V2, S;
static int phase = 0;
double X;
do {
if(phase == 0) {
do {
V1 = DRANQD1;
V2 = DRANQD1;
S = V1 * V1 + V2 * V2;
} while (S >= 1 || S == 0);
X = V1 * sqrt(-2 * log(S) / S);
} else {
X = V2 * sqrt(-2 * log(S) / S);
}
phase = 1 - phase;
} while (X < -1 || X > 1);
synth_out = X;
//printf("%lf\n", synth_out);
}
[-- Attachment #1.2: Type: text/html, Size: 3692 bytes --]
[-- Attachment #2: synth.diff --]
[-- Type: application/octet-stream, Size: 11320 bytes --]
--- /Users/jtrammel/Developer/sox/sox-14.4.2-virgin/src/synth.c 2014-10-05 22:02:30.000000000 -0400
+++ /Users/jtrammel/Developer/sox/sox-14.4.2/src/synth.c 2018-10-05 13:13:26.000000000 -0400
@@ -46,7 +46,8 @@
synth_tpdfnoise,
synth_pinknoise,
synth_brownnoise,
- synth_pluck
+ synth_pluck,
+ synth_gaussianwhitenoise
} type_t;
static lsx_enum_item const synth_type[] = {
@@ -63,6 +64,7 @@
LSX_ENUM_ITEM(synth_, pinknoise)
LSX_ENUM_ITEM(synth_, brownnoise)
LSX_ENUM_ITEM(synth_, pluck)
+ LSX_ENUM_ITEM(synth_, gaussianwhitenoise)
{0, 0}
};
@@ -468,7 +470,7 @@
elapsed_time_s;
break;
case Square:
- phase = (chan->freq + sign(chan->mult) *
+ phase = (chan->freq + sign(chan->mult) *
sqr(p->samples_done * chan->mult)) * elapsed_time_s;
break;
case Exp:
@@ -571,9 +573,622 @@
default: synth_out = 0;
}
} else switch (chan->type) {
- case synth_whitenoise:
+
+
+/*
+ *
+ *
+ * ATTEMPT #1
+ *
+ *
+ *
+ // http://c-faq.com/lib/gaussian.html
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+
+ if(phase == 0) {
+ do {
+ double U1 = DRANQD1 / RAND_MAX;
+ double U2 = DRANQD1 / RAND_MAX;
+
+ V1 = 2 * U1 - 1;
+ V2 = 2 * U2 - 1;
+ S = V1 * V1 + V2 * V2;
+ } while(S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+ } else {
+ X = V2 * sqrt(-2 * log(S) / S);
+ }
+
+ phase = 1 - phase;
+
+ synth_out = X;
+ break;
+
+ }
+
+*/
+
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #2
+ *
+ *
+ *
+
+ // http://c-faq.com/lib/gaussian.html
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+
+ if(phase == 0) {
+ do {
+
+ V1 = 2 * DRANQD1 - 1;
+ V2 = 2 * DRANQD1 - 1;
+ S = V1 * V1 + V2 * V2;
+ } while(S >= 1 || S == 0);
+
+ synth_out = V1 * sqrt(-2 * log(S) / S);
+ } else {
+ synth_out = V2 * sqrt(-2 * log(S) / S);
+ }
+
+ phase = 1 - phase;
+
+ break;
+
+ }
+
+*/
+
+
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #3
+ *
+ *
+ *
+
+ // http://c-faq.com/lib/gaussian.html
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+
+ if(phase == 0) {
+ do {
+
+ double U1 = (double)DRANQD1 / RAND_MAX;
+ double U2 = (double)DRANQD1 / RAND_MAX;
+
+ V1 = 2 * U1 - 1;
+ V2 = 2 * U2 - 1;
+ S = V1 * V1 + V2 * V2;
+ } while(S >= 1 || S == 0);
+
+ synth_out = V1 * sqrt(-2 * log(S) / S);
+
+ } else {
+
+ synth_out = V2 * sqrt(-2 * log(S) / S);
+
+ }
+
+ phase = 1 - phase;
+
+ break;
+
+ }
+
+*/
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #4
+ *
+ *
+ *
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+
+ if(phase == 0) {
+ do {
+
+ double U1 = DRANQD1;
+ double U2 = DRANQD1;
+
+ V1 = 2 * U1 - 1;
+ V2 = 2 * U2 - 1;
+
+ S = V1 * V1 + V2 * V2;
+
+ } while(S >= 1 || S <= -1);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+
+ } else {
+
+ X = V2 * sqrt(-2 * log(S) / S);
+
+ }
+
+ phase = 1 - phase;
+
+ synth_out = X;
+ printf("%lf\n", synth_out);
+
+ break;
+
+ }
+
+*/
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #5
+ *
+ *
+ *
+ case synth_gaussianwhitenoise: {
+
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+
+ if (phase == 0) {
+
+ do {
+
+ double U1 = DRANQD1;
+ double U2 = DRANQD1;
+
+ S = U1 * U1 + U2 * U2;
+
+ } while (S >= 1 || S == 0);
+
+ V1 = U1 * sqrt((-2.0 * log(S)) / S);
+ double ln1 = exp(V1);
+ synth_out = ln1;
+
+ } else {
+
+ V2 = U2 * sqrt((-2.0 * log(S)) / S);
+ double ln2 = exp(V2);
+ synth_out = ln2;
+ }
+
+ phase = 1 - phase;
+ printf("%lf\n", synth_out);
+
+
+ break;
+
+ }
+
+*/
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #6
+ *
+ *
+ *
+
+ case synth_gaussianwhitenoise: {
+
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+ double U1, U2;
+
+
+ do {
+
+ U1 = DRANQD1;
+ U2 = DRANQD1;
+
+ S = U1 * U1 + U2 * U2;
+
+ } while (S >= 1 || S == 0);
+
+ V1 = U1 * sqrt((-2.0 * log(S)) / S);
+ V2 = U2 * sqrt((-2.0 * log(S)) / S);
+
+ if (phase == 0) {
+ synth_out = exp(V1);
+
+ } else {
+
+ synth_out = exp(V2);
+ }
+
+ phase = 1 - phase;
+
+ printf("%lf\n", synth_out);
+
+ break;
+
+ }
+*/
+
+
+
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #7
+ *
+ *
+ *
+ * seems close?
+
+
+ case synth_gaussianwhitenoise: {
+ static double U1, U2, S;
+ static int phase = 0;
+ double X;
+
+ if(phase == 0) {
+ do {
+
+ U1 = DRANQD1;
+ U2 = DRANQD1;
+
+ //V1 = 2 * U1 - 1;
+ //V2 = 2 * U2 - 1;
+ S = U1 * U1 + U2 * U2;
+
+ } while (S >= 1 || S == 0);
+
+ X = U1 * sqrt(-2 * log(S) / S);
+ synth_out = exp(X);
+ printf("%lf\n", synth_out);
+
+ } else {
+
+ X = U2 * sqrt(-2 * log(S) / S);
+ synth_out = exp(X);
+ printf("%lf\n", synth_out);
+ }
+
+ phase = 1 - phase;
+
+ }
+
+*/
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #8
+ *
+ *
+ *
+
+
+case synth_gaussianwhitenoise: {
+ double x1, x2, w, y1;
+ static double y2;
+ static int use_last = 0;
+ static double mean = 0;
+ static double stdev = 4;
+
+ if (use_last) /* use value from previous call
+ {
+ y1 = y2;
+ use_last = 0;
+ } else {
+ do {
+ x1 = DRANQD1;
+ x2 = DRANQD1;
+ w = x1 * x1 + x2 * x2;
+ } while ( w >= 1.0 || w == 0);
+
+ w = sqrt( (-2.0 * log( w ) ) / w );
+ y1 = x1 * w;
+ y2 = x2 * w;
+ use_last = 1;
+ }
+
+ synth_out = mean + y1 * stdev;
+ printf("%lf\n", synth_out);
+
+}
+
+*/
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #9
+ *
+ *
+ *
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+ static int mean = 0;
+ static int stddev = 1;
+
+ if(phase == 0) {
+ do {
+ double U1 = (double)rand() / RAND_MAX;
+ double U2 = (double)rand() / RAND_MAX;
+
+ V1 = 2 * U1 - 1;
+ V2 = 2 * U2 - 1;
+ S = V1 * V1 + V2 * V2;
+
+ } while (S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+
+ } else {
+
+ X = V2 * sqrt(-2 * log(S) / S);
+
+ }
+
+ phase = 1 - phase;
+
+ synth_out = X; //i.e. mean = 0 and stddev = 1
+
+ printf("%lf %lf\n", synth_out, exp(mean + stddev * X));
+ }
+*/
+
+
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #10
+ *
+ *
+ *
+
+// this one works? kind of?
+// figured out where should the phase switch should be?
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+ double X = 1;
+
+ do {
+ if(phase == 0) {
+ do {
+ double U1 = (double)rand() / RAND_MAX;
+ double U2 = (double)rand() / RAND_MAX;
+
+ V1 = 2 * U1 - 1;
+ V2 = 2 * U2 - 1;
+ S = V1 * V1 + V2 * V2;
+
+ } while (S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+
+ } else {
+
+ X = V2 * sqrt(-2 * log(S) / S);
+
+ }
+
+ phase = 1 - phase;
+
+ } while (X < -1 || X > 1);
+
+ //phase = 1 - phase; //loops infinitely
+
+ synth_out = X;
+ printf("%lf\n", synth_out);
+
+ }
+*/
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #11
+ *
+ *
+ *
+
+//this should be the best one too, reuses existing code, and has options
+//but the exp output gives numbers that are not anything like the X numbers: 0.367881 to 2.71828
+//i am officially confused
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+ static int mean = 0;
+ static int stddev = 1;
+
+ do {
+ if(phase == 0) {
+ do {
+ V1 = DRANQD1;
+ V2 = DRANQD1;
+
+ S = V1 * V1 + V2 * V2;
+
+ } while (S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+
+ } else {
+
+ X = V2 * sqrt(-2 * log(S) / S);
+
+ }
+
+ phase = 1 - phase;
+
+ } while (X < -1 || X > 1);
+
+ //phase = 1 - phase; // causes infinite loop
+
+ synth_out = X;
+
+ printf("%lf %lf\n", synth_out, exp(mean + stddev * X));
+
+ }
+
+*/
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #12
+ *
+ *
+ *
+
+// https://stackoverflow.com/questions/2325472/generate-random-numbers-following-a-normal-distribution-in-c-c
+ case synth_gaussianwhitenoise: {
+ static double y1, y2;
+ static double sigma = 82.;
+ static double Mi = 40.;
+
+ y1 = (DRANQD1 + 1.) / ((double)(RAND_MAX) + 1.);
+ y2 = (DRANQD1 + 1.) / ((double)(RAND_MAX) + 1.);
+ synth_out = (cos(2*3.14*y2)*sqrt(-2.*log(y1))) * sigma + Mi;
+
+ break;
+
+ }
+
+
+*/
+
+
+
+
+
+
+
+/*
+ *
+ *
+ * ATTEMPT #13
+ *
+ *
+ *
+ */
+
+//this should be the best one, reuses existing code
+
+ case synth_gaussianwhitenoise: {
+ static double V1, V2, S;
+ static int phase = 0;
+ double X;
+
+ do {
+ if(phase == 0) {
+ do {
+ V1 = DRANQD1;
+ V2 = DRANQD1;
+
+ S = V1 * V1 + V2 * V2;
+
+ } while (S >= 1 || S == 0);
+
+ X = V1 * sqrt(-2 * log(S) / S);
+
+ } else {
+
+ X = V2 * sqrt(-2 * log(S) / S);
+
+ }
+
+ phase = 1 - phase;
+
+ } while (X < -1 || X > 1);
+
+ synth_out = X;
+ printf("%lf\n", synth_out);
+
+ }
+
+
+
+
+
+
+
+
+
+
+
+
+
+ case synth_whitenoise: {
synth_out = DRANQD1;
+ printf("%lf\n", synth_out);
break;
+ }
+
+
case synth_tpdfnoise:
synth_out = .5 * (DRANQD1 + DRANQD1);
@@ -582,15 +1197,15 @@
case synth_pinknoise: { /* "Paul Kellet's refined method" */
#define _ .125 / (65536. * 32768.)
double d = RANQD1;
- chan->c0 = .99886 * chan->c0 + d * (.0555179*_);
- chan->c1 = .99332 * chan->c1 + d * (.0750759*_);
- chan->c2 = .96900 * chan->c2 + d * (.1538520*_);
- chan->c3 = .86650 * chan->c3 + d * (.3104856*_);
- chan->c4 = .55000 * chan->c4 + d * (.5329522*_);
- chan->c5 = -.7616 * chan->c5 - d * (.0168980*_);
+ chan->c0 = .99886 * chan->c0 + d * (.0555179*_);
+ chan->c1 = .99332 * chan->c1 + d * (.0750759*_);
+ chan->c2 = .96900 * chan->c2 + d * (.1538520*_);
+ chan->c3 = .86650 * chan->c3 + d * (.3104856*_);
+ chan->c4 = .55000 * chan->c4 + d * (.5329522*_);
+ chan->c5 = -.7616 * chan->c5 - d * (.0168980*_);
synth_out = chan->c0 + chan->c1 + chan->c2 + chan->c3
- + chan->c4 + chan->c5 + chan->c6 + d * (.5362*_);
- chan->c6 = d * (.115926*_);
+ + chan->c4 + chan->c5 + chan->c6 + d * (.5362*_);
+ chan->c6 = d * (.115926*_);
break;
#undef _
}
@@ -604,10 +1219,10 @@
case synth_pluck: {
double d = chan->buffer[chan->pos];
- chan->hp_last_out =
+ chan->hp_last_out =
(d - chan->hp_last_in) * chan->c3 + chan->hp_last_out * chan->c2;
chan->hp_last_in = d;
-
+
synth_out = range_limit(chan->hp_last_out, -1, 1);
chan->lp_last_out = d = d * chan->c1 + chan->lp_last_out * chan->c0;
[-- Attachment #3: Type: text/plain, Size: 0 bytes --]
[-- Attachment #4: Type: text/plain, Size: 158 bytes --]
_______________________________________________
SoX-devel mailing list
SoX-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/sox-devel
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2018-10-08 22:00 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2018-10-08 21:59 Updating synth.c with a Gaussian white noise effect? James Trammell
Code repositories for project(s) associated with this public inbox
https://80x24.org/mirrors/sox.git
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).