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

Материал из Модулярная арифметики
Перейти к: навигация, поиск

Содержание

Описание

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

Исходник на C

  1. #define _USE_MATH_DEFINES
  2.  
  3. #include <stdio.h>
  4. #include <math.h>
  5. #include <stdlib.h>
  6. #include <time.h>
  7. #include <string.h>
  8.  
  9. #define MAX_MOD_NUM 10
  10. #define NL 1000
  11. #define TRUE 1
  12. #define FALSE 0
  13. #define MAXPRIME 1000000
  14.  
  15. int prime[MAXPRIME];
  16. int plist[MAXPRIME];
  17. int pnum;
  18.  
  19. int find_bit(int p) {
  20. 	int i, k;
  21. 	k = 1;
  22. 	for (i = 0; i < 32; i++) {
  23. 		if (k > p)
  24. 			return i;
  25. 		k *= 2;
  26. 	}
  27. 	return -1;
  28. }
  29.  
  30. __int64 powmod(__int64 a, __int64 k, __int64 n)
  31. {
  32. 	__int64 b = 1;
  33. 	while (k) {
  34. 		if (k%2 == 0) {
  35. 			k /= 2;
  36. 			a = (a*a)%n;
  37. 		} else {
  38. 			k--;
  39. 			b = (b*a)%n;
  40. 		}
  41. 	}
  42. 	return b;
  43. }
  44.  
  45. // Первообразный (примитивный) корень для простого числа $p степени $period
  46.  
  47. int find_primitive_root(int p, int period) {
  48. 	int i, val;
  49. 	int found_count;
  50. 	int proot;
  51. 	int pow_mod;
  52. 	int *checkarr;
  53. 	if ((p-1)%period != 0) {
  54. 		return -1;
  55. 	}
  56. 	checkarr = (int *)calloc(p, sizeof(int));
  57. 	found_count = (p-1)/period;
  58. 	for (proot = 1; proot < p; proot++) {
  59. 		pow_mod = powmod(proot, period, p);
  60. 		if (pow_mod == 1) {
  61. 			checkarr[1] = 1;
  62. 			for (i = 2; i < p; i++)
  63. 				checkarr[i] = 0;
  64.  
  65. 			// Заполняем
  66. 			val = proot;
  67. 			for (i = 1; i < (p-1); i++) {
  68. 				checkarr[val]++;
  69. 				val = (val*proot)%p;
  70. 			}						
  71.  
  72. 			// Проверяем
  73. 			for (i = 1; i < p; i++) {
  74. 				if (checkarr[i] != 0 && checkarr[i] != found_count) {
  75. 					break;
  76. 				}
  77. 			}
  78. 			if (i == p) {
  79. 				free(checkarr);
  80. 				return proot;
  81. 			}
  82. 		}
  83. 	}
  84.  
  85. 	free(checkarr);
  86. 	return -1;
  87. }
  88.  
  89. /*
  90. 	Теоретико числовое дискретное преобразование Фурье
  91. */
  92.  
  93. int dft_discrete
  94. 	(int len,
  95. 	int *x,
  96. 	int gamma,
  97. 	int mod)
  98. {
  99. 	int i, j, pm;
  100. 	int *tmp;
  101. 	tmp = (int *)calloc(len, sizeof(int));
  102.  
  103. 	for (i = 0; i < len; i++) {
  104. 		tmp[i] = 0;
  105. 		for (j = 0; j < len; j++) {
  106. 			pm = powmod(gamma, i*j, mod);
  107. 			tmp[i] += (x[j]*pm)%mod;
  108. 		}
  109. 		tmp[i] = tmp[i]%mod;
  110. 	}
  111. 	for (i = 0; i < len; i++) {
  112. 		x[i] = tmp[i];
  113. 	}
  114.  
  115. 	free(tmp);
  116. 	return 0;
  117. }
  118.  
  119. /* calculates a * *x + b * *y = gcd(a, b) = *d */
  120. void extended_euclid(int a, int b, int *x, int *y, int *d) {
  121.   int q, r, x1, x2, y1, y2;
  122.   if (b == 0) {
  123.     *d = a, *x = 1, *y = 0;
  124.     return;
  125.   }
  126.   x2 = 1, x1 = 0, y2 = 0, y1 = 1;
  127.   while (b > 0) {
  128.     q = a / b, r = a - q * b;
  129.     *x = x2 - q * x1, *y = y2 - q * y1;
  130.     a = b, b = r;
  131.     x2 = x1, x1 = *x, y2 = y1, y1 = *y;
  132.   }
  133.   *d = a, *x = x2, *y = y2;
  134. }
  135.  
  136. /* computes the inverse of a modulo n */
  137. int inverse(int a, int n) {
  138.   int d, x, y;
  139.   extended_euclid(a, n, &x, &y, &d);
  140.   if (d == 1) return x;
  141.   return 0;
  142. }
  143.  
  144. /*
  145. 	Инверсное теоретико числовое дискретное преобразование Фурье
  146. */
  147.  
  148. int inverse_dft_discrete
  149. 	(int len,
  150. 	int *x,
  151. 	int gamma,
  152. 	int mod)
  153. {
  154. 	int i, j, pm, step;
  155. 	int lenInv;
  156. 	int *tmp;
  157. 	tmp = (int *)calloc(len, sizeof(int));
  158.  
  159. 	lenInv = inverse(len, mod);
  160. 	while(lenInv < 0)
  161. 		lenInv += mod;
  162. 	if (lenInv == 0) {
  163. 		printf("Error!\n");
  164. 		return 0;
  165. 	}
  166.  
  167. 	for (i = 0; i < len; i++) {
  168. 		tmp[i] = 0;
  169. 		for (j = 0; j < len; j++) {
  170. 			step = -i*j;
  171. 			while(step < 0)
  172. 				step += len;
  173. 			pm = powmod(gamma, step, mod);
  174. 			tmp[i] += (x[j]*pm)%mod;
  175. 		}
  176. 		tmp[i] = (lenInv*tmp[i])%mod;
  177. 	}
  178. 	for (i = 0; i < len; i++) {
  179. 		x[i] = tmp[i];
  180. 	}
  181.  
  182. 	free(tmp);
  183. 	return 0;
  184. }
  185.  
  186. /*
  187. 	Вычисляет линейную свертку двух последовательностей (примерно одинаковой длины) 
  188. 	используя теоретико числовое ДПФ
  189. 	Длина выхода N = (len1+len2-1)
  190. */
  191.  
  192. int calc_simple_conv_modular (
  193. 	int period,
  194. 	int mod,
  195. 	int *in1, 
  196. 	int len1, 
  197. 	int *in2, 
  198. 	int len2, 
  199. 	int *out)
  200. {
  201. 	int i, N;
  202. 	int gamma;
  203. 	int *x;
  204. 	int *y;
  205. 	int *r;
  206.  
  207. 	if (period < len1+len2-1) {
  208. 		printf("Incorrect modules set!\n");
  209. 		exit(0);
  210. 	}
  211.  
  212. 	// N - длина результата линейной свертки
  213. 	N = period;
  214. 	gamma = find_primitive_root(mod, period);
  215.  
  216. 	x = (int *)calloc(N, sizeof(int));
  217. 	y = (int *)calloc(N, sizeof(int));
  218. 	r = (int *)calloc(N, sizeof(int));
  219.  
  220. 	for (i = 0; i < len1; i++) {
  221. 		x[i] = in1[i];
  222. 	}
  223. 	for (i = len1; i < N; i++) {
  224. 		x[i] = 0;
  225. 	}
  226. 	for (i = 0; i < len2; i++) {
  227. 		y[i] = in2[i];
  228. 	}
  229. 	for (i = len2; i < N; i++) {
  230. 		y[i] = 0;
  231. 	}
  232.  
  233. 	dft_discrete(N, x, gamma, mod);
  234. 	dft_discrete(N, y, gamma, mod);
  235.  
  236. 	for (i = 0; i < N; i++)
  237. 	{
  238. 		r[i] = (x[i]*y[i])%mod;
  239. 	}
  240.  
  241. 	inverse_dft_discrete(N, r, gamma, mod);
  242.  
  243. 	for (i = 0; i < len1+len2-1; i++) {
  244. 		out[i] = r[i];
  245. 	}
  246. 	for (i = len1+len2-1; i < N; i++) {
  247. 		if (r[i] != 0) {
  248. 			printf("Some error here!\n");
  249. 		}
  250. 	}
  251.  
  252. 	free(x);
  253. 	free(y);
  254. 	free(r);
  255. 	return 0;
  256. }
  257.  
  258. /*
  259. 	Простейший подход - действуем по определению свертки
  260. */
  261.  
  262. int calc_fir_filter_deafult(
  263. 	int *vec1,
  264. 	int len1,
  265. 	int *vec2,
  266. 	int len2,
  267. 	int *out) 
  268. {
  269. 	int i, j;
  270. 	for (i = 0; i < len1+len2-1; i++) {
  271. 		out[i] = 0;
  272. 		for (j = (i-len2+1 > 0)?(i-len2+1):0; j <= i; j++) {
  273. 			if (j >= len1)
  274. 				break;
  275. 			out[i] += vec1[j]*vec2[i-j];
  276. 		}
  277. 	}
  278. 	return 0;
  279. }
  280.  
  281. /*
  282. 	Обратное преобразование
  283. */
  284.  
  285. int backward_convert(int len, int *mod, int *vals) 
  286. {
  287. 	int i;
  288. 	__int64 maxvalue = 1;
  289. 	__int64 backward[128];
  290. 	__int64 res = 0;
  291. 	__int64 a, b;
  292.  
  293. 	for (i = 0; i < len; i++) {
  294. 		maxvalue *= mod[i];
  295. 	}
  296.  
  297. 	for (i = 0; i < len; i++)
  298. 	{
  299. 		backward[i] = maxvalue/mod[i];
  300. 		backward[i] *= powmod(backward[i], mod[i]-2, mod[i]);
  301. 	}
  302.  
  303. 	for (i = 0; i < len; i++)
  304. 	{
  305. 		a = (__int64)vals[i];
  306. 		b = backward[i];
  307. 		res += (a*b)%maxvalue;
  308. 	}
  309. 	return (int)(res%maxvalue);
  310. }
  311.  
  312. int find_modular_basis(int period, int *mod, int &modLen, int maxVal) 
  313. {
  314. 	int i, k;
  315. 	int max_basis_value = 1;
  316. 	int proot;
  317. 	int val;
  318.  
  319. 	modLen = 0;
  320. 	for (k = 1; k < 100*period; k++)
  321. 	{
  322. 		val = k*period + 1;
  323. 		if (prime[val] == 1) {
  324. 			proot = find_primitive_root(val, period);
  325. 			if (proot == -1) {
  326. 				printf("Невозможно найти примитивный корень степени %d для простого числа %d!\n", period, val);
  327. 				exit(0);
  328. 			}
  329. 			max_basis_value *= val;
  330. 			mod[modLen] = val;
  331. 			modLen++;
  332. 			if (max_basis_value > maxVal)
  333. 				break;
  334. 		}
  335. 	}
  336. 	if (k == period)
  337. 		return 0;
  338.  
  339. 	// Пробуем исключить самый маленький модуль
  340. 	if (max_basis_value/mod[0] > maxVal) {
  341. 		for (i = 0; i < modLen - 1; i++)
  342. 			mod[i] = mod[i+1];
  343. 		modLen--;
  344. 	}
  345.  
  346. 	return 1;
  347. }
  348.  
  349. /*
  350. 	Свертка через модулярные каналы
  351. 	Простой подход
  352. */
  353.  
  354. int calc_fir_filter_modular_simple(
  355. 	int dataBit,
  356. 	int *vec1,
  357. 	int len1,
  358. 	int *vec2,
  359. 	int len2,
  360. 	int *out) 
  361. {
  362. 	int mod[MAX_MOD_NUM];
  363. 	int vals[MAX_MOD_NUM];
  364. 	int modLen;
  365. 	int i, j, k, N;
  366. 	int *x[MAX_MOD_NUM];
  367. 	int *y[MAX_MOD_NUM];
  368. 	int *r[MAX_MOD_NUM];
  369. 	int *rTmp[MAX_MOD_NUM];
  370. 	int maxValue;
  371. 	int maxBasisValue;
  372. 	int period;
  373.  
  374. 	if (len1 > len2) {
  375. 		printf("Error: len1 > len2\n");
  376. 		exit(0);
  377. 	}
  378.  
  379. 	// Заполняем набор простых чисел
  380.  
  381. 	for (i = 1; i < MAXPRIME; i++)
  382. 		prime[i] = 1;
  383.  
  384. 	for (i = 2; i < MAXPRIME; i++)
  385. 		for (j = 2; i*j < MAXPRIME; j++)
  386. 			prime[i*j] = 0;
  387.  
  388. 	pnum = 0;
  389. 	for (i = 1; i < MAXPRIME; i++)
  390. 		if (prime[i] == 1) {
  391. 			plist[pnum] = i;
  392. 			pnum++;
  393. 		}
  394.  
  395. 	printf("Filter: ");
  396. 	for (i = 0; i < len1; i++) {
  397. 		printf("%d ", vec1[i]);
  398. 	}
  399. 	printf("\n");
  400.  
  401. 	// Максимальное значение для результата свертки
  402. 	maxValue = 0;
  403. 	for (i = 0; i < len1; i++) {
  404. 		maxValue += vec1[i];
  405. 	}
  406. 	maxValue *= (int)(pow(2.0, dataBit) - 0.5);
  407. 	printf("Maximum possible value for result: %d\n", maxValue);
  408.  
  409. 	// Рассчитываем модульный базис
  410. 	i = find_bit(2*len1-1);
  411. 	period = 1 << i;
  412. 	while (!find_modular_basis(period, mod, modLen, maxValue)) {
  413. 		period *= 2;
  414. 	}
  415.  
  416. 	printf("Required length of convolution: %d\n", period);
  417. 	printf("Modular basis:");
  418. 	maxBasisValue = 1;
  419. 	for (i = 0; i < modLen; i++) {
  420. 		printf(" %d", mod[i]);
  421. 		maxBasisValue *= mod[i];
  422. 	}
  423. 	printf("\n");
  424.  
  425. 	printf("Dynamic range of basis: %d (~ %.2lf times larger than required)\n", maxBasisValue, (1.0*maxBasisValue/maxValue));
  426.  
  427. 	N = len1+len2-1;
  428. 	for (i = 0; i < modLen; i++) {
  429. 		x[i] = (int *)calloc(len1, sizeof(double));
  430. 		y[i] = (int *)calloc(len1, sizeof(double));
  431. 		rTmp[i] = (int *)calloc(len1+len1-1, sizeof(double));
  432. 		r[i] = (int *)calloc(N, sizeof(double));
  433. 	}
  434.  
  435. 	for (i = 0; i < len1; i++) {
  436. 		for (j = 0; j < modLen; j++)
  437. 			x[j][i] = vec1[i]%mod[j];
  438. 	}
  439. 	for (i = 0; i < N; i++) {
  440. 		for (j = 0; j < modLen; j++)
  441. 			r[j][i] = 0;
  442. 	}
  443.  
  444. 	for (i = 0; i < len2; i += len1)
  445. 	{
  446. 		for (j = 0; j < len1; j++) {
  447. 			for (k = 0; k < modLen; k++)
  448. 				y[k][j] = vec2[i+j]%mod[k];
  449. 		}
  450. 		for (k = 0; k < modLen; k++) {
  451. 			calc_simple_conv_modular(period, mod[k], x[k], len1, y[k], len1, rTmp[k]);
  452. 		}
  453. 		for (j = 0; j < len1+len1-1; j++) {
  454. 			for (k = 0; k < modLen; k++) {
  455. 				r[k][i+j] = (r[k][i+j] + rTmp[k][j])%mod[k];
  456. 			}
  457. 		}
  458. 	}
  459.  
  460. 	for (i = 0; i < N; i++) {
  461. 		for (j = 0; j < modLen; j++) {
  462. 			printf("r[%d][%d] = %d ", j, i, r[j][i]);
  463. 			vals[j] = r[j][i];
  464. 		}
  465. 		out[i] = backward_convert(modLen, mod, vals);
  466. 		for (j = 0; j < modLen; j++) {
  467. 			if (out[i]%mod[j] != r[j][i]) {
  468. 				printf("Error backward (%d%%%d != %d its equal to %d)!\n", out[i], mod[j], r[j][i], out[i]%mod[j]);
  469. 				exit(0);
  470. 			}
  471. 		}
  472. 		printf("FINALR = %d ", out[i]);
  473. 		printf("\n");
  474. 	}
  475.  
  476. 	for (i = 0; i < modLen; i++) {
  477. 		free(x[i]);
  478. 		free(y[i]);
  479. 		free(rTmp[i]);
  480. 		free(r[i]);
  481. 	}
  482. 	return 0;
  483. }
  484.  
  485.  
  486. int main ()
  487. {
  488. 	int i, z;
  489. 	double maxVal;
  490.  
  491. 	int filter[NL];
  492. 	int filterLen;
  493. 	int testArray[100*NL];
  494. 	int arrayLen;
  495. 	int dataBit;
  496. 	int res_check[10*NL];
  497. 	int res_experimental[10*NL];
  498.  
  499. 	freopen("input.txt", "r", stdin);
  500.  
  501. 	scanf("%d", &z);
  502. 	while(z--) {
  503. 		scanf("%d %d %d", &dataBit, &filterLen, &arrayLen);
  504. 		for (i = 0; i < filterLen; i++) {
  505. 			scanf("%d", &(filter[i]));
  506. 		}
  507. 		for (i = 0; i < arrayLen; i++) {
  508. 			scanf("%d", &(testArray[i]));
  509. 		}
  510. 	}
  511.  
  512. 	maxVal = pow(2.0, 2*dataBit)*(filterLen+arrayLen);
  513. 	printf("Max value: %.0lf\n", maxVal);
  514.  
  515. 	calc_fir_filter_deafult(filter, filterLen, testArray, arrayLen, res_check);
  516. 	calc_fir_filter_modular_simple(dataBit, filter, filterLen, testArray, arrayLen, res_experimental);
  517.  
  518. 	// Check correctness
  519. 	for (i = 0; i < filterLen+arrayLen-1; i++) {
  520. 		if (res_check[i] != res_experimental[i]) {
  521. 			printf("Some error at position %d\n", i);
  522. 			break;
  523. 		}
  524. 	}
  525. 	if (i == filterLen+arrayLen-1) {
  526. 		printf("Congratulations! No errors were found at calcualtions!\n");
  527. 	} else {
  528. 		printf("Some errors were found at calcualtions!\n");
  529. 	}
  530.  
  531. 	return 0;
  532. }

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

Формат входных/выходных данных для 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!

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

Варианты
Действия
Навигация