(Prog) Расчет FIR-фильтра через ДПФ над конечным полем

Описание

Расчет FIR-фильтра (КИХ-фильтра) с использованием теоретико-числового преобразования (ДПФ над конечным полем) и системы остаточных классов. Для склейки ответа используется модулярный аналог Overlap-Add.

Исходник на C

#define _USE_MATH_DEFINES

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

#define MAX_MOD_NUM 10
#define NL 1000
#define TRUE 1
#define FALSE 0
#define MAXPRIME 1000000

int prime[MAXPRIME];
int plist[MAXPRIME];
int pnum;

int find_bit(int p) {
	int i, k;
	k = 1;
	for (i = 0; i < 32; i++) {
		if (k > p)
			return i;
		k *= 2;
	}
	return -1;
}

__int64 powmod(__int64 a, __int64 k, __int64 n)
{
	__int64 b = 1;
	while (k) {
		if (k%2 == 0) {
			k /= 2;
			a = (a*a)%n;
		} else {
			k--;
			b = (b*a)%n;
		}
	}
	return b;
}

// Первообразный (примитивный) корень для простого числа $p степени $period

int find_primitive_root(int p, int period) {
	int i, val;
	int found_count;
	int proot;
	int pow_mod;
	int *checkarr;
	if ((p-1)%period != 0) {
		return -1;
	}
	checkarr = (int *)calloc(p, sizeof(int));
	found_count = (p-1)/period;
	for (proot = 1; proot < p; proot++) {
		pow_mod = powmod(proot, period, p);
		if (pow_mod == 1) {
			checkarr[1] = 1;
			for (i = 2; i < p; i++)
				checkarr[i] = 0;
			
			// Заполняем
			val = proot;
			for (i = 1; i < (p-1); i++) {
				checkarr[val]++;
				val = (val*proot)%p;
			}						
			
			// Проверяем
			for (i = 1; i < p; i++) {
				if (checkarr[i] != 0 && checkarr[i] != found_count) {
					break;
				}
			}
			if (i == p) {
				free(checkarr);
				return proot;
			}
		}
	}

	free(checkarr);
	return -1;
}

/*
	Теоретико числовое дискретное преобразование Фурье
*/

int dft_discrete
	(int len,
	int *x,
	int gamma,
	int mod)
{
	int i, j, pm;
	int *tmp;
	tmp = (int *)calloc(len, sizeof(int));

	for (i = 0; i < len; i++) {
		tmp[i] = 0;
		for (j = 0; j < len; j++) {
			pm = powmod(gamma, i*j, mod);
			tmp[i] += (x[j]*pm)%mod;
		}
		tmp[i] = tmp[i]%mod;
	}
	for (i = 0; i < len; i++) {
		x[i] = tmp[i];
	}

	free(tmp);
	return 0;
}

/* calculates a * *x + b * *y = gcd(a, b) = *d */
void extended_euclid(int a, int b, int *x, int *y, int *d) {
  int q, r, x1, x2, y1, y2;
  if (b == 0) {
    *d = a, *x = 1, *y = 0;
    return;
  }
  x2 = 1, x1 = 0, y2 = 0, y1 = 1;
  while (b > 0) {
    q = a / b, r = a - q * b;
    *x = x2 - q * x1, *y = y2 - q * y1;
    a = b, b = r;
    x2 = x1, x1 = *x, y2 = y1, y1 = *y;
  }
  *d = a, *x = x2, *y = y2;
}

/* computes the inverse of a modulo n */
int inverse(int a, int n) {
  int d, x, y;
  extended_euclid(a, n, &x, &y, &d);
  if (d == 1) return x;
  return 0;
}

/*
	Инверсное теоретико числовое дискретное преобразование Фурье
*/

int inverse_dft_discrete
	(int len,
	int *x,
	int gamma,
	int mod)
{
	int i, j, pm, step;
	int lenInv;
	int *tmp;
	tmp = (int *)calloc(len, sizeof(int));

	lenInv = inverse(len, mod);
	while(lenInv < 0)
		lenInv += mod;
	if (lenInv == 0) {
		printf("Error!\n");
		return 0;
	}

	for (i = 0; i < len; i++) {
		tmp[i] = 0;
		for (j = 0; j < len; j++) {
			step = -i*j;
			while(step < 0)
				step += len;
			pm = powmod(gamma, step, mod);
			tmp[i] += (x[j]*pm)%mod;
		}
		tmp[i] = (lenInv*tmp[i])%mod;
	}
	for (i = 0; i < len; i++) {
		x[i] = tmp[i];
	}

	free(tmp);
	return 0;
}

/*
	Вычисляет линейную свертку двух последовательностей (примерно одинаковой длины) 
	используя теоретико числовое ДПФ
	Длина выхода N = (len1+len2-1)
*/

