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; 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));
91 fprintf(stderr, "min/max: %d %d\n", (int)xg, (int)xg2);
95 sf_count_t size = xg2 - xg;
98 //fprintf(stderr, "round:\n");
99 sf_count_t bestguess = findMaximumSingle(func, xg, xg2, size / 32 + 1);
100 xg = MAX(xg, bestguess - size / 3);
101 xg2 = MIN(bestguess + size / 3, xg2);
104 fprintf(stderr, "guessed: %d\n", (int)xg);
106 if(xg - xg0 < (xg20 - xg0) / 64)
107 fprintf(stderr, "warning: best match very close to left margin, maybe decrease the guess?\n");
109 if(xg20 - xg < (xg20 - xg0) / 64)
110 fprintf(stderr, "warning: best match very close to right margin, maybe increase the guess?\n");
115 int main(int argc, char **argv)
119 memset(&infile_info, 0, sizeof(infile_info));
120 infile_info.format = 0;
121 SNDFILE *infile = sf_open(argv[1], SFM_READ, &infile_info);
124 fprintf(stderr, "Input file has %ld frames, %d Hz, %d channels, format %08x, %d sections and is %s\n",
125 (long) infile_info.frames,
126 infile_info.samplerate,
127 infile_info.channels,
129 infile_info.sections,
130 infile_info.seekable ? "seekable" : "not seekable");
131 chans = infile_info.channels;
133 sf_count_t size = MIN(infile_info.frames, strtod(argv[3], NULL) * infile_info.samplerate);
134 sf_count_t guess = strtod(argv[4], NULL) * infile_info.samplerate;
135 sf_count_t guess2 = strtod(argv[5], NULL) * infile_info.samplerate;
136 int channels = infile_info.channels;
137 size_t fftsize = atoi(argv[2]);
138 size_t ndata = channels * (fftsize/2 + 1);
141 size_t ndata_lowpass = channels * (ndata / channels / 512); // 43 Hz
142 size_t ndata_highpass = channels * (ndata / channels / 4); // 5512 Hz
143 // simplifies the dot product and makes the target function more smooth
145 size_t ndata_lowpass = 0;
146 //size_t ndata_highpass = channels * (ndata / channels / 4); // 5512 Hz
147 size_t ndata_highpass = ndata;
149 if(size < guess + 2 * (sf_count_t) fftsize)
150 err(1, "sound file too small (maybe reduce fftsize?)");
153 err(1, "guess too close to file start (maybe reduce fftsize?)");
155 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);
157 fftw_complex *data_end = fftw_malloc(sizeof(fftw_complex) * ndata);
158 fftw_complex *data_cur = fftw_malloc(sizeof(fftw_complex) * ndata);
160 double *sbuf = malloc(sizeof(double) * (size * channels));
161 if(size != sf_readf_double(infile, sbuf, size))
162 errx(1, "sf_read_double");
165 doFourier(sbuf + (size - fftsize) * chans, chans, fftsize, data_end);
166 fprintf(stderr, "end transformed.\n");
168 double sxx = vectorDot(data_end + ndata_lowpass, data_end + ndata_lowpass, ndata_highpass - ndata_lowpass);
170 errx(1, "ends with silence... use another end point");
176 double diffsum[ndata];
177 fftw_complex save[ndata];
178 for(i = guess; i < size - fftsize; ++i)
180 doFourier(sbuf + i * channels, channels, fftsize, data_cur);
181 for(j = 0; j < ndata; ++j)
183 fftw_complex x = data_cur[j];
186 fftw_complex y = save[j];
188 diffsum[j] += cabs(y - x);
192 fprintf(stderr, "at position %d\n", (int) i);
194 for(j = 0; j < ndata; ++j)
195 printf("%d %.9f %.9f\n", j, sum[j], diffsum[j]);
200 double similarityAt(sf_count_t i)
202 // A trampoline! Wheeeeeeeeeew!
203 doFourier(sbuf + i * channels, channels, fftsize, data_cur);
204 double sxy = vectorDot(data_end + ndata_lowpass, data_cur + ndata_lowpass, ndata_highpass - ndata_lowpass);
205 double syy = vectorDot(data_cur + ndata_lowpass, data_cur + ndata_lowpass, ndata_highpass - ndata_lowpass);
206 double v = syy ? ((sxy*sxy) / (sxx*syy)) : -1;
207 //fprintf(stderr, "Evaluated at %.9f: %f\n", i / (double) infile_info.samplerate, v);
214 for(i = guess; i < size - 2 * fftsize; ++i)
216 printf("%.9f %f\n", i / (double) infile_info.samplerate, similarityAt(i));
222 sf_count_t best = findMaximum(similarityAt, 0, guess - fftsize, guess2 - fftsize, size - 2 * fftsize);
223 fprintf(stderr, "Result: %.9f (sample %ld)\n", (best + fftsize) / (double) infile_info.samplerate, (long) (best + fftsize));
227 // 1. Crossfading end to our start sample
228 fprintf(stderr, "Crossfading...\n");
232 for(j = 0; j < channels; ++j)
233 for(i = 0; i < fftsize; ++i)
235 double f1 = pow(fftsize - i, p);
236 double f2 = pow(i, p);
237 sbuf[(size - fftsize + i) * channels + j] =
239 sbuf[(size - fftsize + i) * channels + j] * f1
241 sbuf[(best + i) * channels + j] * f2
245 // 2. Open sound file
246 fprintf(stderr, "Opening...\n");
247 SNDFILE *outfile = sf_open(argv[6], SFM_WRITE, &infile_info);
251 fprintf(stderr, "Writing...\n");
252 // 1. Write samples from start to just before start of our end FFT
253 sf_writef_double(outfile, sbuf, size);
255 fprintf(stderr, "Closing...\n");
258 printf("%ld %.9f\n", (long) (best + fftsize), (best + fftsize) / (double) infile_info.samplerate);