1 /*
2     Example implementation of the PADsynth basic algorithm
3     By: Nasca O. Paul, Tg. Mures, Romania
4 
5     Ported to pure C by Paul Batchelor
6 
7     This implementation and the algorithm are released under Public Domain
8     Feel free to use it into your projects or your products ;-)
9 
10     This implementation is tested under GCC/Linux, but it's
11     very easy to port to other compiler/OS.
12 */
13 
14 #include <stdlib.h>
15 #include <math.h>
16 #include "soundpipe.h"
sp_osc_create(sp_osc ** osc)17 
18 #ifndef M_PI
19 #define M_PI		3.14159265358979323846
20 #endif
21 
22 int sp_gen_padsynth(sp_data *sp, sp_ftbl *ps, sp_ftbl *amps,
23         SPFLOAT f, SPFLOAT bw)
24 {
25 
26     int i, nh;
27     int N = (int) ps->size;
28     int number_harmonics = (int) amps->size;
29     SPFLOAT *A = amps->tbl;
30     SPFLOAT *smp = ps->tbl;
31 
32     SPFLOAT *freq_amp = malloc((N / 2) * sizeof(SPFLOAT));
33     SPFLOAT *freq_phase = malloc((N / 2) * sizeof(SPFLOAT));
34 
35     for (i=0;i<N/2;i++) freq_amp[i]=0.0;
36 
37     for (nh=1;nh<number_harmonics;nh++) {
38         SPFLOAT bw_Hz;
39         SPFLOAT bwi;
40         SPFLOAT fi;
41         bw_Hz = (pow(2.0, bw/1200.0) - 1.0) * f * nh;
42         bwi = bw_Hz/(2.0*ps->size);
43         fi = f*nh/ps->size;
44         for (i = 0; i < N/2 ; i++) {
45             SPFLOAT hprofile;
46             hprofile = sp_padsynth_profile((i / (SPFLOAT) N) - fi, bwi);
47             freq_amp[i] += hprofile*A[nh];
48         }
49     }
50 
51     for (i=0;i<N/2;i++) {
52         freq_phase[i]= (sp_rand(sp) / (SP_RANDMAX + 1.0)) * 2.0 * M_PI;
53     };
54 
55     sp_padsynth_ifft(N,freq_amp,freq_phase,smp);
56     sp_padsynth_normalize(N,smp);
57 
58     free(freq_amp);
59     free(freq_phase);
60     return SP_OK;
61 }
62 
63 /* This is the profile of one harmonic
64    In this case is a Gaussian distribution (e^(-x^2))
65    The amplitude is divided by the bandwidth to ensure that the harmonic
66    keeps the same amplitude regardless of the bandwidth */
67 
68 SPFLOAT sp_padsynth_profile(SPFLOAT fi, SPFLOAT bwi)
69 {
70     SPFLOAT x =fi/bwi;
71     x *= x;
72 
73 /*
74  * this avoids computing the e^(-x^2) where
75  * it's results are very close to zero
76  */
77     if (x>14.71280603) return 0.0;
78 
79     return exp(-x)/bwi;
80 }
81 
82 int sp_padsynth_ifft(int N, SPFLOAT *freq_amp,
83         SPFLOAT *freq_phase, SPFLOAT *smp)
84 {
85     int i;
86     FFTwrapper *fft;
87     FFTwrapper_create(&fft, N);
88     FFTFREQS fftfreqs;
89     newFFTFREQS(&fftfreqs,N/2);
90 
91     for (i=0; i<N/2; i++){
92         fftfreqs.c[i]=freq_amp[i]*cos(freq_phase[i]);
93         fftfreqs.s[i]=freq_amp[i]*sin(freq_phase[i]);
94     };
95     freqs2smps(fft, &fftfreqs,smp);
96     deleteFFTFREQS(&fftfreqs);
97     FFTwrapper_destroy(&fft);
98     return SP_OK;
99 }
100 
101 /*
102     Simple normalization function. It normalizes the sound to 1/sqrt(2)
103 */
104 
105 int sp_padsynth_normalize(int N, SPFLOAT *smp)
106 {
107     int i;
108     SPFLOAT max=0.0;
109     for (i=0;i<N;i++) if (fabs(smp[i])>max) max=fabs(smp[i]);
110     if (max<1e-5) max=1e-5;
111     for (i=0;i<N;i++) smp[i]/=max*1.4142;
112     return SP_OK;
113 }
114