int calc_simple_conv_modular (
	int period,
	int mod,
	int *in1, 
	int len1, 
	int *in2, 
	int len2, 
	int *out)
{
	int i, N;
	int gamma;
	int *x;
	int *y;
	int *r;
	
	if (period < len1+len2-1) {
		printf("Incorrect modules set!\n");
		exit(0);
	}

	// N - длина результата линейной свертки
	N = period;
	gamma = find_primitive_root(mod, period);

	x = (int *)calloc(N, sizeof(int));
	y = (int *)calloc(N, sizeof(int));
	r = (int *)calloc(N, sizeof(int));

	for (i = 0; i < len1; i++) {
		x[i] = in1[i];
	}
	for (i = len1; i < N; i++) {
		x[i] = 0;
	}
	for (i = 0; i < len2; i++) {
		y[i] = in2[i];
	}
	for (i = len2; i < N; i++) {
		y[i] = 0;
	}

	dft_discrete(N, x, gamma, mod);
	dft_discrete(N, y, gamma, mod);

	for (i = 0; i < N; i++)
	{
		r[i] = (x[i]*y[i])%mod;
	}

	inverse_dft_discrete(N, r, gamma, mod);

	for (i = 0; i < len1+len2-1; i++) {
		out[i] = r[i];
	}
	for (i = len1+len2-1; i < N; i++) {
		if (r[i] != 0) {
			printf("Some error here!\n");
		}
	}

	free(x);
	free(y);
	free(r);
	return 0;
}

/*
	Простейший подход - действуем по определению свертки
*/

int calc_fir_filter_deafult(
	int *vec1,
	int len1,
	int *vec2,
	int len2,
	int *out) 
{
	int i, j;
	for (i = 0; i < len1+len2-1; i++) {
		out[i] = 0;
		for (j = (i-len2+1 > 0)?(i-len2+1):0; j <= i; j++) {
			if (j >= len1)
				break;
			out[i] += vec1[j]*vec2[i-j];
		}
	}
	return 0;
}

/*
	Обратное преобразование
*/

int backward_convert(int len, int *mod, int *vals) 
{
	int i;
	__int64 maxvalue = 1;
	__int64 backward[128];
	__int64 res = 0;
	__int64 a, b;

	for (i = 0; i < len; i++) {
		maxvalue *= mod[i];
	}

	for (i = 0; i < len; i++)
	{
		backward[i] = maxvalue/mod[i];
		backward[i] *= powmod(backward[i], mod[i]-2, mod[i]);
	}
	
	for (i = 0; i < len; i++)
	{
		a = (__int64)vals[i];
		b = backward[i];
		res += (a*b)%maxvalue;
	}
	return (int)(res%maxvalue);
}

int find_modular_basis(int period, int *mod, int &modLen, int maxVal) 
{
	int i, k;
	int max_basis_value = 1;
	int proot;
	int val;
		
	modLen = 0;
	for (k = 1; k < 100*period; k++)
	{
		val = k*period + 1;
		if (prime[val] == 1) {
			proot = find_primitive_root(val, period);
			if (proot == -1) {
				printf("Невозможно найти примитивный корень степени %d для простого числа %d!\n", period, val);
				exit(0);
			}
			max_basis_value *= val;
			mod[modLen] = val;
			modLen++;
			if (max_basis_value > maxVal)
				break;
		}
	}
	if (k == period)
		return 0;

	// Пробуем исключить самый маленький модуль
	if (max_basis_value/mod[0] > maxVal) {
		for (i = 0; i < modLen - 1; i++)
			mod[i] = mod[i+1];
		modLen--;
	}

	return 1;
}

/*
	Свертка через модулярные каналы
	Простой подход
*/

