Один из красивых и интересных примеров применения CUDA представлен здесь: реализация очень точного кроссовера (полосового аудио-фильтра) для домашней акустики на маломощной видеокарте. Человек использует фильтр порядка 8192! Такие порядки нечасто встретишь в подобных алгоритмах :)

КИХ-фильтр — классический пример map/reduce алгоритма: проходим по всему массиву предыдущих значений входного параметра, умножая их на коэффициенты фильтра, потом редуцируем массив, складывая все полученные произведения. Как и любой map-алгоритм, он легко поддаётся параллелизации.

Попробуем повторить пример на сайте NVidia, про кроссовер с фильтрами на 4 полосы аудиочастот. Каждая полоса будет вычленяться своим КИХ-фильтром, а границы полос (т.е. боковые полосы) постараемся сделать как можно тише. Для этого понадобятся фильтры высокого порядка.

Теория

Для удобства разобьём весь алгоритм на этапы:

  1. Вычисление массива коэффициентов фильтра;
  2. Получение новых данных и складывание их в память;
  3. Map — умножение массива данных на массив коэффициентов;
  4. Reduce — сложение всех элементов полученного массива;
  5. Возврат полученного значения.

Коэффициенты фильтра достаточно вычислить один раз в начале, а остальные действия необходимо производить в цикле, с каждым новым входным значением.

Вычисление коэффициентов КИХ-фильтра

Не буду расписывать эту процедуру слишком подробно, благо она достаточно освещена во многих источниках. Вкратце — это просто Фурье-преобразование желаемой АЧХ.

Расчёт коэффициентов достаточно сделать один раз перед началом цикла преобразований. Если хочется управлять кроссовером как эквалайзером — коэффициенты можно периодически пересчитывать, достаточно будет не чаще 10 раз в секунду.

Получение и складирование новых данных

Поскольку мы имеем дело со звуком, новые значения поступают как минимум 44100 раз в секунду. Пересылать эти байты по отдельности было бы весьма накладно, однако операционная система и так объединяет звуковые данные в чанки — кусочки по 1024 элемента, к примеру. Таким образом, 43 раза в секунду мы будем пересылать в видеопамять блоки данных размером 1024*2(стерео)*2(2 байта = 16 бит АЦП) = 4кБ.

Далее, принцип действия КИХ-фильтра — это работа с последними N значениями входного параметра, т.е. скользящее окно. Очень удобно применить структуру данных «кольцевой массив», в которой можно использовать изящную адресацию по модулю. Так, каждое новое значение будем складывать в этот массив на инкрементирующийся адрес, а когда дойдём до конца массива — автоматически перейдём в начало. Забирать значения будем точно так же, модульной адресацией.

С учётом вышесказанного про пакеты данных — при поступлении новых данных будем складировать их сразу кучей, для этого потребуется две операции memcpy — до достижения конца массива, и для копирования остатка на начало массива. Конечно, если новые данные не попадают на границу — можно обойтись одним memcpy.

Map — умножение данных на коэффициенты

Здесь начинается магия CUDA. Вычислить полностью весь массив произведений мы сможем за крайне маленькое время. Нужно на массив входящих значений наложить массив коэффициентов фильтра, и попарно их перемножить. Результат сложим в переменную gpu_out. Функция очень проста, лишь небольшое затруднение может встретиться при расчёте адресации в массивах.

Reduce — сложение всех элементов массива

Строго говоря, эта операция не параллелится напрямую — потому, что здесь происходит модификация одной и той же переменной (аккумулятора). Выйти из положения можно двумя путями — атомарными операциями, либо разбиением всего множества сумм на независимые кусочки. Атомарные операции потребуют довольно большого времени на выполнение, а разбиение даст возможность  попарно сложить всё за O(log N) времени. И всё равно это будет гораздо быстрее, чем пересылка всего массива в ОЗУ и сложение силами CPU.

