.gitignore

.*
!.gitignore
*.o
test_dft
test_fft

README.md

Digital Signal Processing Study

• DFT
• FFT
• etc.

dft.c

/*
 * DFT
 * MIT License Copyright(c) 2019 Hiroshi Shimamoto
 */

#include <math.h>
#include "dft.h"

void dft(double complex *tm, double complex *fr, int nr)
{
	for (int i = 0; i < nr; i++) {
		double complex f = CMPLX(0.0, 0.0);
		for (int k = 0; k < nr; k++) {
			double th = 2 * M_PI * k * i / nr;
			f += tm[k] * cexp(CMPLX(0.0, -th));
		}
		fr[i] = f;
	}
}

void idft(double complex *fr, double complex *tm, int nr)
{
	for (int i = 0; i < nr; i++) {
		double complex t = CMPLX(0.0, 0.0);
		for (int k = 0; k < nr; k++) {
			double th = 2 * M_PI * k * i / nr;
			t += fr[k] * cexp(CMPLX(0.0, th));
		}
		tm[i] = t / nr;
	}
}

dft.h

#ifndef DFT_H
#define DFT_H

#include <complex.h>

void dft(double complex *tm, double complex *fr, int nr);
void idft(double complex *fr, double complex *tm, int nr);

#endif // DFT_H

fft.c

/*
 * FFT
 * MIT License Copyright(c) 2019 Hiroshi Shimamoto
 */

#include <math.h>
#include "fft.h"

static inline int fft_index(int *idx, int nr)
{
	for (int i = 0; i < nr; i++)
		idx[i] = i;
	int o, s, n;
	for (o = 0, s = 1, n = nr >> 1; s < nr; o++, s <<= 1, n >>= 1) {
		for (int i = 0; i < s; i++)
			idx[i + s] = idx[i] + n;
	}
	if (s != nr)
		return -1;
	return o;
}

void fft(double complex *tm, double complex *fr, int nr)
{
	int idx[nr];
	int o = fft_index(idx, nr);

	if (o == -1)
		return; // TODO: panic here

	for (int i = 0; i < nr; i++)
		fr[i] = tm[idx[i]];
	for (int l = 0, p = 2; l < o; l++, p <<= 1) {
		int h = p >> 1;
		complex w = cexp(CMPLX(0.0, -2 * M_PI / p));
		complex ws = CMPLX(1.0, 0.0);
		for (int k = 0; k < h; k++) {
			for (int i = 0; i < nr; i += p) {
				complex a = fr[i + k], b = fr[i + k + h];
				complex wf = ws * b;
				fr[i + k] = a + wf;
				fr[i + k + h] = a - wf;
			}
			ws *= w;
		}
	}
}

void ifft(double complex *fr, double complex *tm, int nr)
{
	int idx[nr];
	int o = fft_index(idx, nr);

	if (o == -1)
		return; // TODO: panic here

	for (int i = 0; i < nr; i++)
		tm[i] = fr[idx[i]];
	for (int l = 0, p = 2; l < o; l++, p <<= 1) {
		int h = p >> 1;
		complex w = cexp(CMPLX(0.0, 2 * M_PI / p));
		complex ws = CMPLX(1.0, 0.0);
		for (int k = 0; k < h; k++) {
			for (int i = 0; i < nr; i += p) {
				complex a = tm[i + k], b = tm[i + k + h];
				complex wf = ws * b;
				tm[i + k] = a + wf;
				tm[i + k + h] = a - wf;
			}
			ws *= w;
		}
	}
	for (int i = 0; i < nr; i++)
		tm[i] /= nr;
}

fft.h

#ifndef FFT_H
#define FFT_H

#include <complex.h>

void fft(double complex *tm, double complex *fr, int nr);
void ifft(double complex *fr, double complex *tm, int nr);

#endif //FDFT_H

Makefile

# Digital Signal Processing
# MIT License Copyright(c) 2019 Hiroshi Shimamoto
CC = gcc
LD = gcc
CFLAGS = -g -O2
LDFLAGS = -g -O2
LIBS = -lm

test_dft_objs = test_dft.o dft.o
test_fft_objs = test_fft.o fft.o

all: test_dft test_fft

test_dft: test_dft.o dft.o
	$(LD) $(LDFLAGS) -o $@ $($(@)_objs) $(LIBS)

test_fft: test_fft.o fft.o
	$(LD) $(LDFLAGS) -o $@ $($(@)_objs) $(LIBS)

%.o: %.c
	$(CC) $(CFLAGS) -c -o $@ $<

clean:
	rm -f *.o
	rm -f test_dft

test_dft.c

/*
 * Test DFT
 * MIT License Copyright(c) 2019 Hiroshi Shimamoto
 */

#include <stdio.h>
#include <math.h>
#include "dft.h"

void dump(double complex *data, int nr)
{
	for (int i = 0; i < nr; i++)
		printf("[%d] %.f%+.fj\n", i, creal(data[i]), cimag(data[i]));
}

int main(int argc, char **argv)
{
	const int nr = 256;
	const int hz = 16;
	double complex tm[nr];

	// create sin curve
	for (int i = 0; i < nr; i++)
		tm[i] = CMPLX(10000.0 * sin(2 * M_PI * i / hz), 0.0);

	double complex fr[nr];

	printf("DFT\n");
	dft(tm, fr, nr);
	dump(fr, nr);

	double complex tm2[nr];

	printf("IDFT\n");
	idft(fr, tm2, nr);
	dump(tm2, nr);

	return 0;
}

test_fft.c

/*
 * Test FFT
 * MIT License Copyright(c) 2019 Hiroshi Shimamoto
 */

#include <stdio.h>
#include <math.h>
#include "fft.h"

void dump(double complex *data, int nr)
{
	for (int i = 0; i < nr; i++)
		printf("[%d] %.f%+.fj\n", i, creal(data[i]), cimag(data[i]));
}

int main(int argc, char **argv)
{
	const int nr = 256;
	const int hz = 16;
	double complex tm[nr];

	// create sin curve
	for (int i = 0; i < nr; i++)
		tm[i] = CMPLX(10000.0 * sin(2 * M_PI * i / hz), 0.0);

	double complex fr[nr];

	printf("FFT\n");
	fft(tm, fr, nr);
	dump(fr, nr);

	double complex tm2[nr];

	printf("IFFT\n");
	ifft(fr, tm2, nr);
	dump(tm2, nr);

	return 0;
}