How to resolve the algorithm Deconvolution/2D+ step by step in the C programming language

Published on 7 June 2024 03:52 AM
#C

How to resolve the algorithm Deconvolution/2D+ step by step in the C programming language

Table of Contents

Problem Statement

This task is a straightforward generalization of Deconvolution/1D to higher dimensions. For example, the one dimensional case would be applicable to audio signals, whereas two dimensions would pertain to images. Define the discrete convolution in

d

{\displaystyle {\mathit {d}}}

dimensions of two functions taking

d

{\displaystyle {\mathit {d}}}

-tuples of integers to real numbers as the function also taking

d

{\displaystyle {\mathit {d}}}

-tuples of integers to reals and satisfying for all

d

{\displaystyle {\mathit {d}}}

-tuples of integers

(

n

0

, … ,

n

d − 1

) ∈

Z

d

{\displaystyle (n_{0},\dots ,n_{d-1})\in \mathbb {Z} ^{d}}

. Assume

F

{\displaystyle {\mathit {F}}}

and

H

{\displaystyle {\mathit {H}}}

(and therefore

G

{\displaystyle {\mathit {G}}}

) are non-zero over only a finite domain bounded by the origin, hence possible to represent as finite multi-dimensional arrays or nested lists

f

{\displaystyle {\mathit {f}}}

,

h

{\displaystyle {\mathit {h}}}

, and

g

{\displaystyle {\mathit {g}}}

. For this task, implement a function (or method, procedure, subroutine, etc.) deconv to perform deconvolution (i.e., the inverse of convolution) by solving for

h

{\displaystyle {\mathit {h}}}

given

f

{\displaystyle {\mathit {f}}}

and

g

{\displaystyle {\mathit {g}}}

. (See Deconvolution/1D for details.) dimension 1: dimension 2: dimension 3:

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Deconvolution/2D+ step by step in the C programming language

This C program implements a deconvolution algorithm. Here's a brief overview of the code:

  1. Constants and Data Structures:

    • It defines PI as atan2(1, 1) * 4 to use in complex number calculations.
    • It also defines a complex number type cplx using the complex.h library.

  2. FFT and IFFT Functions:

    • _fft is a recursive function that implements the Fast Fourier Transform (FFT) algorithm.
    • fft is a wrapper function that calls _fft with the appropriate parameters to perform FFT on a complex array.

  3. Zero Padding Function:

    • pad_two adds zero values to an array to increase its length to the nearest power of two. This is necessary for efficient FFT calculations.

  4. Deconvolution Function:

    • deconv performs deconvolution of two signals g and f in the frequency domain. It applies FFT to both signals, divides them (element-wise), applies IFFT to the result, and returns the deconvolved signal h.

  5. Utility Functions:

    • unpack2, pack2, unpack3, and pack3 are functions to convert between 1D arrays and 2D/3D arrays. They repack data for use with the deconvolution functions.

  6. Main Function:

    • The main function initializes sample data for 2D and 3D deconvolution scenarios.
    • It calls the deconv3 function to perform 3D deconvolution of g and f and h, and prints the results.
    • The commented-out part of the main function shows an example of 2D deconvolution, but it is not fully executed in this code.

In summary, this program allows you to perform multidimensional deconvolution by converting signals to the frequency domain, dividing them, and transforming back to the time domain. The results are stored in the h2 variable, and the corresponding f2 variable stores the result of the deconvolution of g and h.

Source code in the c programming language

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
 
double PI;
typedef double complex cplx;
 
void _fft(cplx buf[], cplx out[], int n, int step)
{
	if (step < n) {
		_fft(out, buf, n, step * 2);
		_fft(out + step, buf + step, n, step * 2);
 
		for (int i = 0; i < n; i += 2 * step) {
			cplx t = cexp(-I * PI * i / n) * out[i + step];
			buf[i / 2]     = out[i] + t;
			buf[(i + n)/2] = out[i] - t;
		}
	}
}
 
void fft(cplx buf[], int n)
{
	cplx out[n];
	for (int i = 0; i < n; i++) out[i] = buf[i];
	_fft(buf, out, n, 1);
}
 
/* pad array length to power of two */
cplx *pad_two(double g[], int len, int *ns)
{
	int n = 1;
	if (*ns) n = *ns;
	else while (n < len) n *= 2;
 
	cplx *buf = calloc(sizeof(cplx), n);
	for (int i = 0; i < len; i++) buf[i] = g[i];
	*ns = n;
	return buf;
}
 