Более того, только на первых этапах сложения заняты все (или многие) потоки, а с каждым новым этапом количество занятых потоков уменьшается в два раза. Здесь кроется ещё одна возможность ускорения: если нам нужно вычислять одновременно несколько КИХ-фильтров (а нам действительно нужно это делать) — последние их этапы можно выполнять параллельно в одном блоке. Таким образом, время редукции 4 КИХ-фильтров будет ненамного больше времени редукции 1 фильтра.

Поскольку время, занимаемое операцией reduce, растёт как логарифм числа taps — а длительность операции map не зависит от количества taps, становится очевидным что количество taps можно увеличивать очень сильно. К примеру, порядок 256 потребует 8 проходов функции сложения, а порядок 8192 (как у автора) — 13. Фильтр становится в 32 раза точнее — а время выполнения reduce вырастает всего на 60%. Время полного прохода так и  вообще вряд ли увеличится больше чем на 20%.

Возврат полученного значения

Возвращать будем тоже чанк, набор из 1024*2 значений, поскольку его удобно пустить сразу в звуковую карту.

Первая версия

Для начала сделаем попроще, лишь бы работало.

__global__ void add(float *a, const int step)
{
	int i = threadIdx.x*2*step;
	a[i]=a[i]+a[i+step];
}

__global__ void mul(const float *in, const float *coeffs, float *out, int count, int pos)
{
	int i = threadIdx.x;
	int j = i-pos-1; if(j<0) j+=count;
	out[i] = in[i] * coeffs[j];
}

int main(void)
{
	#define Ntap 256
	int size = sizeof(float) * Ntap;

	gettimeus();

	float *cpu_in = (float *)malloc(size);		float *gpu_in;		cudaMalloc((void**)&gpu_in, size);
	float *cpu_coeffs = (float *)malloc(size);	float *gpu_coeffs;	cudaMalloc((void**)&gpu_coeffs, size);
	float *cpu_out = (float *)malloc(size);		float *gpu_out;		cudaMalloc((void**)&gpu_out, size);
	float *cpu_depot = (float *)malloc(size);	float *gpu_depot;	cudaMalloc((void**)&gpu_depot, size);

	float new_value=0;

	float FIRCoef[200] = { -0.00005122648355433649, -0.00003699380424890057, -0.00002103274900807157, -0.00000327211239829301, 0.00001635979561118467, 0.00003793506790277917, 0.00006152611486757026, 0.00008720557545659427, 0.00011504622460606516, 0.00014512087703298453, 0.00017750228738663842, 0.00021226304675927511, 0.00024947547555152126, 0.00028921151270481544, 0.00033154260131288549, 0.00037653957062155202, 0.00042427251443229955, 0.00047481066591696940, 0.00052822226886740226, 0.00058457444540803622, 0.00064393306021746997, 0.00070636258130396556, 0.00077192593736691978, 0.00084068437176235220, 0.00091269729308755875, 0.00098802212241961857, 0.00106671413728598260, 0.00114882631247167110, 0.00123440915774948940, 0.00132351055255903150, 0.00141617557759524330, 0.00151244634325898670, 0.00161236181502490700, 0.00171595763591740680, 0.00182326594635640960, 0.00193431520152507880, 0.00204912998614538430, 0.00216773082631934030, 0.00229013399815036590, 0.00241635133327626700, 0.00254639002195676010, 0.00268025241450550950, 0.00281793582128789980, 0.00295943231045147840, 0.00310472850184871420, 0.00325380535610515430, 0.00340663795955351840, 0.00356319530760445960, 0.00372344008919155800, 0.00388732847226631010, 0.00405480988627247420, 0.00422582679549942000, 0.00440031446019303710, 0.00457820068969120190, 0.00475940559824460440, 0.00494384137256164390, 0.00513141204811425580, 0.00532201327585483250, 0.00551553205608981690, 0.00571184643193156650, 0.00591082516523284280, 0.00611232743945058600, 0.00631620261992944100, 0.00652229004873280250, 0.00673041879255496630, 0.00694040725526577120, 0.00715206264461221670, 0.00736518040938842840, 0.00757954383436861520, 0.00779492389319432110, 0.00801107921961348270, 0.00822775582904587170, 0.00844468624500578810, 0.00866158806964697350, 0.00887816258501568710, 0.00909409420605807040, 0.00930905111071003360, 0.00952268624522066740, 0.00973463694634357430, 0.00994452172002656570, 0.01015193468863220300, 0.01035644080784727000, 0.01055757583582475800, 0.01075485223319555200, 0.01094776627472374100, 0.01113579679964778900, 0.01131838792838221400, 0.01149491924755947600, 0.01166468213598349800, 0.01182688707507986400, 0.01198071078821300400, 0.01212535305940168200, 0.01226003171318920500, 0.01238383822962717400, 0.01249543932773643600, 0.01259274209407373700, 0.01267279770133641400, 0.01273234773794170200, 0.01277006407712861500, 0.01279300545121231500, 0.01277006407712861500, 0.01273234773794170200, 0.01267279770133641400, 0.01259274209407373700, 0.01249543932773643600, 0.01238383822962717400, 0.01226003171318920500, 0.01212535305940168200, 0.01198071078821300400, 0.01182688707507986400, 0.01166468213598349800, 0.01149491924755947600, 0.01131838792838221400, 0.01113579679964778900, 0.01094776627472374100, 0.01075485223319555200, 0.01055757583582475800, 0.01035644080784727000, 0.01015193468863220300, 0.00994452172002656570, 0.00973463694634357430, 0.00952268624522066740, 0.00930905111071003360, 0.00909409420605807040, 0.00887816258501568710, 0.00866158806964697350, 0.00844468624500578810, 0.00822775582904587170, 0.00801107921961348270, 0.00779492389319432110, 0.00757954383436861520, 0.00736518040938842840, 0.00715206264461221670, 0.00694040725526577120, 0.00673041879255496630, 0.00652229004873280250, 0.00631620261992944100, 0.00611232743945058600, 0.00591082516523284280, 0.00571184643193156650, 0.00551553205608981690, 0.00532201327585483250, 0.00513141204811425580, 0.00494384137256164390, 0.00475940559824460440, 0.00457820068969120190, 0.00440031446019303710, 0.00422582679549942000, 0.00405480988627247420, 0.00388732847226631010, 0.00372344008919155800, 0.00356319530760445960, 0.00340663795955351840, 0.00325380535610515430, 0.00310472850184871420, 0.00295943231045147840, 0.00281793582128789980, 0.00268025241450550950, 0.00254639002195676010, 0.00241635133327626700, 0.00229013399815036590, 0.00216773082631934030, 0.00204912998614538430, 0.00193431520152507880, 0.00182326594635640960, 0.00171595763591740680, 0.00161236181502490700, 0.00151244634325898670, 0.00141617557759524330, 0.00132351055255903150, 0.00123440915774948940, 0.00114882631247167110, 0.00106671413728598260, 0.00098802212241961857, 0.00091269729308755875, 0.00084068437176235220, 0.00077192593736691978, 0.00070636258130396556, 0.00064393306021746997, 0.00058457444540803622, 0.00052822226886740226, 0.00047481066591696940, 0.00042427251443229955, 0.00037653957062155202, 0.00033154260131288549, 0.00028921151270481544, 0.00024947547555152126, 0.00021226304675927511, 0.00017750228738663842, 0.00014512087703298453, 0.00011504622460606516, 0.00008720557545659427, 0.00006152611486757026, 0.00003793506790277917, 0.00001635979561118467, -0.00000327211239829301, -0.00002103274900807157, -0.00003699380424890057, -0.00005122648355433649, -0.00006380142958806893 };

	for(int i=0; i<Ntap; i++) cpu_in[i]=0;
	for(int i=0; i<Ntap; i++) cpu_coeffs[i]=0;
	for(int i=0; i<200; i++) cpu_coeffs[i]=FIRCoef[i];

	cudaMemcpy(gpu_in,		cpu_in,		size,					cudaMemcpyHostToDevice);
	cudaMemcpy(gpu_coeffs,  cpu_coeffs, sizeof(float) * Ntap,	cudaMemcpyHostToDevice);
	#define pi 3.1415926535

	//
	// Main cycle
	//

	int i, j, step;
	int freq=1040, discr=44100;
	float energy=0, val;
	long long time=-gettimeus();
	long long time2, time_mem_to, time_mem_from, time_calc_mul, time_calc_add;	
	int N=16;//Ntap;
	for(j=0; j<discr; j++) {
		time2=-gettimeus();
		new_value=sin(freq*j*2.0*pi/discr);
		val=new_value;
		i=j%Ntap;

		time_mem_to=-gettimeus();
		cudaMemcpy(gpu_in+i, &new_value, sizeof(float), cudaMemcpyHostToDevice);		//push new value
		time_mem_to+=gettimeus();

		time_calc_mul=-gettimeus();
		mul<<<1, N>>>(gpu_in, gpu_coeffs, gpu_out, N, i);							//calc multiplies, O(1)
		time_calc_mul+=gettimeus();

		time_calc_add=-gettimeus();
		step=1; while(step<N) { add<<<1, N/step/2>>>(gpu_out, step); step*=2; }	//calc sums, O(log n)
		time_calc_add+=gettimeus();

		time_mem_from=-gettimeus();
		cudaMemcpy(&new_value, gpu_out, sizeof(float), cudaMemcpyDeviceToHost);			//pop new value
		time_mem_from+=gettimeus();
		time2+=gettimeus();

		printf("%f %fn", val, new_value);
		energy+=fabs(new_value);
	}
	time+=gettimeus();
	printf("%f %lld s; all: %lld, time_mem_to: %lld, time_calc_mul: %lld, time_calc_add: %lld, time_mem_from: %lldn", energy, time/1000000, time2, time_mem_to, time_calc_mul, time_calc_add, time_mem_from);

	free(cpu_in);		cudaFree(gpu_in);
	free(cpu_coeffs);	cudaFree(gpu_coeffs);
	free(cpu_out);		cudaFree(gpu_out);
	scanf("");
}

О чём здесь речь? Передаём в ГПУ таблицу коэффициентов фильтра, а массив входящих значений зануляем. Далее в цикле отправляем новый отчёт, и проводим вычисления — map с умножением на фильтр и reduce со сложением. Получаем обратно вычисленное значение, переходим к следующей итерации цикла.

В качестве контрольного значения я вычисляю «энергию» звукового сигнала (сумма всех значений по модулю), и тестирую алгоритм на чистом тоне 1040Гц. Результаты теста, прямо говоря, удручают — 44100 отчётов, т.е. 1 секунда моно-записи, обрабатывается в течение 21 секунды. До реал-тайма здесь далеко :)

Однако, я включил в цикл простенький профайлер, и что он нам говорит? «all: 463, time_mem_to: 15, time_calc_mul: 12, time_calc_add: 90, time_mem_from: 402«. Огромное время уходит на копирование значения из видеопамяти в ОЗУ, на 4 байта — 400мкс, это 10 кБ/с. Жуть какая-то, наши прошлые измерения производительности memcpy говорили о 4ГБ/с!

Это говорит об одном — гораздо выгоднее передавать большие объёмы данных. В нашем случае можно передавать сразу целый чанк, например целую секунду — 44100 отчёта. Заодно, сделаем это и для копирования в ГПУ.

Вторая версия

__global__ void add(float *a, const int step, const int tr)
{
    int i = threadIdx.x*tr;
	a[i]=a[i]+a[i+step];
}

__global__ void mul(const float *in, const float *coeffs, float *out, int count, int pos)
{
	int i = threadIdx.x;
	int j = i-pos-1; if(j<0) j+=count;
	out[i] = in[i] * coeffs[j];
}