int calc_fir_filter_modular_simple(
	int dataBit,
	int *vec1,
	int len1,
	int *vec2,
	int len2,
	int *out) 
{
	int mod[MAX_MOD_NUM];
	int vals[MAX_MOD_NUM];
	int modLen;
	int i, j, k, N;
	int *x[MAX_MOD_NUM];
	int *y[MAX_MOD_NUM];
	int *r[MAX_MOD_NUM];
	int *rTmp[MAX_MOD_NUM];
	int maxValue;
	int maxBasisValue;
	int period;

	if (len1 > len2) {
		printf("Error: len1 > len2\n");
		exit(0);
	}

	// Заполняем набор простых чисел
	
	for (i = 1; i < MAXPRIME; i++)
		prime[i] = 1;

	for (i = 2; i < MAXPRIME; i++)
		for (j = 2; i*j < MAXPRIME; j++)
			prime[i*j] = 0;

	pnum = 0;
	for (i = 1; i < MAXPRIME; i++)
		if (prime[i] == 1) {
			plist[pnum] = i;
			pnum++;
		}

	printf("Filter: ");
	for (i = 0; i < len1; i++) {
		printf("%d ", vec1[i]);
	}
	printf("\n");

	// Максимальное значение для результата свертки
	maxValue = 0;
	for (i = 0; i < len1; i++) {
		maxValue += vec1[i];
	}
	maxValue *= (int)(pow(2.0, dataBit) - 0.5);
	printf("Maximum possible value for result: %d\n", maxValue);

	// Рассчитываем модульный базис
	i = find_bit(2*len1-1);
	period = 1 << i;
	while (!find_modular_basis(period, mod, modLen, maxValue)) {
		period *= 2;
	}

	printf("Required length of convolution: %d\n", period);
	printf("Modular basis:");
	maxBasisValue = 1;
	for (i = 0; i < modLen; i++) {
		printf(" %d", mod[i]);
		maxBasisValue *= mod[i];
	}
	printf("\n");

	printf("Dynamic range of basis: %d (~ %.2lf times larger than required)\n", maxBasisValue, (1.0*maxBasisValue/maxValue));
	
	N = len1+len2-1;
	for (i = 0; i < modLen; i++) {
		x[i] = (int *)calloc(len1, sizeof(double));
		y[i] = (int *)calloc(len1, sizeof(double));
		rTmp[i] = (int *)calloc(len1+len1-1, sizeof(double));
		r[i] = (int *)calloc(N, sizeof(double));
	}

	for (i = 0; i < len1; i++) {
		for (j = 0; j < modLen; j++)
			x[j][i] = vec1[i]%mod[j];
	}
	for (i = 0; i < N; i++) {
		for (j = 0; j < modLen; j++)
			r[j][i] = 0;
	}

	for (i = 0; i < len2; i += len1)
	{
		for (j = 0; j < len1; j++) {
			for (k = 0; k < modLen; k++)
				y[k][j] = vec2[i+j]%mod[k];
		}
		for (k = 0; k < modLen; k++) {
			calc_simple_conv_modular(period, mod[k], x[k], len1, y[k], len1, rTmp[k]);
		}
		for (j = 0; j < len1+len1-1; j++) {
			for (k = 0; k < modLen; k++) {
				r[k][i+j] = (r[k][i+j] + rTmp[k][j])%mod[k];
			}
		}
	}

	for (i = 0; i < N; i++) {
		for (j = 0; j < modLen; j++) {
			printf("r[%d][%d] = %d ", j, i, r[j][i]);
			vals[j] = r[j][i];
		}
		out[i] = backward_convert(modLen, mod, vals);
		for (j = 0; j < modLen; j++) {
			if (out[i]%mod[j] != r[j][i]) {
				printf("Error backward (%d%%%d != %d its equal to %d)!\n", out[i], mod[j], r[j][i], out[i]%mod[j]);
				exit(0);
			}
		}
		printf("FINALR = %d ", out[i]);
		printf("\n");
	}

	for (i = 0; i < modLen; i++) {
		free(x[i]);
		free(y[i]);
		free(rTmp[i]);
		free(r[i]);
	}
	return 0;
}


int main ()
{
	int i, z;
	double maxVal;

	int filter[NL];
	int filterLen;
	int testArray[100*NL];
	int arrayLen;
	int dataBit;
	int res_check[10*NL];
	int res_experimental[10*NL];

	freopen("input.txt", "r", stdin);
	
	scanf("%d", &z);
	while(z--) {
		scanf("%d %d %d", &dataBit, &filterLen, &arrayLen);
		for (i = 0; i < filterLen; i++) {
			scanf("%d", &(filter[i]));
		}
		for (i = 0; i < arrayLen; i++) {
			scanf("%d", &(testArray[i]));
		}
	}
	
	maxVal = pow(2.0, 2*dataBit)*(filterLen+arrayLen);
	printf("Max value: %.0lf\n", maxVal);

	calc_fir_filter_deafult(filter, filterLen, testArray, arrayLen, res_check);
	calc_fir_filter_modular_simple(dataBit, filter, filterLen, testArray, arrayLen, res_experimental);

	// Check correctness
	for (i = 0; i < filterLen+arrayLen-1; i++) {
		if (res_check[i] != res_experimental[i]) {
			printf("Some error at position %d\n", i);
			break;
		}
	}
	if (i == filterLen+arrayLen-1) {
		printf("Congratulations! No errors were found at calcualtions!\n");
	} else {
		printf("Some errors were found at calcualtions!\n");
	}

	return 0;
}

Входные/выходные данные

Формат входных/выходных данных для FIR-фильтра

Пример (входные данные)

1
4 4 8
3 0 1 2
4 3 5 1 5 6 3 2

Пример (выходные данные)

Max value: 3072
Filter: 3 0 1 2
Maximum possible value for result: 90
Required length of convolution: 8
Modular basis: 17 41
Dynamic range of basis: 697 (~ 7.74 times larger than required)
r[0][0] = 12 r[1][0] = 12 FINALR = 12
r[0][1] = 9 r[1][1] = 9 FINALR = 9
r[0][2] = 2 r[1][2] = 19 FINALR = 19
r[0][3] = 14 r[1][3] = 14 FINALR = 14
r[0][4] = 9 r[1][4] = 26 FINALR = 26
r[0][5] = 12 r[1][5] = 29 FINALR = 29
r[0][6] = 16 r[1][6] = 16 FINALR = 16
r[0][7] = 5 r[1][7] = 22 FINALR = 22
r[0][8] = 15 r[1][8] = 15 FINALR = 15
r[0][9] = 8 r[1][9] = 8 FINALR = 8
r[0][10] = 4 r[1][10] = 4 FINALR = 4
Congratulations! No errors were found at calcualtions!