Skip to content

Instantly share code, notes, and snippets.

@hshimamoto
Last active October 21, 2019 21:52
Show Gist options
  • Save hshimamoto/1b07a381a67c06f3d389d6b7482eb7a6 to your computer and use it in GitHub Desktop.
Save hshimamoto/1b07a381a67c06f3d389d6b7482eb7a6 to your computer and use it in GitHub Desktop.
.*
!.gitignore
*.o
test_dft
test_fft

Digital Signal Processing Study

  • DFT
  • FFT
  • etc.
/*
* 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;
}
}
#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
* 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;
}
#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
# 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
* 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
* 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment