В процессе рихтовки вдруг ни с того, ни с сего на порядок (в 10 раз) снизилась скорость работы фрагмента программы
#include <malloc.h>
#include <math.h>
#include "MultiDraw.h"
#include "PrepCurve.h"
#include "IIR_LP-8_16k_1532k.h"
// Задача:
// Есть записанный на частоте дискретизации freqSou сэмпл ноты с частотой noteSou.
// Нужно пересчитать его для другой ноты noteDest для частоты дискретизации freqDest.
// Для этого нота сначала передискретизируется на высокую частоту freqP (равную 1.536 МГц),
// там пропускается через ФНЧ и оттуда передискретизируется на freqDest.
// Для передискретизации вычисляется, сколько периодов freqP помещается в один период freqSou :
// periodSou = freqP / freqSou * noteSou / noteDest;
// а дальше по отдельности вычисляются целая и дробная части этого периода.
// одним отсчетом частоты freqSou заполняется либо (int)periodSou отсчетов freqP, либо ((int)periodSou + 1) отсчет.
// Соотношение между этими случаями регулируется величиной дробной части периода, т.е. происходит своеобразный dithering,
// в результате которого средняя выходная частота получается точной, а период немного колеблется вокруг среднего значения.
// Сигнал с высокой частотой дискреизации пропускается через ФНЧ 8-го порядка и передискретизируется на нужную частоту
// дискретизации - 48 кГц. С последним проблем не возниукает, так как частоты кратны.
// Для тестов проводится переинтерполяция большого буфера (4 секунды звучания), но в процессе работы нужно, чтобы
// задержка сигнала была минимальной. Поэтому все операции нужно проводить на с буфером, а буквально по одному отсчету.
// И здесь у меня возникла проблема: после выделения алгоритма в отдельную функцию, скорость его работы уменьшилась раз в 10.
// На Intel c частотой ~2.5-3.0 ГГц время выполнения составляет порядка 18 и 180 мс соответственно (Visual Studio 2019)
// Хотелось бы выяснить причину этого явления и устранить ее.
//
// PS. Для сравнения с другими конфигурациями цикл ФНЧ выполняется примерно 700 мс.
class QPTime {
public:
QPTime() {
BOOL b = QueryPerformanceFrequency(&Frequency);
GQPTdividor = (Frequency.LowPart + 500) / 1000;
}
~QPTime() {}
void printConst() {
BOOL b = QueryPerformanceFrequency(&Frequency);
outf << "bool: " << b << ", Freq: " << Frequency.HighPart << " " << Frequency.LowPart << ", dividor: " << GQPTdividor << endl;
b = QueryPerformanceCounter(&time);
outf << "bool: " << b << ", time: " << time.HighPart << " " << time.LowPart << endl;
}
unsigned int GetQPtime() { // возвращает текущее время в мс. Вероятно, исходный счетчик работает в единицай 1/3117265 секунды.
LARGE_INTEGER time;
BOOL b = QueryPerformanceCounter(&time);
if (b) return time.LowPart / GQPTdividor; // / 3117;
return 0;
}
protected:
LARGE_INTEGER Frequency, time;
int GQPTdividor; // константа для вычисления времени в мс.
};
constexpr double Pi2 = 6.283185307179586476925286766559;
constexpr double freqDest = 48000.0; // конечная частота дискретизации
constexpr double freqP = 1536000.0; // помежуточная частота дискретизации
constexpr double dTime = 4.0; // время воспроизведения 4 секунды
constexpr int offset = 2; // -6; // 2; // 16; // 20; // дополнительное смещение в промежуточном массиве: оптимальный коэффициент 20-21
constexpr double knorm = 0.915; // 0.920; //0.915; // 0.912 - нормировочный коэффициент, откуда - непонятно. Так работает ФНЧ.
int size0; // рабочий размер исходного массива (выделяется с запасом)
int size1; // рабочий размер выходного массива (выделяется с запасом)
int sizeP; // рабочий размер промежуточного массива (выделяется с запасом)
double* buf0; // буфер для образца
double* buf1; // выходной буфер
double* bufP; // промежуточный буфер
double* buf2; // 2-й вариант выходного буфера
double* bufN; // то, что должно быть в идеале
double* x; // шкала времени freqDest
double* y; // шкала времени freqP
double* dif1; // разница, умноженная на 100
double* dif2; // разница, умноженная на 100
double* P2; // промежуточный массив после вычитания из него образца и умноженный на 100
void allocateMemory() {
buf0 = new(double[size0 + 200]);
buf1 = new(double[size1 + 200]);
bufP = new(double[sizeP + 100000]);
buf2 = new(double[size1 + 200]);
bufN = new(double[size1 + 200]);
x = new(double[size1 + 200]);
y = new(double[sizeP + 10000]);
dif1 = new(double[size1 + 200]);
dif2 = new(double[size1 + 200]);
P2 = new(double[sizeP + 100000]);
}
class fillP { // класс для выдачи очередного сэмпла массива 1536к, из массива buf0 с шагом (целая часть = stepSou, дробная часть = dstepSou)
public:
fillP() {}
~fillP() {}
void reset(double noteSou_I, double noteDest_I, double freqSou_I) { // частота образца ноты, частота нужной ноты, частота дискретизации образца
mi_idx0 = 0; // индекс во входном буфере buf0
mi_idxStep = 0; // счетчик до stepSou, количество "вставок" buf0[i - 1]
mi_mInsert = 0; // количество "вставок" buf0[i - 1] в bufP[j] (нужен только для отладки)
md_fracSumDStep = 0.0; // дробная часть суммы шагов
mi_j = 0; // вспомогательный счетчик цикла а bufP для проверки
double periodSou = freqP / freqSou_I * noteSou_I / noteDest_I; // период действия одного сэмпла в тактах 1536к : (noteSou*freqP)/(noteDest*freqSou)
mi_stepSou = (int)periodSou; // целая часть периода
md_dstepSou = periodSou - mi_stepSou; // дробная часть периода
}
void printM() { outf << " m: " << mi_mInsert << " of " << size0 << ", step0: " << mi_stepSou << ", dtstep: " << md_dstepSou << endl; }
inline double next() { // выдает следующий сэмпл на частоте freqP без фильтрации - основной исследуемый вариант, без массива на высокой частоте (вариант 3)
if (mi_idxStep++ <mi_stepSou) { // результат текущего шага в пределах от 0-1 до <stepSou
return buf0[mi_idx0];
} else {
md_fracSumDStep += md_dstepSou;
mi_idx0++;
if (md_fracSumDStep >= 1.0) { // дополнительный и последний результат текущего шага
md_fracSumDStep -=1.0;
mi_mInsert++;
mi_idxStep = 0;
return buf0[mi_idx0 - 1];
} else { // это по сути уже первый результат следующего шага
mi_idxStep = 1;
return buf0[mi_idx0];
}
}
}
inline void next1() { // осущесьвляет переинтерполяцию через массивы - вариант для сравнения с вложенным циклом (вариант 4)
for (mi_idx0 = 0; mi_idx0 < size0; mi_idx0++) {
for (int mi_idxStep = 0; mi_idxStep < mi_stepSou; mi_idxStep++) {
bufP[mi_j++] = buf0[mi_idx0]; // в промежуточном массиве отсчеты сдвинуты в среднем на (periodSou/2.0 - 0.5) отсчетов ~= 40
}
md_fracSumDStep += md_dstepSou;
if (md_fracSumDStep >= 1.0) {
bufP[mi_j++] = buf0[mi_idx0];
md_fracSumDStep -= 1.0;
mi_mInsert++;
}
}
}
inline void next2() { // осущесьвляет переинтерполяцию через массивы - вариант для сравнения с одним циклом (вариант 5)
while (mi_j < sizeP) {
if (mi_idxStep++ < mi_stepSou) {
bufP[mi_j++] = buf0[mi_idx0];
}
else {
md_fracSumDStep += md_dstepSou;
if (md_fracSumDStep >= 1.0) {
bufP[mi_j++] = buf0[mi_idx0];
md_fracSumDStep -= 1.0;
mi_mInsert++;
}
mi_idxStep = 0;
mi_idx0++;
}
}
}
inline void next3() { // работает через вариант 3, но внутри класса (вариант 6)
for (mi_j = 0; mi_j < sizeP; mi_j++) {
bufP[mi_j] = next();
}
}
protected:
int mi_idx0; // индекс во входном буфере buf0
int mi_idxStep; // счетчик до stepSou, количество "вставок" buf0[i - 1] (нужен только для отладки)
int mi_mInsert; // количество "вставок" buf0[i - 1] в bufP[j] (нужен только для отладки)
double md_fracSumDStep; // дробная часть суммы шагов
int mi_j; // вспомогательный счетчик цикла а bufP для проверки
int mi_stepSou; // целая часть периода действия входного сэмпла (столько раз он повторяется на промежуточной частоте)
double md_dstepSou; // дробная часть периода действия входного сэмпла (столько раз он повторяется на промежуточной частоте)
};
double noteSou = 1000.0; // нота (частота) образца на частоте дискретизации freqSou
double noteDest = 834.78; // 1333.33; //834.780; ////495.0; // 437.0; // выходная (требуемая) нота - это noteSou на эффектиыной частоте дискретизации freqSou*noteDest/noteSou
double freqSou = 44100.0; // начальная частота дискретизации
double periodDest = freqP / freqDest; // период действия одного сэмпла в тактах 1536к
double dT1 = dTime * freqDest; // полный размер буфера для выходного сигнала
double dTP = dTime * freqP; // полный размер буфера для промежуточного сигнала
double dT0 = dTime * freqSou * (noteDest / noteSou); // полный размер буфера для исходного сигнала
double periodSou = freqP / freqSou * noteSou / noteDest; // период действия одного сэмпла в тактах 1536к // (noteSou*freqP)/(noteDest*freqSou)
int stepSou = (int)periodSou;
double dstepSou = periodSou - stepSou;
QPTime QPT;
fillP FP;
MultiDraw MD;
MultiDrawL MDL;
void prepCurves() {
// сначала настраиваем систему контроля времени
QPT.printConst();
// а теперь переходим к вычислениям
outf << "per0: " << periodSou << ", per1: " << periodDest << endl;
int stepDest = (int)(freqP / freqDest);
outf << "stepSou: " << stepSou << ", dstepSou: " << dstepSou << ", stepDest: " << stepDest << endl;
size0 = (int)dT0;
size1 = (int)dT1;
sizeP = (int)dTP;
outf << "size0: " << size0 << ", size1: " << size1 << ", sizeP: " << sizeP << endl;
allocateMemory();
int t0 = QPT.GetQPtime();
// подготовка исходного массива (если он есть готовый, например, существующий сэмпл - не нужно)
for (int i = 0; i < size0; i++) {
buf0[i] = sin(Pi2 * i * (noteSou / freqSou));
}
int t1 = QPT.GetQPtime();
// очистка промежуточного массива и заполнение его данными
for (int i = 0; i < sizeP; i++) {
bufP[i] = 0.0;
}
int t2 = QPT.GetQPtime();
double dj = 0;
int m = 0;
int j = 0; // (1)
int currVariant = 0;
#define _variant_p_6 // здесь мы указываем вариант, который бужем исследовать: от _variant_p_1 до _variant_p_6
#ifdef _variant_p_1
currVariant = 1; // двойной цикл заполнения массива bufP (для сравнения)
for (int i = 0; i < size0; i++) {
for (int k = 0; k < stepSou; k++) {
bufP[j++] = buf0[i]; // в промежуточном массиве отсчеты сдвинуты в среднем на (periodSou/2.0 - 0.5) отсчетов ~= 40
}
dj += dstepSou;
if (dj >= 1.0) {
bufP[j++] = buf0[i];
dj -= 1.0;
m++;
}
} // (1)
#endif
#ifdef _variant_p_2
currVariant = 2; // одиночный цикл заполнения массива bufP (для сравнения)
int k = 0; // (2)
int i = 0;
j = 0;
while (j < sizeP) {
if (k++ < stepSou) {
bufP[j++] = buf0[i];
} else {
dj += dstepSou;
if (dj >= 1.0) {
bufP[j++] = buf0[i];
dj -= 1.0;
m++;
}
k = 0;
i++;
}
} // (2)
#endif
FP.reset(noteSou, noteDest, freqSou);
#ifdef _variant_p_3
currVariant = 3; // это должно быть основным вариантом: переинтерполяцмя bufSou => bufP на лету
for (j = 0; j < sizeP; j++){
bufP[j] = FP.next();
}
#endif
#ifdef _variant_p_4
currVariant = 4; // метод класса с заполнением массива двойным циклом (для сравнения)
FP.next1();
#endif
#ifdef _variant_p_5
currVariant = 5; // метод класса с заполнением массива одиночным циклом (для сравнения)
FP.next2();
#endif
#ifdef _variant_p_6
currVariant = 6; // метод класса с заполнением массива, используя основной метод класса - заполнение на лету (для сравнения)
FP.next3();
#endif
FP.printM();
int t3 = QPT.GetQPtime();
// фильтруем ФНЧ 8 порядка на частоте 16000 гц промежуточный массив, оцифрованный на частоте 1536000 Гц
for (int i = 0; i < sizeP; i++) {
bufP[i] = iir(bufP[i]); // фильтр тоже дает какую-то сдвижку - предположительно 4 отсчета (для 8 порядка фильтра)
}
int t4 = QPT.GetQPtime();
// вычисляем два варианта выходного массаива, нормируем их ?!? - непонятноЮ откуда коэффициент 0.912, но - какой есть (перенести offset в начальное значение j)
j = offset;
for (int i = 0; i < size1; i++) {
buf2[i] = bufP[j++]; // здесь - 20 отсчетов
for (int k = 1; k < stepDest; k++) {
buf2[i] += bufP[j++]; // еще на 16 отсчетов, всего 16+20 = 36 отсчета (назад)
}
buf2[i] *= (knorm/stepDest);
}
int t5 = QPT.GetQPtime();
// готовим массивы Х-координаты (нужно только для графика), а также готовим образец массива-эталона
for (int i = 0; i < size1; i++) {
x[i] = Pi2 * i * (noteDest / freqDest);
bufN[i] = sin(Pi2 * i * (noteDest / freqDest));
for (int k = 0; k < stepDest; k++) {
y[i * stepDest + k] = Pi2 * (i + ((double)k/stepDest)) * (noteDest / freqDest);
}
}
for (int i = 0; i < 30; i++) {
outf << i << " \t" << x[i] << " \t" << buf1[i] << " \t" << buf2[i] << " \t" << bufN[i] << endl;
}
int maxl = 120;
// вычисляем массивы разности с эталоном - для анализа
for (int i = 0; i < size1; i++) {
// dif1[i] = (buf1[i + 3] - bufN[i]) * 100.0; // к имеющимся добавляем еще 96 отсчетов
dif2[i] = (buf2[i + 3] - bufN[i]) * 100.0; // итого: при заполнении bufP: 40 отсчетов, buf1 и buf2 - еще -36, фильтр еще +4, итого: +8
}
// вычисляем разницу с эталоном на пормежуточной частоте 1536000 Гц
for (int i = 0; i < sizeP; i++) {
P2[i] = (bufP[i + 96+16 + offset]*knorm - sin(Pi2 * i * (noteDest / freqP)))*100.0;
}
int addi = 3; // +110/*-4*/; // похоже, результат примерно на 1% (4 точки на ~440 осчетв выше по частоте)
// посылаем графики на построение
// MDL.Add(maxl, x, &buf1[addi+115]); // size1 черный
MDL.Add(maxl, x, &buf2[addi+115]); // size1 синий
MDL.Add(maxl, x, bufN); // size1 зеленый
MDL.Add(maxl *periodDest, y, &bufP[addi*stepDest + 16 + offset]); // sizeP голубой // сдвиг на 96 + 3520 + 36
// MDL.Add(maxl, x, dif1); // size1 зеленый
MDL.Add(maxl, x, &dif2[110]); // size1 зеленый
MDL.Add(maxl * periodDest, y, &P2[(addi+110) * stepDest + 16 + offset - 96 - 16 ]); // sizeP голубой
outf << "\n (in 1 ms unit) sin(19k): " << t1 - t0 << ", zero P(1M5): " << t2 - t1 << ", fill P(1M5): " << t3 - t2 << ", filter(1M5): " << t4-t3 << ", buf1/buf2(48r): " << t5-t4 << ", Variant: " << currVariant << '\n' << endl;
}
Отрисовка волновой формы и искажений, а также фильтр 8-го порядка - в других файлах. Думаю, для вопроса они не существенны. По сути интересует только фрагмент в строках 102-118. Почему он работает на порядок медленнее аналогичных фрагментов, расположенных рядом?