int main(void)
{
	#define Ntap 256
    int size = sizeof(float) * Ntap;
	int d_size = sizeof(float) * 44100;

	gettimeus();

	float *cpu_in = (float *)malloc(size);		float *gpu_in;		cudaMalloc((void**)&gpu_in, size);
    float *cpu_coeffs = (float *)malloc(size);	float *gpu_coeffs;	cudaMalloc((void**)&gpu_coeffs, size);
    float *cpu_out = (float *)malloc(size);		float *gpu_out;		cudaMalloc((void**)&gpu_out, size);
    float *cpu_depot = (float *)malloc(d_size);	float *gpu_depot;	cudaMalloc((void**)&gpu_depot, d_size);
	float new_value=0;

	float FIRCoef[200] = ; //Возьмите в первом листинге

    for(int i=0; i<Ntap; i++) cpu_in[i]=0;
    for(int i=0; i<Ntap; i++) cpu_coeffs[i]=0;
	for(int i=0; i<200; i++) cpu_coeffs[i]=FIRCoef[i];

    cudaMemcpy(gpu_in,		cpu_in,		size,					cudaMemcpyHostToDevice);
    cudaMemcpy(gpu_coeffs,  cpu_coeffs, sizeof(float) * Ntap,	cudaMemcpyHostToDevice);
	#define pi 3.1415926535

	//
	// Main cycle
	//	
	int i, j, step;
	int freq=1040, discr=44100;
	float energy=0, val;
	long long time=-gettimeus();
	long long time2, time_mem_to, time_mem_from, time_calc_mul, time_calc_add;	
	int N=Ntap;

	for(j=0; j<discr; j++) cpu_depot[j]=sin(freq*j*2.0*pi/discr);
    cudaMemcpy(gpu_depot, cpu_depot, d_size, cudaMemcpyHostToDevice);
	for(j=0; j<discr; j++) {
		time2=-gettimeus();
		i=j%Ntap;
		time_mem_to=-gettimeus();
		cudaMemcpy(gpu_in+i, gpu_depot+j, sizeof(float), cudaMemcpyDeviceToDevice);	//push new value
		time_mem_to+=gettimeus();

		time_calc_mul=-gettimeus();
		mul<<<1, N>>>(gpu_in, gpu_coeffs, gpu_out, N, i);							//calc multiplies, O(1)
		time_calc_mul+=gettimeus();

		time_calc_add=-gettimeus();
		step=1; while(step<N) { add<<<1, N/step/2>>>(gpu_out, step); step*=2; }		//calc sums, O(log n)
		time_calc_add+=gettimeus();

		time_mem_from=-gettimeus();
		cudaMemcpy(gpu_depot+j, gpu_out, sizeof(float), cudaMemcpyDeviceToDevice);	//push new value
		time_mem_from+=gettimeus();
		time2+=gettimeus();
	}
    cudaMemcpy(cpu_depot, gpu_depot, d_size, cudaMemcpyDeviceToHost);
	for(j=0; j<discr; j++) energy+=fabs(cpu_depot[j]);  //184.95

	time+=gettimeus();
	printf("%f %f s; all: %lld, time_mem_to: %lld, time_calc_mul: %lld, time_calc_add: %lld, time_mem_from: %lldn", energy, time/1000000.0, time2, time_mem_to, time_calc_mul, time_calc_add, time_mem_from);	
	//
	// End of main cycle
	//
    free(cpu_in);		cudaFree(gpu_in);
    free(cpu_coeffs);	cudaFree(gpu_coeffs);
    free(cpu_out);		cudaFree(gpu_out);
    scanf("");
}

Тут уже гораздо веселее, запись результата занимает всего 10мкс. Как это сделано? Заводим массив gpu_depot, который передаёт значения в обработку, и принимает результат. Во время обработки мы копируем значение из этого массива в gpu_in. Поскольку оба массива находятся в видеопамяти, такое копирование занимает крайне малое время.

После обработки копируем результат обратно в gpu_depot за те же 10мкс, и лишь после того как весь чанк обработан — передаём его в ОЗУ. Теперь один отчёт обрабатывается всего за 47мкс, а весь секундный чанк — за 2.3 секунды. Это всё ещё совсем не реалтайм, нужно ещё улучшать алгоритм.