1 /*
2   timpulse.c
3   David Rowe Dec 2019
4 
5   Generate a synthetic speech signal from a sum of sinusoids.  Generates a known test
6   signals for phaseNN and ampNN projects.
7 */
8 
9 #include <assert.h>
10 #include <math.h>
11 #include <stdlib.h>
12 #include <stdio.h>
13 #include <getopt.h>
14 
15 #define FS 8000
16 
main(int argc,char * argv[])17 int main(int argc, char *argv[]) {
18     short buf[FS] = {0};
19     float f0 = 60.0;
20     float n0 = 0.0;
21     int   Nsecs = 1;
22     int   randf0 = 0;
23     int   filter = 0;
24     int   rande = 0;
25 
26     int o = 0;
27     int opt_idx = 0;
28     while( o != -1 ) {
29         static struct option long_opts[] = {
30             {"help",   no_argument,       0, 'h'},
31             {"n0",     required_argument, 0, 'n'},
32             {"f0",     required_argument, 0, 'f'},
33             {"secs",   required_argument, 0, 's'},
34             {"randf0", no_argument, 0, 'r'},
35             {"rande",  required_argument, 0, 'e'},
36             {"filter", no_argument, 0, 'i'},
37             {0, 0, 0, 0}
38         };
39 
40         o = getopt_long(argc,argv,"hn:f:s:r",long_opts,&opt_idx);
41 
42         switch(o) {
43         case 'n':
44 	    n0 = atof(optarg);
45             break;
46         case 'f':
47             f0 = atof(optarg);
48 	    break;
49         case 's':
50             Nsecs = atoi(optarg);
51 	    break;
52         case 'r':
53             randf0 = 1;
54 	    break;
55         case 'i':
56             filter = 1;
57 	    break;
58         case 'e':
59             rande = atoi(optarg);
60 	    break;
61         case '?':
62         case 'h':
63 	    fprintf(stderr,
64 		    "usage: %s\n"
65 		    "[--f0 f0Hz]          fixed F0\n"
66                     "[--n0 samples]       time offset\n"
67                     "[--secs Nsecs]       number of seconds to generate\n"
68 	            "[--randf0]           choose a random F0 every second\n"
69 	            "[--rande Ndiscrete]  choose a random frame energy every second, Ndiscrete values\n"
70 		    "\n", argv[0]);
71 	    exit(1);
72 	break;
73         }
74     }
75 
76     int t = 0;
77     float A = 100.0;
78 
79     /* optionally filter with 2nd order system */
80     float alpha = 0.25*M_PI, gamma=0.99;
81     float a[2] = {-2.0*gamma*cos(alpha), gamma*gamma};
82     float mem[2] = {0};
83 
84     for (int j=0; j<Nsecs; j++) {
85 	if (rande) {
86 	    float AdB_min = 20.0*log10(100.0);
87 	    float AdB_step = 6.0;
88 	    float num_values = rande;
89 
90 	    // discrete RV between 0..1
91 	    float r = (float)rand()/RAND_MAX;
92 	    r = floor(r*num_values);
93 
94 	    float AdB = AdB_min + r*AdB_step;
95 	    A = pow(10.0,AdB/20.0);
96 	    fprintf(stderr, "r: %f AdB: %f A: %f\n", r, AdB, A);
97 	}
98 	if (randf0) {
99 	    float pitch_period = FS/400.0 + (FS/80.0 - FS/400.0)*rand()/RAND_MAX;
100 	    f0 = (float)FS/pitch_period;
101 	    //fprintf(stderr, "P: %f f0: %f\n", pitch_period, f0);
102 	}
103 	float Wo = 2.0*M_PI*f0/FS;
104 	int L = M_PI/Wo;
105 	float e = 0.0;
106 	for(int i=0; i<FS; i++) {
107 	    buf[i] = 0;
108 	    // 1/sqrt(L) term makes power constant across Wo
109 	    for(int m=1; m<L; m++)
110 		buf[i] += (A/sqrt(L))*cos(m*Wo*(t + n0));
111 	    e += pow(buf[i], 2.0);
112 	    t++;
113 	}
114 	//fprintf(stderr, "e (dB): %f\n", 10*log10(e));
115 	if (filter) {
116 	    for(int i=0; i<FS; i++) {
117 		float x = (float)buf[i];
118 		float y = (x - mem[0]*a[0] - mem[1]*a[1]);
119 		mem[1] = mem[0]; mem[0] = y;
120 		buf[i] = (short)y;
121 	    }
122 	}
123 
124 	fwrite(buf, sizeof(short), FS, stdout);
125     }
126 
127     return 0;
128 }
129