void deconv(double g[], int lg, double f[], int lf, double out[], int row_len) {
	int ns = 0;
	cplx *g2 = pad_two(g, lg, &ns);
	cplx *f2 = pad_two(f, lf, &ns);
 
	fft(g2, ns);
	fft(f2, ns);
 
	cplx h[ns];
	for (int i = 0; i < ns; i++) h[i] = g2[i] / f2[i];
	fft(h, ns);

	for (int i = 0; i < ns; i++) {
		if (cabs(creal(h[i])) < 1e-10) 
			h[i] = 0;
	}

	for (int i = 0; i > lf - lg - row_len; i--)
		out[-i] = h[(i + ns) % ns]/32;
	free(g2);
	free(f2);
}

double* unpack2(void *m, int rows, int len, int to_len)
{
	double *buf = calloc(sizeof(double), rows * to_len);
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < len; j++)
			buf[i * to_len + j] = ((double(*)[len])m)[i][j];
	return buf;
}

void pack2(double * buf, int rows, int from_len, int to_len, void *out)
{
	for (int i = 0; i < rows; i++)
		for (int j = 0; j < to_len; j++)
			((double(*)[to_len])out)[i][j] = buf[i * from_len + j] / 4;
}

void deconv2(void *g, int row_g, int col_g, void *f, int row_f, int col_f, void *out) {
	double *g2 = unpack2(g, row_g, col_g, col_g);
	double *f2 = unpack2(f, row_f, col_f, col_g);

	double ff[(row_g - row_f + 1) * col_g];
	deconv(g2, row_g * col_g, f2, row_f * col_g, ff, col_g);
	pack2(ff, row_g - row_f + 1, col_g, col_g - col_f + 1, out);

	free(g2);
	free(f2);
}

double* unpack3(void *m, int x, int y, int z, int to_y, int to_z)
{
	double *buf = calloc(sizeof(double), x * to_y * to_z);
	for (int i = 0; i < x; i++)
		for (int j = 0; j < y; j++) {
			for (int k = 0; k < z; k++)
				buf[(i * to_y + j) * to_z + k] =
					((double(*)[y][z])m)[i][j][k];
		}
	return buf;
}

void pack3(double * buf, int x, int y, int z, int to_y, int to_z, void *out)
{
	for (int i = 0; i < x; i++)
		for (int j = 0; j < to_y; j++)
			for (int k = 0; k < to_z; k++)
				((double(*)[to_y][to_z])out)[i][j][k] =
					buf[(i * y + j) * z + k] / 4;
}

void deconv3(void *g, int gx, int gy, int gz, void *f, int fx, int fy, int fz, void *out) {
	double *g2 = unpack3(g, gx, gy, gz, gy, gz);
	double *f2 = unpack3(f, fx, fy, fz, gy, gz);

	double ff[(gx - fx + 1) * gy * gz];
	deconv(g2, gx * gy * gz, f2, fx * gy * gz, ff, gy * gz);
	pack3(ff, gx - fx + 1, gy, gz, gy - fy + 1, gz - fz + 1, out);

	free(g2);
	free(f2);
}

