11 #include <sys/types.h>
13 #define MIN(a,b) ((a)<(b) ? (a) : (b))
14 #define MAX(a,b) ((a)>(b) ? (a) : (b))
16 void doFourier(double *input, int channels, sf_count_t len, fftw_complex *output)
18 static sf_count_t len0;
20 static double *windowedData;
21 static fftw_plan planlos;
22 static int planned = 0;
27 windowedData = fftw_malloc(sizeof(double) * channels * len);
30 fprintf(stderr, "Planning...\n");
31 planlos = fftw_plan_many_dft_r2c(1, &l, channels,
32 windowedData, NULL, channels, 1,
33 output, NULL, channels, 1,
34 FFTW_MEASURE | FFTW_PRESERVE_INPUT | FFTW_UNALIGNED);
38 assert(channels0 == channels);
40 sf_count_t nmax = channels * len;
42 for(i = 0; i < nmax; ++i)
44 double ang = (i/channels) * 2 * M_PI / (len - 1);
45 windowedData[i] = input[i] * (
52 fftw_execute_dft_r2c(planlos, windowedData, output);
55 double vectorDot(fftw_complex *v1, fftw_complex *v2, size_t length)
59 for(i = 0; i < length; ++i)
60 sum += cabs(v1[i]) * cabs(v2[i]);
64 sf_count_t findMaximumSingle(double (*func) (sf_count_t), sf_count_t x0, sf_count_t x1, sf_count_t step)
66 sf_count_t bestpos = x1;
67 double best = func(x1);
71 for(i = x0 + step; i < x1; i += step)
84 sf_count_t findMaximum(double (*func) (sf_count_t), sf_count_t x0, sf_count_t xg, sf_count_t xg2, sf_count_t x1)
88 xg0 = xg = MAX(x0, MIN(xg, x1));
89 xg20 = xg2 = MAX(x0, MIN(xg2, x1));
93 sf_count_t size = xg2 - xg;
96 fprintf(stderr, "round:\n");
97 sf_count_t bestguess = findMaximumSingle(func, xg, xg2, (xg2 - xg) / 16 + 1);
98 xg = MAX(x0, bestguess - size / 3);
99 xg2 = MIN(bestguess + size / 3, x1);
102 if(xg - xg0 < (xg20 - xg0) / 16)
103 fprintf(stderr, "warning: best match very close to left margin, maybe decrease the guess?\n");
105 if(xg20 - xg < (xg20 - xg0) / 16)
106 fprintf(stderr, "warning: best match very close to right margin, maybe increase the guess?\n");
111 int main(int argc, char **argv)
115 memset(&infile_info, 0, sizeof(infile_info));
116 infile_info.format = 0;
117 SNDFILE *infile = sf_open(argv[1], SFM_READ, &infile_info);
120 fprintf(stderr, "Input file has %ld frames, %d Hz, %d channels, format %08x, %d sections and is %s\n",
121 (long) infile_info.frames,
122 infile_info.samplerate,
123 infile_info.channels,
125 infile_info.sections,
126 infile_info.seekable ? "seekable" : "not seekable");
127 chans = infile_info.channels;
129 sf_count_t size = MIN(infile_info.frames, strtod(argv[3], NULL) * infile_info.samplerate);
130 sf_count_t guess = strtod(argv[4], NULL) * infile_info.samplerate;
131 sf_count_t guess2 = strtod(argv[5], NULL) * infile_info.samplerate;
132 int channels = infile_info.channels;
133 size_t fftsize = atoi(argv[2]);
134 size_t ndata = channels * (fftsize/2 + 1);
137 size_t ndata_lowpass = channels * (ndata / channels / 512); // 43 Hz
138 size_t ndata_highpass = channels * (ndata / channels / 4); // 5512 Hz
139 // simplifies the dot product and makes the target function more smooth
141 size_t ndata_lowpass = 0;
142 //size_t ndata_highpass = channels * (ndata / channels / 4); // 5512 Hz
143 size_t ndata_highpass = ndata;
145 if(size < guess + 2 * (sf_count_t) fftsize)
146 err(1, "sound file too small (maybe reduce fftsize?)");
149 err(1, "guess too close to file start (maybe reduce fftsize?)");
151 fprintf(stderr, "Using %ld frames of %d channels, guess at %ld, %d FFT size -> %d FFT result size\n", (long) size, channels, (long) guess, (int) fftsize, (int) ndata);
153 fftw_complex *data_end = fftw_malloc(sizeof(fftw_complex) * ndata);
154 fftw_complex *data_cur = fftw_malloc(sizeof(fftw_complex) * ndata);
156 double *sbuf = malloc(sizeof(double) * (size * channels));
157 if(size != sf_readf_double(infile, sbuf, size))
158 errx(1, "sf_read_double");
161 doFourier(sbuf + (size - fftsize) * chans, chans, fftsize, data_end);
162 fprintf(stderr, "end transformed.\n");
164 double sxx = vectorDot(data_end + ndata_lowpass, data_end + ndata_lowpass, ndata_highpass - ndata_lowpass);
166 errx(1, "ends with silence... use another end point");
172 double diffsum[ndata];
173 fftw_complex save[ndata];
174 for(i = guess; i < size - fftsize; ++i)
176 doFourier(sbuf + i * channels, channels, fftsize, data_cur);
177 for(j = 0; j < ndata; ++j)
179 fftw_complex x = data_cur[j];
182 fftw_complex y = save[j];
184 diffsum[j] += cabs(y - x);
188 fprintf(stderr, "at position %d\n", (int) i);
190 for(j = 0; j < ndata; ++j)
191 printf("%d %.9f %.9f\n", j, sum[j], diffsum[j]);
196 double similarityAt(sf_count_t i)
198 // A trampoline! Wheeeeeeeeeew!
199 doFourier(sbuf + i * channels, channels, fftsize, data_cur);
200 double sxy = vectorDot(data_end + ndata_lowpass, data_cur + ndata_lowpass, ndata_highpass - ndata_lowpass);
201 double syy = vectorDot(data_cur + ndata_lowpass, data_cur + ndata_lowpass, ndata_highpass - ndata_lowpass);
202 double v = syy ? ((sxy*sxy) / (sxx*syy)) : -1;
203 fprintf(stderr, "Evaluated at %.9f: %f\n", i / (double) infile_info.samplerate, v);
210 for(i = guess; i < size - 2 * fftsize; ++i)
212 printf("%.9f %f\n", i / (double) infile_info.samplerate, similarityAt(i));
218 sf_count_t best = findMaximum(similarityAt, 0, guess - fftsize, guess2 - fftsize, size - 2 * fftsize);
219 fprintf(stderr, "Result: %.9f (sample %ld)\n", (best + fftsize) / (double) infile_info.samplerate, (long) (best + fftsize));
223 // 1. Crossfading end to our start sample
224 fprintf(stderr, "Crossfading...\n");
228 for(j = 0; j < channels; ++j)
229 for(i = 0; i < fftsize; ++i)
231 double f1 = pow(fftsize - i, p);
232 double f2 = pow(i, p);
233 sbuf[(size - fftsize + i) * channels + j] =
235 sbuf[(size - fftsize + i) * channels + j] * f1
237 sbuf[(best + i) * channels + j] * f2
241 // 2. Open sound file
242 fprintf(stderr, "Opening...\n");
243 SNDFILE *outfile = sf_open(argv[6], SFM_WRITE, &infile_info);
247 fprintf(stderr, "Writing...\n");
248 // 1. Write samples from start to just before start of our end FFT
249 sf_writef_double(outfile, sbuf, size);
251 fprintf(stderr, "Closing...\n");
254 printf("%ld %.9f\n", (long) (best + fftsize), (best + fftsize) / (double) infile_info.samplerate);