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

Содержание

1 Описание
2 Исходник на C
3 Входные/выходные данные
- 3.1 Пример (входные данные)
4 Пример (выходные данные)

Описание

Расчет 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!

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

Содержание

Описание

Исходник на C

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

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

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

Навигация

Персональные инструменты

Пространства имён

Варианты

Просмотры

Действия

Поиск

Навигация

Инструменты