int main()
{
	PI = atan2(1,1) * 4;
	double h[2][3][4] = {
		{{-6, -8, -5,  9}, {-7, 9, -6, -8}, { 2, -7,  9,  8}},
		{{ 7,  4,  4, -6}, { 9, 9,  4, -4}, {-3,  7, -2, -3}}
	};
	int hx = 2, hy = 3, hz = 4;
	double f[3][2][3] = {	{{-9,  5, -8}, { 3,  5,  1}},
				{{-1, -7,  2}, {-5, -6,  6}},
				{{ 8,  5,  8}, {-2, -6, -4}} };
	int fx = 3, fy = 2, fz = 3;
	double g[4][4][6] = {
		{	{ 54,  42,  53, -42,  85, -72}, { 45,-170,  94, -36,  48,  73},
			{-39,  65,-112, -16, -78, -72}, {  6, -11,  -6,  62,  49,   8} },
		{ 	{-57,  49, -23,   52, -135,  66},{-23, 127, -58,   -5, -118,  64},
			{ 87, -16,  121,  23,  -41, -12},{-19,  29,   35,-148,  -11,  45} },
		{	{-55, -147, -146, -31,  55,  60},{-88,  -45,  -28,  46, -26,-144},
			{-12, -107,  -34, 150, 249,  66},{ 11,  -15,  -34,  27, -78, -50} },
		{	{ 56,  67, 108,   4,  2,-48},{ 58,  67,  89,  32, 32, -8},
			{-42, -31,-103, -30,-23, -8},{  6,   4, -26, -10, 26, 12}
		}
	};
	int gx = 4, gy = 4, gz = 6;

	double h2[gx - fx + 1][gy - fy + 1][gz - fz + 1];
	deconv3(g, gx, gy, gz, f, fx, fy, fz, h2);
	printf("deconv3(g, f):\n");
	for (int i = 0; i < gx - fx + 1; i++) {
		for (int j = 0; j < gy - fy + 1; j++) {
			for (int k = 0; k < gz - fz + 1; k++)
				printf("%g ", h2[i][j][k]);
			printf("\n");
		}
		if (i < gx - fx) printf("\n");
	}

	double f2[gx - hx + 1][gy - hy + 1][gz - hz + 1];
	deconv3(g, gx, gy, gz, h, hx, hy, hz, f2);
	printf("\ndeconv3(g, h):\n");
	for (int i = 0; i < gx - hx + 1; i++) {
		for (int j = 0; j < gy - hy + 1; j++) {
			for (int k = 0; k < gz - hz + 1; k++)
				printf("%g ", f2[i][j][k]);
			printf("\n");
		}
		if (i < gx - hx) printf("\n");
	}
}

/* two-D case; since task doesn't require showing it, it's commented out */
/*
int main()
{
	PI = atan2(1,1) * 4;
	double h[][6] = { 	{-8, 1, -7, -2, -9, 4},
				{4, 5, -5, 2, 7, -1},
				{-6, -3, -3, -6, 9, 5} };
	int hr = 3, hc = 6;

	double f[][5] = {	{-5, 2, -2, -6, -7},
				{9, 7, -6, 5, -7},
				{1, -1, 9, 2, -7},
				{5, 9, -9, 2, -5},
				{-8, 5, -2, 8, 5} };
	int fr = 5, fc = 5;
	double g[][10] = {
			{40, -21, 53, 42, 105, 1, 87, 60, 39, -28},
			{-92, -64, 19, -167, -71, -47, 128, -109, 40, -21},
			{58, 85, -93, 37, 101, -14, 5, 37, -76, -56},
			{-90, -135, 60, -125, 68, 53, 223, 4, -36, -48},
			{78, 16, 7, -199, 156, -162, 29, 28, -103, -10},
			{-62, -89, 69, -61, 66, 193, -61, 71, -8, -30},
			{48, -6, 21, -9, -150, -22, -56, 32, 85, 25}	};
	int gr = 7, gc = 10;

	double h2[gr - fr + 1][gc - fc + 1];
	deconv2(g, gr, gc, f, fr, fc, h2);
	for (int i = 0; i < gr - fr + 1; i++) {
		for (int j = 0; j < gc - fc + 1; j++)
			printf(" %g", h2[i][j]);
		printf("\n");
	}

	double f2[gr - hr + 1][gc - hc + 1];
	deconv2(g, gr, gc, h, hr, hc, f2);
	for (int i = 0; i < gr - hr + 1; i++) {
		for (int j = 0; j < gc - hc + 1; j++)
			printf(" %g", f2[i][j]);
		printf("\n");
	}
}*/


deconv3(g, f):
-6 -8 -5 9 
-7 9 -6 -8 
2 -7 9 8 

7 4 4 -6 
9 9 4 -4 
-3 7 -2 -3 

deconv3(g, h):
-9 5 -8 
3 5 1 

-1 -7 2 
-5 -6 6 

8 5 8 
-2 -6 -4

You may also check:How to resolve the algorithm Factorial step by step in the Standard ML programming language
You may also check:How to resolve the algorithm Monads/List monad step by step in the Phix programming language
You may also check:How to resolve the algorithm Break OO privacy step by step in the OCaml programming language
You may also check:How to resolve the algorithm Higher-order functions step by step in the Python programming language
You may also check:How to resolve the algorithm Find limit of recursion step by step in the Erlang programming language