1 module soundtab.audio.wavefile; 2 3 immutable int PERIOD_FFT_SIZE = 256; 4 immutable int LPC_SIZE = 20; 5 struct PeriodInfo { 6 float startTime = 0; 7 float periodTime = 0; 8 float ampPlus = 0; 9 float ampMinus = 0; 10 float energy = 0; 11 float[PERIOD_FFT_SIZE / 2] fftAmp; 12 float[PERIOD_FFT_SIZE / 2] fftPhase; 13 float[LPC_SIZE] lpc; 14 float[LPC_SIZE] lsp; 15 @property float middleTime() { 16 return startTime + periodTime / 2; 17 } 18 @property float endTime() { 19 return startTime + periodTime; 20 } 21 } 22 23 class WaveFile { 24 string filename; 25 int channels; 26 int sampleRate; 27 int frames; 28 float[] data; 29 30 float[] marks; 31 float[] negativeMarks; 32 33 PeriodInfo[] periods; 34 35 float middleFrequency = 0; 36 float[] frequencies; // multipliers relative to middle frequency, to get real frequency, use frequency[i]*middleFrequency 37 float minFrequencyK = 0; 38 float maxFrequencyK = 0; 39 float[LPC_SIZE] lpc; 40 float[LPC_SIZE] lsp; 41 float[] excitation; 42 43 this() { 44 } 45 46 this(WaveFile v, bool forceMono = false) { 47 channels = v.channels; 48 sampleRate = v.sampleRate; 49 frames = v.frames; 50 if (channels == 1 || !forceMono) { 51 // just duplicate 52 data = v.data.dup; 53 } else { 54 // take left channel only 55 data = new float[frames]; 56 for (int i = 0; i < frames; i++) 57 data[i] = v.data[i * channels]; 58 channels = 1; 59 } 60 if (v.marks.length) 61 marks = v.marks.dup; 62 middleFrequency = v.middleFrequency; 63 if (v.frequencies.length) 64 frequencies = v.frequencies.dup; // multipliers relative to middle frequency, to get real frequency, use frequency[i]*middleFrequency 65 minFrequencyK = v.minFrequencyK; 66 maxFrequencyK = v.maxFrequencyK; 67 } 68 69 void fill(float value) { 70 data[0..$] = value; 71 } 72 73 void setMarks(float[] _marks, float[] _negativeMarks = null) { 74 marks = _marks; 75 negativeMarks = _negativeMarks; 76 } 77 78 void fillPeriodsFromMarks() { 79 periods = null; 80 if (marks.length > 1) { 81 periods.length = marks.length - 1; 82 for (int i = 0; i + 1 < marks.length; i++) { 83 fillPeriodInfo(periods[i], marks[i], marks[i + 1]); 84 } 85 } 86 if (marks.length > 3) { 87 // LSPs on edges are calculated inaccurate - copy from inner frames 88 periods[0].lpc[0..$] = periods[1].lpc[0..$]; 89 periods[0].lsp[0..$] = periods[1].lsp[0..$]; 90 periods[$-1].lpc[0..$] = periods[$-2].lpc[0..$]; 91 periods[$-1].lsp[0..$] = periods[$-2].lsp[0..$]; 92 } 93 // calc avg lsp/lpc 94 int startIndex = cast(int)(periods.length / 3); 95 int endIndex = cast(int)(periods.length * 2 / 3); 96 lsp[0..$] = 0; 97 for (int i = startIndex; i < endIndex; i++) { 98 for (int j = 0; j < LPC_SIZE; j++) { 99 lsp[j] += periods[i].lsp[j]; 100 } 101 } 102 for (int j = 0; j < LPC_SIZE; j++) { 103 lsp[j] /= endIndex - startIndex; 104 } 105 lspToLpc(lsp, lpc); 106 calcExcitation(); 107 } 108 109 void calcExcitation() { 110 import dlangui.core.logger; 111 excitation.length = data.length; 112 // inverse filtering 113 float dataEnergy = 0; 114 float excitationEnergy = 0; 115 for (int i = 0; i < data.length; i++) { 116 float predicted = 0; 117 for (int j = 0; j < LPC_SIZE; j++) { 118 int index = i - j - 1; 119 float v = index >= 0 ? data[index] : 0; 120 predicted -= v * lpc[j]; 121 } 122 //predicted = -predicted; 123 excitation[i] = data[i] - predicted; 124 dataEnergy += data[i] * data[i]; 125 excitationEnergy += excitation[i] * excitation[i]; 126 } 127 dataEnergy /= frames; 128 excitationEnergy /= frames; 129 Log.d("Excitation: energy = ", excitationEnergy, ", data energy = ", dataEnergy); 130 float[] tmpReconstructed; 131 tmpReconstructed.length = excitation.length; 132 dataEnergy = 0; 133 for (int i = 0; i < data.length; i++) { 134 float predicted = 0; 135 for (int j = 0; j < LPC_SIZE; j++) { 136 int index = i - j - 1; 137 float v = index >= 0 ? tmpReconstructed[index] : 0; 138 predicted -= v * lpc[j]; 139 } 140 //predicted = -predicted; 141 tmpReconstructed[i] = predicted + excitation[i]; 142 dataEnergy += tmpReconstructed[i] * tmpReconstructed[i]; 143 } 144 dataEnergy /= frames; 145 Log.d("Reconstructed data energy = ", dataEnergy); 146 } 147 148 immutable bool sqEnergy = true; 149 float[] amplitudes; 150 /// returns sign of biggest amplitude 151 int fillAmplitudesFromPeriods() { 152 amplitudes = null; 153 if (periods.length > 1) { 154 float[] ampPlus = new float[frames]; 155 float[] ampMinus = new float[frames]; 156 float[] energy = new float[frames]; 157 int prevPeriodMiddleFrame = 0; 158 float prevPeriodAmpPlus = periods[0].ampPlus; 159 float prevPeriodAmpMinus = periods[0].ampMinus; 160 float prevPeriodEnergy = periods[0].energy; 161 float avgPeriod = 0; 162 for (int i = 0; i < periods.length; i++) { 163 avgPeriod += periods[i].periodTime; 164 int middleFrame = timeToFrame(periods[i].middleTime); 165 float currentAmpPlus = periods[i].ampPlus; 166 float currentAmpMinus = periods[i].ampMinus; 167 float currentEnergy = periods[i].energy; 168 int frameLen = middleFrame - prevPeriodMiddleFrame; 169 for (int j = 0; j < frameLen; j++) { 170 ampPlus[prevPeriodMiddleFrame + j] = prevPeriodAmpPlus + (currentAmpPlus - prevPeriodAmpPlus) * j / frameLen; 171 ampMinus[prevPeriodMiddleFrame + j] = prevPeriodAmpMinus + (currentAmpMinus - prevPeriodAmpMinus) * j / frameLen; 172 energy[prevPeriodMiddleFrame + j] = prevPeriodEnergy + (currentEnergy - prevPeriodEnergy) * j / frameLen; 173 } 174 prevPeriodMiddleFrame = middleFrame; 175 prevPeriodAmpPlus = currentAmpPlus; 176 prevPeriodAmpMinus = currentAmpMinus; 177 prevPeriodEnergy = currentEnergy; 178 } 179 avgPeriod /= periods.length; 180 // fill till end of wave 181 for (int i = prevPeriodMiddleFrame; i < frames; i++) { 182 ampPlus[i] = prevPeriodAmpPlus; 183 ampMinus[i] = prevPeriodAmpMinus; 184 energy[i] = prevPeriodEnergy; 185 } 186 int lowpassFilterSize = cast(int)(10 * (1/avgPeriod)) | 1; 187 float[] lowpassFirFilter = blackmanWindow(lowpassFilterSize); 188 float[] ampPlusFiltered = new float[frames]; 189 float[] ampMinusFiltered = new float[frames]; 190 float[] energyFiltered = new float[frames]; 191 applyFirFilter(ampPlus, ampPlusFiltered, lowpassFirFilter); 192 applyFirFilter(ampMinus, ampMinusFiltered, lowpassFirFilter); 193 applyFirFilter(energy, energyFiltered, lowpassFirFilter); 194 amplitudes = new float[frames]; 195 float maxAmpPlus = ampPlusFiltered[0]; 196 float maxAmpMinus = ampMinusFiltered[0]; 197 float maxEnergy = energyFiltered[0]; 198 for (int i = 1; i < frames; i++) { 199 if (maxAmpPlus < ampPlusFiltered[i]) 200 maxAmpPlus = ampPlusFiltered[i]; 201 if (maxAmpMinus < ampMinusFiltered[i]) 202 maxAmpMinus = ampMinusFiltered[i]; 203 if (maxEnergy < energyFiltered[i]) 204 maxEnergy = energyFiltered[i]; 205 } 206 float maxAmp = maxAmpPlus > maxAmpMinus ? maxAmpPlus : maxAmpMinus; 207 float kPlus = 0; 208 float kMinus = 0; 209 float kEnergy = 0; 210 //if (maxAmpPlus > maxAmpMinus) 211 // amplitudes = ampPlusFiltered; 212 //else 213 // amplitudes = ampMinusFiltered; 214 kPlus = maxAmpPlus / (maxAmpPlus + maxAmpMinus); 215 kMinus = maxAmpMinus / (maxAmpPlus + maxAmpMinus); 216 for (int i = 0; i < frames; i++) { 217 amplitudes[i] = kPlus * ampPlusFiltered[i] + kMinus * ampMinusFiltered[i] + kEnergy * maxAmp / maxEnergy; 218 //amplitudes[i] = energyFiltered[i] * maxAmp / maxEnergy; 219 } 220 return maxAmpPlus > maxAmpMinus ? 1 : -1; 221 } 222 return 1; 223 } 224 225 void normalizeAmplitude() { 226 if (amplitudes && amplitudes.length == frames) { 227 for (int i = 0; i < frames; i++) { 228 data.ptr[i] /= amplitudes.ptr[i]; 229 } 230 } 231 } 232 233 void correctMarksForNormalizedAmplitude() { 234 correctMarksForNormalizedAmplitude(marks, 1); 235 //smoothTimeMarks(marks); 236 //smoothTimeMarks(marks); 237 correctMarksForNormalizedAmplitude(negativeMarks, -1); 238 //smoothTimeMarks(negativeMarks); 239 //smoothTimeMarks(negativeMarks); 240 } 241 242 void correctMarksForNormalizedAmplitude(ref float[] marks, int sign) { 243 if (marks.length < 3) 244 return; 245 for (int i = 1; i + 1 < marks.length; i++) { 246 float period = (marks[i + 1] - marks[i - 1]) / 2; 247 float freq = 1 / period; 248 float zeroPhaseTimePosSeconds; 249 float amp; 250 findNearPhase0(marks[i], freq, sign, zeroPhaseTimePosSeconds, amp); 251 float diff = zeroPhaseTimePosSeconds - marks[i]; 252 if (diff < period / 5 && diff > -period / 5) 253 marks[i] = zeroPhaseTimePosSeconds; 254 } 255 for (int i = 0; i < 8; i++) 256 smoothTimeMarksShifted(marks, negativeMarks); 257 } 258 259 void fillPeriodInfo(ref PeriodInfo period, float start, float end) { 260 import std.math : sqrt; 261 period.startTime = start; 262 period.periodTime = end - start; 263 int startSample = timeToFrame(start); 264 int endSample = timeToFrame(end); 265 float maxValue = 0; 266 float minValue = 0; 267 float energy = 0; 268 float energyDiv = 0; 269 for (int i = startSample; i < endSample; i++) { 270 if (i >= 0 && i < frames) { 271 float sample = data.ptr[i * channels]; 272 float absSample = sample < 0 ? -sample : sample; 273 if (maxValue < sample) 274 maxValue = sample; 275 if (minValue > sample) 276 minValue = sample; 277 if (i == startSample) { 278 float part = start * sampleRate - startSample; 279 part = 1 - part; 280 if (sqEnergy) 281 energy += absSample * absSample * part * part; 282 else 283 energy += absSample * part; 284 energyDiv = energyDiv + part; 285 } else if (i == endSample) { 286 float part = end * sampleRate - endSample; 287 if (sqEnergy) 288 energy += absSample * absSample * part * part; 289 else 290 energy += absSample * part; 291 energyDiv = energyDiv + part; 292 } else { 293 if (sqEnergy) 294 energy += absSample * absSample; 295 else 296 energy += absSample; 297 energyDiv = energyDiv + 1; 298 } 299 } 300 } 301 energy = energy / energyDiv; 302 if (sqEnergy) 303 energy = sqrt(energy); 304 period.ampPlus = maxValue; 305 period.ampMinus = -minValue; 306 period.energy = energy; 307 double[] src = new double[PERIOD_FFT_SIZE]; 308 float step = (end - start) / PERIOD_FFT_SIZE; 309 float t = start; 310 for (int i = 0; i < PERIOD_FFT_SIZE; i++) { 311 src[i] = getSampleInterpolated(t); 312 t += step; 313 } 314 import std.numeric; 315 import std.math; 316 import std.complex; 317 Complex!double[] res = new Complex!double[PERIOD_FFT_SIZE]; 318 static Fft periodFFT = null; 319 if (!periodFFT) 320 periodFFT = new Fft(PERIOD_FFT_SIZE); 321 periodFFT.fft(src[0..$], res[0..$]); 322 for (int i = 0; i < PERIOD_FFT_SIZE / 2; i++) { 323 double re = res[i].re; 324 double im = res[i].im; 325 double amp = sqrt((re*re + im*im) / PERIOD_FFT_SIZE); 326 double phase = 0; 327 if (amp > 0.0001) { 328 phase = atan2(re, im); 329 } else { 330 amp = 0; 331 } 332 period.fftAmp[i] = amp; 333 period.fftPhase[i] = phase; 334 } 335 import dlangui.core.logger; 336 //Log.d("period[", period.startTime, "] fft amps=", period.fftAmp); 337 338 float periodCenter = period.middleTime * sampleRate; // center frame 339 float periodTime = period.periodTime * sampleRate; // period frames 340 int windowSize = cast(int)(5 * periodTime); 341 float[] samples = new float[windowSize]; 342 float[] window = blackmanWindow(windowSize); 343 getSamplesInterpolated(cast(int)periodCenter, 1, samples); 344 for(int i = 0; i < samples.length; i++) { 345 samples[i] *= window[i]; 346 } 347 float error = calcLPC(samples, period.lpc[0..$]); 348 lpcToLsp(period.lpc[0..$], period.lsp[0..$]); 349 lspToLpc(period.lsp[0..$], period.lpc[0..$]); 350 Log.d("LPC256: ", period.lpc, " error=", error, "\n lsp=", period.lsp[0..$]); 351 /* 352 getSamplesInterpolated(periodCenter, periodTime / 200, samples); 353 for(int i = 0; i < samples.length; i++) { 354 samples[i] *= window[i]; 355 } 356 error = calcLPC(samples, period.lpc[0..$]); 357 lpcToLsp(period.lpc[0..$], period.lsp[0..$]); 358 lspToLpc(period.lsp[0..$], period.lpc[0..$]); 359 Log.d(" 200: ", period.lpc, " error=", error, "\n lsp=", period.lsp[0..$]); 360 361 getSamplesInterpolated(periodCenter, periodTime / 128, samples); 362 for(int i = 0; i < samples.length; i++) { 363 samples[i] *= window[i]; 364 } 365 error = calcLPC(samples, period.lpc[0..$]); 366 lpcToLsp(period.lpc[0..$], period.lsp[0..$]); 367 lspToLpc(period.lsp[0..$], period.lpc[0..$]); 368 Log.d(" 128: ", period.lpc, " error=", error, "\n lsp=", period.lsp[0..$]); 369 370 getSamplesInterpolated(periodCenter, periodTime / 64, samples); 371 for(int i = 0; i < samples.length; i++) { 372 samples[i] *= window[i]; 373 } 374 error = calcLPC(samples, period.lpc[0..$]); 375 lpcToLsp(period.lpc[0..$], period.lsp[0..$]); 376 lspToLpc(period.lsp[0..$], period.lpc[0..$]); 377 Log.d(" 064: ", period.lpc, " error=", error, "\n lsp=", period.lsp[0..$]); 378 379 */ 380 } 381 382 void smoothLSP() { 383 float[LPC_SIZE][] lsps; 384 lsps.length = periods.length; 385 for(int i = 0; i < periods.length; i++) { 386 lsps[i][0..$] = periods[i].lsp[0..$]; 387 } 388 for(int i = 0; i < periods.length; i++) { 389 for (int j = 0; j < LPC_SIZE; j++) 390 periods[i].lsp[j] = (lsps[i][j] + lsps[i > 0 ? i - 1 : i][j] + lsps[i + 1 < periods.length ? i + 1 : i][j]) / 3; 391 lspToLpc(periods[i].lsp[0..$], periods[i].lpc[0..$]); 392 } 393 } 394 395 void smoothMarks() { 396 smoothTimeMarks(marks); 397 } 398 399 void generateFrequenciesFromMarks() { 400 frequencies = null; 401 middleFrequency = 0; 402 minFrequencyK = 0; 403 maxFrequencyK = 0; 404 if (!data.length || marks.length < 3) 405 return; 406 frequencies.length = data.length; 407 double s = 0; 408 // calc avg frequency 409 float lastFreq = 1 / (marks[1] - marks[0]); 410 double minFreq = lastFreq; 411 double maxFreq = lastFreq; 412 int lastFreqPos = 0; 413 for (int i = 1; i < marks.length; i++) { 414 float freq = 1 / (marks[i] - marks[i - 1]); 415 if (minFreq > freq) 416 minFreq = freq; 417 if (maxFreq < freq) 418 maxFreq = freq; 419 s += freq; 420 int pos = timeToFrame((marks[i] + marks[i - 1]) / 2); 421 for (int j = lastFreqPos; j < pos; j++) { 422 float f = lastFreq + (freq - lastFreq) * (j - lastFreqPos) / (pos - lastFreqPos); 423 frequencies[j] = f; 424 } 425 lastFreq = freq; 426 lastFreqPos = pos; 427 } 428 for (int j = lastFreqPos; j < frequencies.length; j++) { 429 frequencies[j] = lastFreq; 430 } 431 middleFrequency = s / (marks.length - 1); 432 433 434 435 //int lowpassFilterSize = 3 * timeToFrame(1/middleFrequency) | 1; 436 //float[] lowpassFirFilter = blackmanWindow(lowpassFilterSize); 437 //float[] frequenciesFiltered = new float[frames]; 438 //applyFirFilter(frequencies, frequenciesFiltered, lowpassFirFilter); 439 //frequencies = frequenciesFiltered; 440 441 442 443 /* 444 445 for (int i = 1; i < marks.length; i++) { 446 double period = marks[i] - marks[i - 1]; 447 double freq = 1 / period; 448 if (i == 1 || minFreq > freq) 449 minFreq = freq; 450 if (i == 1 || maxFreq < freq) 451 maxFreq = freq; 452 } 453 middleFrequency = (maxFreq + minFreq) / 2; 454 minFrequencyK = minFreq / middleFrequency; 455 maxFrequencyK = maxFreq / middleFrequency; 456 int firstFilled = 0; 457 int lastFilled = 0; 458 for (int i = 1; i + 1 < marks.length; i++) { 459 float prevFreq = 1 / (marks[i] - marks[i - 1]); 460 float nextFreq = 1 / (marks[i + 1] - marks[i]); 461 int prevPeriodMiddle = timeToFrame((marks[i] + marks[i - 1]) / 2); 462 int nextPeriodMiddle = timeToFrame((marks[i + 1] + marks[i]) / 2); 463 float prevFreqK = prevFreq / middleFrequency; 464 float nextFreqK = nextFreq / middleFrequency; 465 // interpolate coeffs 466 for (int j = prevPeriodMiddle; j < nextPeriodMiddle; j++) { 467 frequencies[j] = prevFreqK + (nextFreqK - prevFreqK) * (j - prevPeriodMiddle) / (nextPeriodMiddle / prevPeriodMiddle); 468 } 469 if (!firstFilled) 470 firstFilled = prevPeriodMiddle; 471 lastFilled = nextPeriodMiddle - 1; 472 } 473 frequencies[0 .. firstFilled] = frequencies[firstFilled]; 474 frequencies[lastFilled .. $] = frequencies[lastFilled - 1]; 475 */ 476 } 477 478 void removeDcOffset(float minTime, float maxTime) { 479 int startFrame = timeToFrame(minTime); 480 int endFrame = timeToFrame(maxTime); 481 double s = 0; 482 for (int i = startFrame; i < endFrame; i++) 483 s += data.ptr[i]; 484 s /= (endFrame - startFrame); 485 for (int i = 0; i < data.length; i++) 486 data.ptr[i] -= s; 487 } 488 489 int timeToFrame(float time) { 490 return cast(int)(time * sampleRate); 491 } 492 float frameToTime(int frame) { 493 return (cast(float)frame / sampleRate); 494 } 495 float frameToTime(float frame) { 496 return (frame / sampleRate); 497 } 498 void limitFrameIndex(ref int index) { 499 if (index >= frames) 500 index = frames; 501 if (index < 0) 502 index = 0; 503 } 504 float getSample(int index, int channel = 0) { 505 if (index < 0 || index >= frames) 506 return 0; 507 return data.ptr[index * channels + channel]; 508 } 509 void getSamples(int startIndex, int channel, float[] buf) { 510 int len = cast(int)buf.length; 511 for (int i = 0; i < len; i++) { 512 int index = startIndex + i; 513 buf.ptr[i] = index >= 0 && index < frames ? data.ptr[index] : 0.0f; 514 } 515 } 516 void getSamplesZpad4(int startIndex, int channel, float[] buf) { 517 int len = cast(int)buf.length; 518 for (int i = 0; i * 4 < len; i++) { 519 int index = startIndex + i; 520 float sample = getSample(index); 521 buf.ptr[i * 4] = sample; 522 buf.ptr[i * 4 + 1] = 0; 523 buf.ptr[i * 4 + 2] = 0; 524 buf.ptr[i * 4 + 3] = 0; 525 } 526 } 527 /// linearly interpolated sample by time 528 float getSampleInterpolated(float time, int channel = 0) { 529 if (channel >= channels) 530 channel = channel % channels; 531 int index = cast(int)(time * sampleRate); 532 float deltaTime = time - cast(float)index / sampleRate; 533 float s0 = getSample(index, channel); 534 float s1 = getSample(index + 1, channel); 535 return s0 * (1 - deltaTime) + s1 * deltaTime; 536 } 537 538 void getSamplesInterpolated(float frameMiddle, float step, float[] buf) { 539 int len = cast(int)buf.length; 540 float x = frameMiddle - (len / 2) * step; 541 for (int i = 0; i < len; i++) { 542 int index = cast(int)x; 543 float delta = x - index; // delta is 0..1 544 float s0 = getSample(index, 0); 545 float s1 = getSample(index + 1, 0); 546 buf.ptr[i] = s0 * (1 - delta) + s1 * delta; 547 x += step; 548 } 549 } 550 551 WaveFile getRange(int start, int end, bool forceMono = false) { 552 int len = end - start; 553 if (len <= 0 || start < 0 || end > frames) 554 return null; 555 WaveFile res = new WaveFile(); 556 res.channels = forceMono ? 1 : channels; 557 res.sampleRate = sampleRate; 558 res.frames = len; 559 res.data = new float[len * res.channels]; 560 if (res.channels == channels) 561 res.data[0..$] = data[start * channels .. (start + len) * channels]; 562 else { 563 // copy only one channel 564 for (int i = 0; i < res.frames; i++) { 565 res.data[i] = data[start * channels + i * channels]; 566 } 567 } 568 return res; 569 } 570 571 WaveFile upsample4x(int start, int end) { 572 if (start < 0) 573 start = 0; 574 if (end > frames) 575 end = frames; 576 int dstart = 16; 577 int dend = 16; 578 if (start - dstart < 0) 579 dstart = start; 580 if (end + dend > frames) 581 dend = frames - end; 582 WaveFile tmp = getRange(start - dstart, end + dend); 583 WaveFile resampled = tmp.upsample4x(); 584 return resampled.getRange(dstart * 4, resampled.frames - dend * 4); 585 } 586 587 // frequency domain z padding 588 WaveFile upsample4x() { 589 import std.math : PI, sin; 590 591 float maxsrc = 0; 592 float maxdst = 0; 593 for (int i = 0; i < frames; i++) 594 if (maxsrc < data[i]) 595 maxsrc = data[i]; 596 597 WaveFile res = new WaveFile(); 598 res.channels = 1; 599 res.sampleRate = sampleRate * 4; 600 res.frames = frames * 4; 601 res.data = new float[res.frames]; 602 res.data[0..$] = 0; 603 // upsampling with zero stuffing 604 for (int i = 0; i < frames; i++) { 605 res.data[i * 4] = data[i] * 4; 606 res.data[i * 4 + 1] = 0; 607 res.data[i * 4 + 2] = 0; 608 res.data[i * 4 + 3] = 0; 609 } 610 // filtering 611 int windowSize = 31; 612 float sincN = 5.8f; 613 float[] window = blackmanWindow(windowSize); 614 float[] sinc = new float[windowSize]; 615 for (int i = 0; i < windowSize; i++) { 616 float x = i - windowSize / 2; 617 x = x / sincN; 618 float v = i == windowSize / 2 ? 1 : sin(PI * x) / (PI * x); 619 sinc[i] = v * window[i]; 620 } 621 normalizeFirFilter(sinc); 622 float[] tmp = res.data.dup; 623 applyFirFilter(tmp, res.data, sinc); 624 625 626 for (int i = 0; i < res.frames; i++) 627 if (maxdst < res.data[i]) 628 maxdst = res.data[i]; 629 import dlangui.core.logger; 630 Log.d("max value before upsampling: ", maxsrc, " after: ", maxdst); 631 632 return res; 633 } 634 635 // frequency domain z padding 636 WaveFile upsample4xFreqDomainZpad() { 637 import std.numeric; 638 import std.complex; 639 immutable int BATCH_SIZE = 256; 640 immutable int overlap = BATCH_SIZE/4; 641 Fft fft = new Fft(BATCH_SIZE); 642 Fft ifft = new Fft(BATCH_SIZE * 4); 643 WaveFile res = new WaveFile(); 644 res.channels = 1; 645 res.sampleRate = sampleRate * 4; 646 res.frames = frames * 4; 647 res.data = new float[res.frames]; 648 res.data[0..$] = 0; 649 float[] srcBuf = new float[BATCH_SIZE]; 650 Complex!double[] invBuf = new Complex!double[BATCH_SIZE]; 651 Complex!double[] invBufResampled = new Complex!double[BATCH_SIZE * 4]; 652 Complex!double[] fftFrame = new Complex!double[BATCH_SIZE * 4]; 653 654 float maxsrc = 0; 655 float maxdst = 0; 656 for (int i = 0; i < frames; i++) 657 if (maxsrc < data[i]) 658 maxsrc = data[i]; 659 660 for (int x = -overlap; x < frames; x += BATCH_SIZE - overlap * 2) { 661 getSamples(x, 0, srcBuf); 662 //immutable int SMOOTH_RANGE = BATCH_SIZE / 10; 663 //for (int i = 0; i < SMOOTH_RANGE; i++) { 664 // float k = 1.0f - i / cast(float)SMOOTH_RANGE; 665 // float s1 = srcBuf[i * 4]; 666 // float s2 = srcBuf[BATCH_SIZE - i * 4 - 4]; 667 // float sm = (s1 + s2) / 2; 668 // srcBuf[i * 4] = s1 * (1 - k) + sm * k; 669 // srcBuf[BATCH_SIZE - i * 4 - 4] = s2 * (1 - k) + sm * k; 670 //} 671 fft.fft(srcBuf, invBuf); 672 invBufResampled[0 .. $] = Complex!double(0,0); 673 for (int i = 0; i < BATCH_SIZE / 2; i++) { 674 invBufResampled[i] = invBuf[i]; 675 } 676 for (int i = 0; i < BATCH_SIZE / 2 - 1; i++) { 677 invBufResampled[$ - i - 1] = invBuf[$ - i - 1]; 678 } 679 invBufResampled[$ / 2] = invBuf[$ / 2]; 680 ifft.inverseFft(invBufResampled, fftFrame); 681 for (int i = 0; i < BATCH_SIZE * 4 - overlap * 8; i++) { 682 int dstIndex = x * 4 + i + overlap * 4; 683 if (dstIndex >= 0 && dstIndex < res.frames) 684 res.data.ptr[dstIndex] = fftFrame[i + overlap * 4].re * 4; 685 } 686 } 687 688 for (int i = 0; i < res.frames; i++) 689 if (maxdst < res.data[i]) 690 maxdst = res.data[i]; 691 692 import dlangui.core.logger; 693 694 Log.d("max value before upsampling: ", maxsrc, " after: ", maxdst); 695 696 return res; 697 } 698 699 // time domain z padding upsampling 700 WaveFile upsample4xTimeDomainZpad() { 701 import std.numeric; 702 import std.complex; 703 immutable int BATCH_SIZE = 256; 704 immutable int overlap = BATCH_SIZE/4; 705 Fft fft = new Fft(BATCH_SIZE); 706 Fft ifft = new Fft(BATCH_SIZE); 707 WaveFile res = new WaveFile(); 708 res.channels = 1; 709 res.sampleRate = sampleRate * 4; 710 res.frames = frames * 4; 711 res.data = new float[res.frames]; 712 res.data[0..$] = 0; 713 float[] srcBuf = new float[BATCH_SIZE]; 714 Complex!double[] fftFrame = new Complex!double[BATCH_SIZE]; 715 Complex!double[] invBuf = new Complex!double[BATCH_SIZE]; 716 717 float maxsrc = 0; 718 float maxdst = 0; 719 for (int i = 0; i < frames; i++) 720 if (maxsrc < data[i]) 721 maxsrc = data[i]; 722 723 for (int x = -overlap / 4; x < frames; x += BATCH_SIZE / 4 - overlap / 2) { 724 getSamplesZpad4(x, 0, srcBuf); 725 //immutable int SMOOTH_RANGE = BATCH_SIZE / 10; 726 //for (int i = 0; i < SMOOTH_RANGE; i++) { 727 // float k = 1.0f - i / cast(float)SMOOTH_RANGE; 728 // float s1 = srcBuf[i * 4]; 729 // float s2 = srcBuf[BATCH_SIZE - i * 4 - 4]; 730 // float sm = (s1 + s2) / 2; 731 // srcBuf[i * 4] = s1 * (1 - k) + sm * k; 732 // srcBuf[BATCH_SIZE - i * 4 - 4] = s2 * (1 - k) + sm * k; 733 //} 734 fft.fft(srcBuf, invBuf); 735 for (int i = BATCH_SIZE / 8 - 5; i <= BATCH_SIZE * 7 / 8 + 5 + 1; i++) { 736 invBuf[i] = Complex!double(0,0); 737 } 738 // smoother lowpass filter 739 /* 740 immutable int SMOOTH_DIST = BATCH_SIZE / 30; 741 for (int i = 1; i < SMOOTH_DIST; i++) { 742 float k = SMOOTH_DIST / i; 743 invBuf[BATCH_SIZE / 8 - i].re *= k; 744 invBuf[BATCH_SIZE / 8 - i].im *= k; 745 invBuf[BATCH_SIZE * 7 / 8 + i].re *= k; 746 invBuf[BATCH_SIZE * 7 / 8 + i].im *= k; 747 } 748 */ 749 ifft.inverseFft(invBuf, fftFrame); 750 for (int i = overlap; i < BATCH_SIZE - overlap; i++) { 751 int index = x * 4 + i; 752 if (index >= 0 && index < res.frames) 753 res.data.ptr[index] = fftFrame[i].re * 4; 754 } 755 } 756 757 for (int i = 0; i < res.frames; i++) 758 if (maxdst < res.data[i]) 759 maxdst = res.data[i]; 760 761 import dlangui.core.logger; 762 763 Log.d("max value before upsampling: ", maxsrc, " after: ", maxdst); 764 765 return res; 766 } 767 768 /// returns autocorrelation best frequency at center of file 769 float calcBaseFrequency() { 770 return calcLocalFrequency(frameToTime(frames / 2), 30); 771 } 772 773 int getMaxAmplitudeSign() { 774 float maxPositive = 0; 775 float maxNegative = 0; 776 foreach(v; data) { 777 if (v > 0 && maxPositive < v) 778 maxPositive = v; 779 if (v < 0 && maxNegative < -v) 780 maxNegative = -v; 781 } 782 return maxPositive > maxNegative ? 1 : -1; 783 } 784 785 float findZeroCrossingNear(float timePosition) { 786 int frame = timeToFrame(timePosition); 787 788 int frame0, frame1; 789 frame0 = frame1 = frame; 790 for (int i = 0; i < 10000; i++) { 791 float sample = getSample(frame + i); 792 float sample1 = getSample(frame + i + 1); 793 if ((sample >= 0 && sample1 <=0) || (sample <= 0 && sample1 >=0)) { 794 frame0 = frame + i; 795 frame1 = frame + i + 1; 796 break; 797 } 798 sample = getSample(frame - i); 799 sample1 = getSample(frame - i - 1); 800 if ((sample >= 0 && sample1 <=0) || (sample <= 0 && sample1 >=0)) { 801 frame0 = frame - i - 1; 802 frame1 = frame - i; 803 break; 804 } 805 } 806 807 if (frame0 < frame1) { 808 float sample0 = getSample(frame0); 809 float sample1 = getSample(frame1); 810 if (sample0 != sample1) { 811 float time0 = frameToTime(frame0); 812 float time1 = frameToTime(frame1); 813 // linear equation, find crossing 0 814 float t = time0 - sample0 * (time1 - time0) / (sample1 - sample0); 815 return t; 816 } 817 return frameToTime(frame0); 818 } 819 // cannot correct 820 return timePosition; 821 } 822 823 float findMaxDiffPosition(int startFrame, int endFrame) { 824 float maxDiffPlus = 0; 825 float maxDiffMinus = 0; 826 int maxDiffPositionPlus = startFrame; 827 int maxDiffPositionMinus = startFrame; 828 for (int i = startFrame + 1; i < endFrame; i++) { 829 float diff = getSample(i) - getSample(i - 1); 830 if (diff > 0) { 831 if (maxDiffPlus < diff) { 832 maxDiffPlus = diff; 833 maxDiffPositionPlus = i; 834 } 835 } else { 836 if (maxDiffMinus < -diff) { 837 maxDiffMinus = -diff; 838 maxDiffPositionMinus = i; 839 } 840 } 841 } 842 if (maxDiffPlus > maxDiffMinus) 843 return frameToTime(maxDiffPositionPlus); 844 return frameToTime(maxDiffPositionMinus); 845 } 846 847 /// sign == 1 if biggest amplitude is positive, -1 if negative 848 float[] findZeroCrossingPositions(int sign = 1) { 849 import dlangui.core.logger; 850 851 float freqPosition = frameToTime(frames / 2); 852 float freqBigRange = calcLocalFrequency(freqPosition, 40); 853 float freqShortRange = calcLocalFrequency(freqPosition, freqBigRange / 4); 854 Log.d("Frequencies: ", freqBigRange, " ", freqShortRange); 855 856 freqBigRange = calcLocalFrequency(freqPosition, 40); 857 freqShortRange = calcLocalFrequency(freqPosition, freqBigRange / 4); 858 859 float zeroPhaseTime; 860 float amplitude; 861 //findNearPhase0(freqPosition, freqShortRange, 1, zeroPhaseTime, amplitude); 862 float zeroPhaseTime2; 863 float amplitude2; 864 float freq; 865 866 // find max difference position 867 freqPosition = findMaxDiffPosition(timeToFrame(freqPosition - 1/freqShortRange), timeToFrame(freqPosition + 1/freqShortRange)); 868 // find near zero crossing 869 freqPosition = findZeroCrossingNear(freqPosition); 870 871 float initialPosition = freqPosition; 872 float initialFreq = freqShortRange; 873 float[] zpositions; 874 float[] zpositionsBefore; 875 zpositions ~= freqPosition; 876 877 float maxtime = frames / cast(float)sampleRate - 1f / initialFreq; 878 float mintime = 1f / initialFreq; 879 880 881 Log.d("Scanning time range ", initialPosition, " .. ", maxtime); 882 freqPosition = initialPosition; 883 freq = initialFreq; 884 float step = 1/freq; 885 float prevPosition = freqPosition; 886 freqPosition += step; 887 for (int i = 0; i < 10000; i++) { 888 if (freqPosition > maxtime) 889 break; 890 freqPosition = findZeroCrossingNear(freqPosition); 891 step = freqPosition - prevPosition; 892 prevPosition = freqPosition; 893 zpositions ~= freqPosition; 894 freqPosition += step; 895 } 896 // till end - just use fixed freq 897 maxtime = frames / cast(float)sampleRate - 0.1f / initialFreq; 898 for (int i = 0; i < 3; i++) { 899 if (freqPosition > maxtime) 900 break; 901 zpositions ~= freqPosition; 902 freqPosition += step; 903 } 904 Log.d("Scanning time range ", mintime, " .. ", initialPosition); 905 freqPosition = initialPosition; 906 freq = initialFreq; 907 step = 1/freq; 908 prevPosition = freqPosition; 909 freqPosition -= step; 910 for (int i = 0; i < 10000; i++) { 911 if (freqPosition < mintime) 912 break; 913 freqPosition = findZeroCrossingNear(freqPosition); 914 step = prevPosition - freqPosition; 915 prevPosition = freqPosition; 916 zpositionsBefore ~= freqPosition; 917 freqPosition -= step; 918 } 919 // till beginning - just use fixed freq 920 mintime = 0.1f / initialFreq; 921 for (int i = 0; i < 3; i++) { 922 if (freqPosition < mintime) 923 break; 924 zpositionsBefore ~= freqPosition; 925 freqPosition -= step; 926 } 927 928 //Log.d("zpositions: ", zpositions, " before: ", zpositionsBefore); 929 930 // combine positions before and after 931 float[] zpositionsAll; 932 for (int i = cast(int)zpositionsBefore.length - 1; i >= 0; i--) 933 zpositionsAll ~= zpositionsBefore[i]; 934 zpositionsAll ~= zpositions; 935 936 float[] freqs; 937 for (int i = 1; i < zpositionsAll.length; i++) { 938 freqs ~= 1 / (zpositionsAll[i] - zpositionsAll[i - 1]); 939 } 940 941 //Log.d("zpositionsAll: ", zpositionsAll); 942 Log.d("freqs: ", freqs); 943 return zpositionsAll; 944 } 945 946 947 948 /// sign == 1 if biggest amplitude is positive, -1 if negative 949 float[] findZeroPhasePositions(int sign = 1) { 950 import dlangui.core.logger; 951 952 float freqPosition = frameToTime(frames / 2); 953 float freqBigRange = calcLocalFrequency(freqPosition, 40); 954 float freqShortRange = calcLocalFrequency(freqPosition, freqBigRange / 4); 955 Log.d("Frequencies: ", freqBigRange, " ", freqShortRange); 956 957 freqBigRange = calcLocalFrequency(freqPosition, 40); 958 freqShortRange = calcLocalFrequency(freqPosition, freqBigRange / 4); 959 960 float zeroPhaseTime; 961 float amplitude; 962 //findNearPhase0(freqPosition, freqShortRange, 1, zeroPhaseTime, amplitude); 963 float zeroPhaseTime2; 964 float amplitude2; 965 float freq; 966 findNearPhase0FreqAutocorrection(freqPosition, freqShortRange, 1/freqShortRange, sign, zeroPhaseTime2, amplitude2, freq); 967 Log.d("Zero phase near ", freqPosition, " at freq ", freqShortRange, " : pos=", zeroPhaseTime2, " amp=", amplitude2, " freq=", freq); 968 freqPosition = zeroPhaseTime2; 969 970 971 float initialPosition = freqPosition; 972 float initialFreq = freq; 973 float[] zpositions; 974 float[] zpositionsBefore; 975 zpositions ~= freqPosition; 976 977 float maxtime = frames / cast(float)sampleRate - 2f / initialFreq; 978 float mintime = 2f / initialFreq; 979 980 981 Log.d("Scanning time range ", initialPosition, " .. ", maxtime); 982 freqPosition = initialPosition; 983 freq = initialFreq; 984 float step = 1/freq; 985 freqPosition += step; 986 for (int i = 0; i < 10000; i++) { 987 if (freqPosition > maxtime) 988 break; 989 float oldFreq = freq; 990 findNearPhase0FreqAutocorrection(freqPosition, oldFreq, 0.2 / oldFreq, sign, zeroPhaseTime2, amplitude2, freq); 991 Log.d("Zero phase near ", freqPosition, " at freq ", oldFreq, " : pos=", zeroPhaseTime2, " amp=", amplitude2, " step=", step, " freq=", freq); 992 step = 1/freq; 993 freqPosition = zeroPhaseTime2; 994 zpositions ~= freqPosition; 995 freqPosition += step; 996 } 997 // till end - just use fixed freq 998 maxtime = frames / cast(float)sampleRate - 0.1f / initialFreq; 999 for (int i = 0; i < 3; i++) { 1000 if (freqPosition > maxtime) 1001 break; 1002 zpositions ~= freqPosition; 1003 freqPosition += step; 1004 } 1005 Log.d("Scanning time range ", mintime, " .. ", initialPosition); 1006 freqPosition = initialPosition; 1007 freq = initialFreq; 1008 step = 1/freq; 1009 freqPosition -= step; 1010 for (int i = 0; i < 10000; i++) { 1011 if (freqPosition < mintime) 1012 break; 1013 float oldFreq = freq; 1014 findNearPhase0FreqAutocorrection(freqPosition, oldFreq, 0.2 / oldFreq, sign, zeroPhaseTime2, amplitude2, freq); 1015 Log.d("Zero phase near ", freqPosition, " at freq ", oldFreq, " : pos=", zeroPhaseTime2, " amp=", amplitude2, " step=", step, " freq=", freq); 1016 step = 1/freq; 1017 freqPosition = zeroPhaseTime2; 1018 zpositionsBefore ~= freqPosition; 1019 freqPosition -= step; 1020 } 1021 // till beginning - just use fixed freq 1022 mintime = 0.1f / initialFreq; 1023 for (int i = 0; i < 3; i++) { 1024 if (freqPosition < mintime) 1025 break; 1026 zpositionsBefore ~= freqPosition; 1027 freqPosition -= step; 1028 } 1029 1030 //Log.d("zpositions: ", zpositions, " before: ", zpositionsBefore); 1031 1032 float[] zpositionsAll; 1033 for (int i = cast(int)zpositionsBefore.length - 1; i >= 0; i--) 1034 zpositionsAll ~= zpositionsBefore[i]; 1035 zpositionsAll ~= zpositions; 1036 1037 float[] freqs; 1038 for (int i = 1; i < zpositionsAll.length; i++) { 1039 freqs ~= 1 / (zpositionsAll[i] - zpositionsAll[i - 1]); 1040 } 1041 1042 //Log.d("zpositionsAll: ", zpositionsAll); 1043 Log.d("freqs: ", freqs); 1044 return zpositionsAll; 1045 } 1046 1047 float calcLocalFrequency(float position, float minFreq) { 1048 import dlangui.core.logger; 1049 int pos = timeToFrame(position); 1050 int windowSize = timeToFrame(1 / minFreq) * 2; 1051 windowSize &= 0xFFFFFE; 1052 float[] window = blackmanWindow(windowSize); 1053 float[] frame = new float[windowSize + 1]; 1054 float[] corr = new float[windowSize + 1]; 1055 float[] windowcorr = new float[windowSize + 1]; 1056 getSamples(pos - windowSize / 2, 0, frame); 1057 for (int i = 0; i < frame.length; i++) { 1058 frame[i] *= window[i]; 1059 } 1060 correlation(frame, frame, corr); 1061 correlation(window, window, windowcorr); 1062 1063 int p = 1; 1064 for (; p < corr.length; p++) { 1065 if (corr[p] < 0) 1066 break; 1067 } 1068 //Log.d("Negative correlation at offset ", p); 1069 float maxcorr = corr[p]; 1070 int maxcorrpos = p; 1071 for (; p < corr.length; p++) { 1072 if (maxcorr < corr[p]) { 1073 maxcorr = corr[p]; 1074 maxcorrpos = p; 1075 } 1076 } 1077 float a, b, c; 1078 float correction1 = windowcorr[maxcorrpos-1]; 1079 float correction2 = windowcorr[maxcorrpos]; 1080 float correction3 = windowcorr[maxcorrpos+1]; 1081 correction1 = correction2 = correction3 = 1; 1082 calcParabola(maxcorrpos - 1, 1083 corr[maxcorrpos-1] / correction1, 1084 corr[maxcorrpos] / correction2, 1085 corr[maxcorrpos+1] / correction3, a, b, c); 1086 float x0 = -b / (2 * a); 1087 //Log.d("Max correlation = ", maxcorr, " at offset ", maxcorrpos, " corr0: ", corr[0], " approx best position = ", x0); 1088 float period = frameToTime(x0); 1089 if (period > 0) 1090 return 1 / period; 1091 return 0; 1092 } 1093 1094 void findNearPhase0FreqAutocorrection(float positionTimeSeconds, float freqHerz, float maxCorrection, int sign, ref float zeroPhaseTimePosSeconds, ref float amplitude, ref float newFreq) { 1095 import dlangui.core.logger; 1096 float startPosition = positionTimeSeconds; 1097 findNearPhase0(positionTimeSeconds, freqHerz, sign, zeroPhaseTimePosSeconds, amplitude); 1098 float autocorrelatedFrequency = calcLocalFrequency(zeroPhaseTimePosSeconds, freqHerz / 2); 1099 //Log.d("Initial phase detection: freq=", freqHerz, " pos=", positionTimeSeconds, " => ", zeroPhaseTimePosSeconds, " amp=", amplitude, " autocorrFreq=", autocorrelatedFrequency); 1100 if (autocorrelatedFrequency >= freqHerz * 0.8 && autocorrelatedFrequency <= freqHerz * 1.2) { 1101 freqHerz = autocorrelatedFrequency; 1102 positionTimeSeconds = zeroPhaseTimePosSeconds; 1103 findNearPhase0(positionTimeSeconds, freqHerz, sign, zeroPhaseTimePosSeconds, amplitude); 1104 //Log.d(" position updated (1) for freq=", freqHerz, " pos=", positionTimeSeconds, " => ", zeroPhaseTimePosSeconds, " amp=", amplitude); 1105 } 1106 if (zeroPhaseTimePosSeconds >= startPosition - maxCorrection && zeroPhaseTimePosSeconds <= startPosition + maxCorrection) { 1107 positionTimeSeconds = zeroPhaseTimePosSeconds; 1108 newFreq = freqHerz; 1109 } else { 1110 zeroPhaseTimePosSeconds = startPosition; 1111 newFreq = freqHerz; 1112 } 1113 // calculate again at zero phase position 1114 //findNearPhase0(positionTimeSeconds, freqHerz, sign, zeroPhaseTimePosSeconds, amplitude); 1115 //Log.d(" position updated (2) for freq=", freqHerz, " pos=", positionTimeSeconds, " => ", zeroPhaseTimePosSeconds, " amp=", amplitude); 1116 //newFreq = freqHerz; 1117 } 1118 1119 void findNearPhase0(float positionTimeSeconds, float freqHerz, int sign, ref float zeroPhaseTimePosSeconds, ref float amplitude) { 1120 import std.math : sqrt, atan2, PI; 1121 immutable int TABLE_LEN = 256; 1122 // immutable int TABLE_LEN = 768; 1123 float[TABLE_LEN] buf; 1124 float x = positionTimeSeconds * sampleRate; 1125 float periodSamples = sampleRate / freqHerz; 1126 // get interpolated one period of freq 1127 getSamplesInterpolated(x, periodSamples / TABLE_LEN, buf[0..$]); 1128 double sumSin = 0; 1129 double sumCos = 0; 1130 //double sumSin2 = 0; 1131 //double sumCos2 = 0; 1132 //float stepWidth = periodSamples / 256; 1133 //for (int i = 0; i < TABLE_LEN; i++) { 1134 // sumSin += sign * buf.ptr[i] * SIN_SYNC_TABLE_768.ptr[i]; 1135 // sumCos += sign * buf.ptr[i] * SIN_SYNC_TABLE_768.ptr[i]; 1136 // //sumSin2 += SIN_SYNC_TABLE_768.ptr[i] * SIN_SYNC_TABLE_768.ptr[i]; 1137 // //sumCos2 += SIN_SYNC_TABLE_768.ptr[i] * SIN_SYNC_TABLE_768.ptr[i]; 1138 //} 1139 for (int i = 0; i < TABLE_LEN; i++) { 1140 sumSin += sign * buf.ptr[i] * SIN_TABLE_256.ptr[i]; 1141 sumCos += sign * buf.ptr[i] * COS_TABLE_256.ptr[i]; 1142 //sumSin2 += SIN_TABLE_256.ptr[i] * SIN_TABLE_256.ptr[i]; 1143 //sumCos2 += COS_TABLE_256.ptr[i] * SIN_TABLE_256.ptr[i]; 1144 } 1145 sumSin /= TABLE_LEN; 1146 sumCos /= TABLE_LEN; 1147 //sumSin /= (periodSamples * periodSamples); 1148 //sumCos /= (periodSamples * periodSamples); 1149 // calc amplitude 1150 amplitude = sqrt(sumSin * sumSin + sumCos * sumCos); // / periodSamples; 1151 // normalize 1152 //sumSin /= amplitude; 1153 //sumCos /= amplitude; 1154 float phase = atan2(sumSin, sumCos) - PI/2; 1155 if (phase < -PI) 1156 phase += PI * 2; 1157 //float phase2 = atan2(sumSin2, sumCos2) - PI/2; 1158 float zeroPhaseX = x + periodSamples * phase / (2 * PI); 1159 zeroPhaseTimePosSeconds = zeroPhaseX / sampleRate; 1160 } 1161 1162 float sinAmpAt(float positionTimeSeconds, float freqHerz, int sign) { 1163 import std.math : sqrt, atan2, PI; 1164 float[256] buf; 1165 float x = positionTimeSeconds * sampleRate; 1166 float periodSamples = sampleRate / freqHerz; 1167 // get interpolated one period of freq 1168 getSamplesInterpolated(x, periodSamples / 256, buf[0..$]); 1169 float sumSin = 0; 1170 //float stepWidth = periodSamples / 256; 1171 for (int i = 0; i < 256; i++) { 1172 sumSin += sign * buf.ptr[i] * SIN_TABLE_256.ptr[i]; 1173 } 1174 sumSin /= 256; 1175 //sumSin /= (periodSamples * periodSamples); 1176 //sumCos /= (periodSamples * periodSamples); 1177 // calc amplitude 1178 //float amplitude = sqrt(sumSin * sumSin); // / periodSamples; 1179 return sumSin; //amplitude; 1180 } 1181 1182 /// create WaveFile which is FIR filtered copy of current wave 1183 WaveFile firFilter(float[] filter) { 1184 WaveFile res = new WaveFile(this, true); 1185 applyFirFilter(data, res.data, filter); 1186 return res; 1187 } 1188 1189 /// create WaveFile which is FIR inverse filtered copy of current wave 1190 WaveFile firFilterInverse(float[] filter) { 1191 WaveFile res = new WaveFile(this, true); 1192 applyFirFilterInverse(data, res.data, filter); 1193 return res; 1194 } 1195 1196 } 1197 1198 /* Autocorrelation LPC coeff generation algorithm invented by 1199 N. Levinson in 1947, modified by J. Durbin in 1959. */ 1200 1201 /* Input : elements of time doamin data (with window applied) 1202 Output: lpc coefficients, excitation energy */ 1203 float calcLPC(float[] data, float[] lpci) { 1204 int n = cast(int)data.length; 1205 int m = cast(int)lpci.length; 1206 double[] aut = new double[m + 1]; 1207 double[] lpc = new double[m]; 1208 double error; 1209 double epsilon; 1210 int i,j; 1211 1212 /* autocorrelation, p+1 lag coefficients */ 1213 j=m+1; 1214 while(j--){ 1215 double d=0; /* double needed for accumulator depth */ 1216 for(i=j; i<n; i++) 1217 d += cast(double)data[i] * data[i-j]; 1218 aut[j]=d; 1219 } 1220 1221 /* Generate lpc coefficients from autocorr values */ 1222 1223 /* set our noise floor to about -100dB */ 1224 error = aut[0] * (1. + 1e-10); 1225 epsilon = 1e-9*aut[0]+1e-10; 1226 1227 for(i=0; i<m; i++){ 1228 double r= -aut[i+1]; 1229 1230 if (error < epsilon) { 1231 lpc[i .. $] = 0; 1232 break; 1233 } 1234 1235 /* Sum up this iteration's reflection coefficient; note that in 1236 Vorbis we don't save it. If anyone wants to recycle this code 1237 and needs reflection coefficients, save the results of 'r' from 1238 each iteration. */ 1239 1240 for (j=0; j<i; j++) 1241 r -= lpc[j] * aut[i - j]; 1242 r /= error; 1243 1244 /* Update LPC coefficients and total error */ 1245 1246 lpc[i]=r; 1247 for(j=0; j<i/2; j++) { 1248 double tmp = lpc[j]; 1249 1250 lpc[j] += r * lpc[i-1-j]; 1251 lpc[i-1-j] += r*tmp; 1252 } 1253 if(i&1) 1254 lpc[j]+=lpc[j]*r; 1255 1256 error *= 1. - r*r; 1257 } 1258 1259 /* slightly damp the filter */ 1260 { 1261 double g = .99; 1262 double damp = g; 1263 for(j=0;j<m;j++){ 1264 lpc[j]*=damp; 1265 damp*=g; 1266 } 1267 } 1268 1269 for(j=0; j<m; j++) 1270 lpci[j] = cast(float)lpc[j]; 1271 1272 /* we need the error value to know how big an impulse to hit the 1273 filter with later */ 1274 1275 return error; 1276 } 1277 1278 void applyFirFilterInverse(float[] src, float[] dst, float[] filter) { 1279 assert(src.length == dst.length); 1280 applyFirFilter(src, dst, filter); 1281 int len = cast(int)src.length; 1282 for(int i = 0; i < len; i++) { 1283 dst.ptr[i] = src.ptr[i] - dst.ptr[i]; 1284 } 1285 } 1286 1287 void applyFirFilter(float[] src, float[] dst, float[] filter) { 1288 assert(src.length == dst.length); 1289 int filterLen = cast(int)filter.length; 1290 int filterMiddle = filterLen / 2; 1291 int len = cast(int)src.length; 1292 for (int x = 0; x < len; x++) { 1293 double filterSum = 0; 1294 double resultSum = 0; 1295 for (int i = 0; i < filterLen; i++) { 1296 int index = i - filterMiddle + x; 1297 if (index >= 0 && index < len) { 1298 float sample = src.ptr[index]; 1299 float flt = filter.ptr[i]; 1300 resultSum += sample * flt; 1301 filterSum += flt; 1302 } 1303 } 1304 dst.ptr[x] = resultSum / filterSum; 1305 } 1306 } 1307 1308 // calc parabola coefficients for points (x1, y1), (x1 + 1, y2), (x1 + 2, y3) 1309 void calcParabola(int x1, float y1, float y2, float y3, ref float a, ref float b, ref float c) { 1310 a = (y3 + y1) / 2 - y2; 1311 b = y2 - y1 - a * (2 * x1 + 1); 1312 c = (x1 + 1) *y1 - x1 * y2 + a * x1 * (x1 + 1); 1313 } 1314 1315 float[] makeLowpassBlackmanFirFilter(int N) { 1316 float[] res = blackmanWindow(N).dup; 1317 double s = 0; 1318 foreach(v; res) 1319 s += v; 1320 foreach(ref v; res) 1321 v = -v / s; 1322 res[N / 2] += 2; 1323 1324 1325 s = 0; 1326 foreach(v; res) 1327 s += v; 1328 1329 return res; 1330 } 1331 1332 void normalizeFirFilter(float[] coefficients) { 1333 double s = 0; 1334 foreach(v; coefficients) 1335 s += v; 1336 foreach(ref v; coefficients) 1337 v /= s; 1338 } 1339 1340 __gshared float[][int] BLACKMAN_WINDOW_CACHE; 1341 1342 // generate blackman window in array [0..N] (value at N/2 == 1) 1343 float[] blackmanWindow(int N) { 1344 if (auto existing = N in BLACKMAN_WINDOW_CACHE) { 1345 return *existing; 1346 } 1347 import std.math : cos, PI; 1348 float[] res = new float[N + 1]; 1349 for (int i = 1; i <= N + 1; i++) { 1350 res[i - 1] = 0.42f - 0.5f * cos(2 * PI * i / (N + 2)) + 0.08 * cos(4 * PI * i / (N + 2)); 1351 } 1352 BLACKMAN_WINDOW_CACHE[N] = res; 1353 return res; 1354 } 1355 1356 void correlation(float[] a, float[] b, float[] res) { 1357 assert(a.length == b.length); 1358 assert(res.length >= a.length); 1359 for (int diff = 0; diff < a.length; diff++) { 1360 float sum = 0; 1361 for (int i = 0; i < a.length - diff; i++) 1362 sum += a[i] * b[i + diff]; 1363 res[diff] = sum; 1364 } 1365 } 1366 1367 void smoothTimeMarks(ref float[] marks) { 1368 if (marks.length < 3) 1369 return; 1370 float[] tmp = new float[marks.length]; 1371 tmp[0] = marks[0]; 1372 tmp[$ - 1] = marks[$ - 1]; 1373 for (int i = 1; i + 1 < marks.length; i++) { 1374 tmp[i] = (marks[i] + marks[i - 1] + marks[i + 1]) / 3; 1375 } 1376 marks = tmp; 1377 } 1378 1379 void smoothTimeMarksShifted(ref float[] marks, ref float[] negativeMarks) { 1380 float[] tmpPositive = smoothTimeMarks(marks, negativeMarks); 1381 float[] tmpNegative = smoothTimeMarks(negativeMarks, marks); 1382 marks = tmpPositive; 1383 negativeMarks = tmpNegative; 1384 } 1385 1386 /// smooth phase marks using marks shifted by half period 1387 float[] smoothTimeMarks(float[] marks, float[] marksShifted) { 1388 if (marks.length < 3 || marksShifted.length < 3) 1389 return marks; 1390 int i = 0; 1391 int i0 = 0; 1392 if (marks[i + i0] < marksShifted[i]) { 1393 i = 1; 1394 i0 = -1; 1395 } 1396 float[] tmp = marks.dup; 1397 for (; i < marks.length && i + i0 + 1 < marksShifted.length; i++) { 1398 if (marks[i] > marksShifted[i + i0] && marks[i] < marksShifted[i + i0 + 1]) 1399 tmp[i] = (marks[i] + marksShifted[i + i0] + marksShifted[i + i0 + 1]) / 3; 1400 } 1401 if (i0 == -1) { 1402 // smooth first item 1403 tmp[0] = (tmp[0] + tmp[1] - (tmp[2] - tmp[1])) / 2; 1404 } 1405 if (tmp.length - 1 + i0 < marksShifted.length - 1) { 1406 // smooth last item 1407 tmp[$ - 1] = (tmp[$ - 1] + tmp[$ - 2] + (tmp[$ - 2] - tmp[$ - 3])) / 2; 1408 } 1409 return tmp; 1410 } 1411 1412 // generate several periods sin table multiplied by blackman window, with phase 0 at middle point 1413 float[] generateSyncSinTable(int len, int periods) { 1414 float[] singlePeriod = generateSinTable(len); 1415 float[] full; 1416 for (int i = 0; i < periods; i++) 1417 full ~= singlePeriod; 1418 float[] window = blackmanWindow(len * periods); 1419 double s = 0; 1420 for(int i = 0; i < full.length; i++) { 1421 full[i] = full[i] * window[i]; 1422 s += full[i]; 1423 } 1424 s /= full.length; 1425 for(int i = 0; i < full.length; i++) { 1426 full[i] -= s; 1427 } 1428 return full; 1429 } 1430 1431 // generate several periods sin table multiplied by blackman window, with phase 0 at middle point 1432 float[] generateSyncCosTable(int len, int periods) { 1433 float[] singlePeriod = generateCosTable(len); 1434 float[] full; 1435 for (int i = 0; i < periods; i++) 1436 full ~= singlePeriod; 1437 float[] window = blackmanWindow(len * periods); 1438 double s = 0; 1439 for(int i = 0; i < full.length; i++) { 1440 full[i] = full[i] * window[i]; 1441 s += full[i]; 1442 } 1443 s /= full.length; 1444 for(int i = 0; i < full.length; i++) { 1445 full[i] -= s; 1446 } 1447 return full; 1448 } 1449 1450 // generate sin table with phase 0 at middle point 1451 float[] generateSinTable(int len) { 1452 float[] res = new float[len]; 1453 import std.math : sin, cos, PI; 1454 for (int i = 0; i < len; i++) { 1455 double x = (i - len / 2) * 2 * PI / len; 1456 res[i] = sin(x); 1457 } 1458 return res; 1459 } 1460 1461 // generate sin table with phase 0 at middle point 1462 float[] generateCosTable(int len) { 1463 float[] res = new float[len]; 1464 import std.math : sin, cos, PI; 1465 for (int i = 0; i < len; i++) { 1466 double x = (i - len / 2) * 2 * PI / len; 1467 res[i] = cos(x); 1468 } 1469 return res; 1470 } 1471 1472 __gshared float[] SIN_TABLE_256; 1473 __gshared float[] COS_TABLE_256; 1474 __gshared float[] SIN_SYNC_TABLE_768; 1475 __gshared float[] COS_SYNC_TABLE_768; 1476 1477 __gshared static this() { 1478 SIN_TABLE_256 = generateSinTable(256); 1479 COS_TABLE_256 = generateCosTable(256); 1480 SIN_SYNC_TABLE_768 = generateSyncSinTable(256, 3); 1481 COS_SYNC_TABLE_768 = generateSyncCosTable(256, 3); 1482 } 1483 1484 1485 immutable float LPC_SCALING = 1.0f; 1486 immutable float FREQ_SCALE = 1.0f; 1487 //#define X2ANGLE(x) (acos(x)) 1488 immutable float LSP_DELTA1 = 0.2f; 1489 immutable float LSP_DELTA2 = 0.05f; 1490 bool SIGN_CHANGE(float a, float b) { return a*b < 0; } 1491 1492 float X2ANGLE(float x) { 1493 import std.math : acos; 1494 return acos(x); 1495 } 1496 1497 float cheb_poly_eva(float *coef, float x, int m) 1498 { 1499 int k; 1500 float b0, b1, tmp; 1501 1502 /* Initial conditions */ 1503 b0=0; /* b_(m+1) */ 1504 b1=0; /* b_(m+2) */ 1505 x*=2; 1506 /* Calculate the b_(k) */ 1507 for (k=m;k>0;k--) { 1508 tmp=b0; /* tmp holds the previous value of b0 */ 1509 b0=x*b0-b1+coef[m-k]; /* b0 holds its new value based on b0 and b1 */ 1510 b1=tmp; /* b1 holds the previous value of b0 */ 1511 } 1512 return(-b1 + 0.5f * x*b0+coef[m]); 1513 } 1514 1515 1516 /* float *a lpc coefficients */ 1517 /* int lpcrdr order of LPC coefficients (10) */ 1518 /* float *freq LSP frequencies in the x domain */ 1519 /* int nb number of sub-intervals (4) */ 1520 /* float delta grid spacing interval (0.02) */ 1521 int lpc_to_lsp (const float *a, int lpcrdr, float *freq, int nb, float delta) 1522 { 1523 import std.math : fabs; 1524 float temp_xr,xl,xr,xm=0; 1525 float psuml,psumr,psumm,temp_psumr /*,temp_qsumr*/; 1526 int i,j,m,flag,k; 1527 float[64] Q; /* ptrs for memory allocation */ 1528 float[64] P; 1529 float *Q16 = null; /* ptrs for memory allocation */ 1530 float *P16 = null; 1531 float *px; /* ptrs of respective P'(z) & Q'(z) */ 1532 float *qx; 1533 float *p; 1534 float *q; 1535 float *pt; /* ptr used for cheb_poly_eval() 1536 whether P' or Q' */ 1537 int roots=0; /* DR 8/2/94: number of roots found */ 1538 flag = 1; /* program is searching for a root when, 1539 1 else has found one */ 1540 m = lpcrdr/2; /* order of P'(z) & Q'(z) polynomials */ 1541 1542 /* Allocate memory space for polynomials */ 1543 Q[0..$] = 0; 1544 P[0..$] = 0; 1545 //Q = (FLOAT_DMEM*) calloc(1,sizeof(FLOAT_DMEM)*(m+1)); 1546 //P = (FLOAT_DMEM*) calloc(1,sizeof(FLOAT_DMEM)*(m+1)); 1547 1548 /* determine P'(z)'s and Q'(z)'s coefficients where 1549 P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */ 1550 1551 px = P.ptr; /* initialise ptrs */ 1552 qx = Q.ptr; 1553 p = px; 1554 q = qx; 1555 1556 *px++ = LPC_SCALING; 1557 *qx++ = LPC_SCALING; 1558 for(i=0;i<m;i++) { 1559 *px++ = (a[i]+a[lpcrdr-1-i]) - *p++; 1560 *qx++ = (a[i]-a[lpcrdr-1-i]) + *q++; 1561 } 1562 px = P.ptr; 1563 qx = Q.ptr; 1564 for(i=0;i<m;i++) { 1565 *px = 2**px; 1566 *qx = 2**qx; 1567 px++; 1568 qx++; 1569 } 1570 1571 px = P.ptr; /* re-initialise ptrs */ 1572 qx = Q.ptr; 1573 /* now that we have computed P and Q convert to 16 bits to 1574 speed up cheb_poly_eval */ 1575 P16 = P.ptr; 1576 Q16 = Q.ptr; 1577 1578 /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z). 1579 Keep alternating between the two polynomials as each zero is found */ 1580 1581 xr = 0; /* initialise xr to zero */ 1582 xl = FREQ_SCALE; /* start at point xl = 1 */ 1583 1584 for(j=0;j<lpcrdr;j++){ 1585 if(j&1) /* determines whether P' or Q' is eval. */ 1586 pt = Q16; 1587 else 1588 pt = P16; 1589 1590 psuml = cheb_poly_eva(pt,xl,m); /* evals poly. at xl */ 1591 flag = 1; 1592 while(flag && (xr >= -FREQ_SCALE)){ 1593 float dd; 1594 /* Modified by JMV to provide smaller steps around x=+-1 */ 1595 1596 dd=delta*(1.0f - 0.9f *xl*xl); 1597 if (fabs(psuml)<.2) 1598 dd *= 0.5f; 1599 1600 xr = xl - dd; /* interval spacing */ 1601 psumr = cheb_poly_eva(pt,xr,m);/* poly(xl-delta_x) */ 1602 temp_psumr = psumr; 1603 temp_xr = xr; 1604 1605 /* if no sign change increment xr and re-evaluate poly(xr). Repeat til 1606 sign change. 1607 if a sign change has occurred the interval is bisected and then 1608 checked again for a sign change which determines in which 1609 interval the zero lies in. 1610 If there is no sign change between poly(xm) and poly(xl) set interval 1611 between xm and xr else set interval between xl and xr and repeat till 1612 root is located within the specified limits */ 1613 1614 if(SIGN_CHANGE(psumr,psuml)) 1615 { 1616 roots++; 1617 1618 psumm=psuml; 1619 for(k=0;k<=nb;k++){ 1620 xm = 0.5f*(xl+xr); /* bisect the interval */ 1621 psumm=cheb_poly_eva(pt,xm,m); 1622 /*if(psumm*psuml>0.)*/ 1623 if(!SIGN_CHANGE(psumm,psuml)) 1624 { 1625 psuml=psumm; 1626 xl=xm; 1627 } else { 1628 psumr=psumm; 1629 xr=xm; 1630 } 1631 } 1632 1633 /* once zero is found, reset initial interval to xr */ 1634 if (xm>1.0) { xm = 1.0; } /* <- fixed a possible NAN issue here...? */ 1635 else if (xm<-1.0) { xm = -1.0; } 1636 freq[j] = X2ANGLE(xm); 1637 xl = xm; 1638 flag = 0; /* reset flag for next search */ 1639 } 1640 else{ 1641 psuml=temp_psumr; 1642 xl=temp_xr; 1643 } 1644 } 1645 } 1646 return(roots); 1647 } 1648 1649 1650 int lpcToLsp(const float[] src, float[] dst) // idxi=input field index 1651 { 1652 return lpcToLsp(src.ptr, dst.ptr, cast(int)src.length, cast(int)dst.length); 1653 } 1654 1655 int lpcToLsp(const float *src, float *dst, int Nsrc, int Ndst) // idxi=input field index 1656 { 1657 int nLpc = Nsrc; 1658 1659 /* LPC to LSPs (x-domain) transform */ 1660 int roots; 1661 roots = lpc_to_lsp (src, nLpc, dst, 10, LSP_DELTA1); 1662 if (roots!=nLpc) { 1663 roots = lpc_to_lsp (src, nLpc, dst, 10, LSP_DELTA2); // nLpc was Nsrc 1664 if (roots!=nLpc) { 1665 int i; 1666 for (i=roots; i<nLpc; i++) { 1667 dst[i]=0.0; 1668 } 1669 } 1670 } 1671 1672 return 1; 1673 } 1674 1675 //float X2ANGLE(float x) { 1676 // import std.math : acos; 1677 // return (acos(.00006103515625*(x)) * LSP_SCALING); 1678 //} 1679 float ANGLE2X(float a) { 1680 import std.math; 1681 return (cos(a)); 1682 } 1683 1684 void lspToLpc(const float[] src, float[] dst) { 1685 int lpcrdr = cast(int)src.length; 1686 lsp_to_lpc(src.ptr, dst.ptr, lpcrdr); 1687 } 1688 1689 void lsp_to_lpc(const float *freq, float *ak, int lpcrdr) 1690 /* float *freq array of LSP frequencies in the x domain */ 1691 /* float *ak array of LPC coefficients */ 1692 /* int lpcrdr order of LPC coefficients */ 1693 1694 1695 { 1696 int i,j; 1697 float xout1,xout2,xin1,xin2; 1698 float[] Wp; 1699 float * pw, n1, n2, n3, n4; 1700 float[] x_freq; 1701 int m = lpcrdr>>1; 1702 1703 Wp.length = 4*m+2; 1704 Wp[0..$] = 0; /* initialise contents of array */ 1705 1706 /* Set pointers up */ 1707 1708 pw = Wp.ptr; 1709 xin1 = 1.0; 1710 xin2 = 1.0; 1711 1712 x_freq.length = lpcrdr; 1713 for (i=0;i<lpcrdr;i++) 1714 x_freq[i] = ANGLE2X(freq[i]); 1715 1716 /* reconstruct P(z) and Q(z) by cascading second order 1717 polynomials in form 1 - 2xz(-1) +z(-2), where x is the 1718 LSP coefficient */ 1719 1720 for(j=0; j <= lpcrdr ;j++) { 1721 int i2=0; 1722 for(i=0; i<m; i++, i2 += 2) { 1723 n1 = pw + (i*4); 1724 n2 = n1 + 1; 1725 n3 = n2 + 1; 1726 n4 = n3 + 1; 1727 xout1 = xin1 - 2.0f*x_freq[i2] * *n1 + *n2; 1728 xout2 = xin2 - 2.0f*x_freq[i2+1] * *n3 + *n4; 1729 *n2 = *n1; 1730 *n4 = *n3; 1731 *n1 = xin1; 1732 *n3 = xin2; 1733 xin1 = xout1; 1734 xin2 = xout2; 1735 } 1736 xout1 = xin1 + *(n4 + 1); 1737 xout2 = xin2 - *(n4 + 2); 1738 if (j > 0) 1739 ak[j-1] = (xout1 + xout2)*0.5f; 1740 *(n4+1) = xin1; 1741 *(n4+2) = xin2; 1742 1743 xin1 = 0.0; 1744 xin2 = 0.0; 1745 } 1746 1747 } 1748 1749 1750 /* Autocorrelation in the time domain (for ACF LPC method) */ 1751 /* x is signal to correlate, n is number of samples in input buffer, 1752 *outp is array to hold #lag output coefficients, lag is # of output coefficients (= max. lag) */ 1753 void autoCorr(const float *x, const int n, float *outp, int lag) 1754 { 1755 int i; 1756 while (lag) { 1757 outp[--lag] = 0.0; 1758 for (i=lag; i < n; i++) { 1759 outp[lag] += x[i] * x[i - lag]; 1760 } 1761 } 1762 } 1763 1764 /* LPC analysis via acf (=implementation of Durbin recursion)*/ 1765 int calcLpcAcf(float * r, float *a, int p, float *gain, float *k) 1766 { 1767 int i, m = 0; 1768 float e; 1769 int errF = 1; 1770 float k_m; 1771 float[64] al; 1772 1773 if (!a) return 0; 1774 if (!r) return 0; 1775 1776 if ((r[0] == 0.0) || (r[0] == -0.0)) { 1777 for (i=0; i < p; i++) { 1778 a[i] = 0.0; 1779 } 1780 return 0; 1781 } 1782 1783 //al = (FLOAT_DMEM*)malloc(sizeof(FLOAT_DMEM)*(p)); 1784 1785 // Initialisation, Gl. 158 1786 e = r[0]; 1787 if (e==0.0) { 1788 for (i=0; i<=p; i++) { 1789 a[i] = 0.0; 1790 if (k) 1791 k[m] = 0.0; 1792 } 1793 } else { 1794 // Iterations: m = 1... p, Gl. 159 1795 for (m=1; m<=p; m++) { 1796 // Gl. 159 (a) 1797 float sum = 1.0f * r[m]; 1798 for (i=1; i<m; i++) { 1799 sum += a[i-1] * r[m-i]; 1800 } 1801 k_m = (-1.0f / e ) * sum; 1802 // copy refl. coeff. 1803 if (k) k[m-1] = k_m; 1804 // Gl. 159 (b) 1805 a[m-1] = k_m; 1806 for (i=1; i<=m/2; i++) { 1807 float x = a[i-1]; 1808 a[i-1] += k_m * a[m-i-1]; 1809 if ((i < (m/2))||((m&1)==1)) a[m-i-1] += k_m * x; 1810 } 1811 // update error 1812 e *= (1.0f - k_m*k_m); 1813 if (e==0.0) { 1814 for (i=m; i<=p; i++) { 1815 a[i] = 0.0; 1816 if (k) k[m] = 0.0; 1817 } 1818 break; 1819 } 1820 } 1821 } 1822 1823 if (gain) 1824 *gain=e; 1825 return 1; 1826 } 1827 1828 // return value: gain 1829 float calcLpc(const float *x, int Nsrc, float * lpc, int nCoeff, float *refl) 1830 { 1831 float gain = 0.0; 1832 float[64] acf; 1833 //if (acf == NULL) acf = (FLOAT_DMEM *)malloc(sizeof(FLOAT_DMEM)*(nCoeff+1)); 1834 autoCorr(x, Nsrc, acf.ptr, nCoeff + 1); 1835 calcLpcAcf(acf.ptr, lpc, nCoeff, &gain, refl); 1836 return gain; 1837 }